text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ApplicationGatewaysApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the backend health of the specified application gateway in a resource group.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: ApplicationGatewayBackendHealth
"""
function applicationGatewaysBackendHealth(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", ApplicationGatewayBackendHealth, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}/backendhealth", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Creates or updates the specified application gateway.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: parameters::ApplicationGateway (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ApplicationGateway
"""
function applicationGatewaysCreateOrUpdate(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ApplicationGateway, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified application gateway.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function applicationGatewaysDelete(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified application gateway.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ApplicationGateway
"""
function applicationGatewaysGet(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ApplicationGateway, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all application gateways in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ApplicationGatewayListResult
"""
function applicationGatewaysList(_api::ApplicationGatewaysApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ApplicationGatewayListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all the application gateways in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ApplicationGatewayListResult
"""
function applicationGatewaysListAll(_api::ApplicationGatewaysApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ApplicationGatewayListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/applicationGateways", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all available web application firewall rule sets.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ApplicationGatewayAvailableWafRuleSetsResult
"""
function applicationGatewaysListAvailableWafRuleSets(_api::ApplicationGatewaysApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ApplicationGatewayAvailableWafRuleSetsResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/applicationGatewayAvailableWafRuleSets", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Starts the specified application gateway.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function applicationGatewaysStart(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}/start", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Stops the specified application gateway in a resource group.
Param: resourceGroupName::String (required)
Param: applicationGatewayName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function applicationGatewaysStop(_api::ApplicationGatewaysApi, resourceGroupName::String, applicationGatewayName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/applicationGateways/{applicationGatewayName}/stop", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "applicationGatewayName", applicationGatewayName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export applicationGatewaysBackendHealth, applicationGatewaysCreateOrUpdate, applicationGatewaysDelete, applicationGatewaysGet, applicationGatewaysList, applicationGatewaysListAll, applicationGatewaysListAvailableWafRuleSets, applicationGatewaysStart, applicationGatewaysStop


text:
# # [Rotating Beam with a Swept Tip](@id rotating)
#
# In this example we analyze a rotating beam with a swept tip.  The parameters for this
# example come from "Finite element solution of nonlinear intrinsic equations for curved
# composite beams" by Hodges, Shang, and Cesnik.
#
# ![](../assets/rotating-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`rotating.ipynb`](@__NBVIEWER_ROOT_URL__/examples/rotating.ipynb).
#-

using GXBeam, LinearAlgebra

sweep = 45 * pi/180
rpm = 0:25:750

## straight section of the beam
L_b1 = 31.5 ## inch
r_b1 = [2.5, 0, 0]
nelem_b1 = 13
lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1)

## swept section of the beam
L_b2 = 6 ## inch
r_b2 = [34, 0, 0]
nelem_b2 = 3
cs, ss = cos(sweep), sin(sweep)
frame_b2 = [cs ss 0; -ss cs 0; 0 0 1]
lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2;
    frame = frame_b2)

## combine elements and points into one array
nelem = nelem_b1 + nelem_b2
points = vcat(xp_b1, xp_b2[2:end])
start = 1:nelem_b1 + nelem_b2
stop = 2:nelem_b1 + nelem_b2 + 1
lengths = vcat(lengths_b1, lengths_b2)
midpoints = vcat(xm_b1, xm_b2)
Cab = vcat(Cab_b1, Cab_b2)

## cross section
w = 1 ## inch
h = 0.063 ## inch

## material properties
E = 1.06e7 ## lb/in^2
ν = 0.325
ρ = 2.51e-4 ## lb sec^2/in^4

## shear and torsion correction factors
ky = 1.2000001839588001
kz = 14.625127919304001
kt = 65.85255016982444

A = h*w
Iyy = w*h^3/12
Izz = w^3*h/12
J = Iyy + Izz

## apply corrections
Ay = A/ky
Az = A/kz
Jx = J/kt

G = E/(2*(1+ν))

compliance = fill(Diagonal([1/(E*A), 1/(G*Ay), 1/(G*Az), 1/(G*Jx), 1/(E*Iyy),
    1/(E*Izz)]), nelem)

mass = fill(Diagonal([ρ*A, ρ*A, ρ*A, ρ*J, ρ*Iyy, ρ*Izz]), nelem)

## create assembly
assembly = Assembly(points, start, stop;
    compliance = compliance,
    mass = mass,
    frames = Cab,
    lengths = lengths,
    midpoints = midpoints)

## create dictionary of prescribed conditions
prescribed_conditions = Dict(
    ## root section is fixed
    1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0)
    )

nonlinear_states = Vector{AssemblyState{Float64}}(undef, length(rpm))
linear_states = Vector{AssemblyState{Float64}}(undef, length(rpm))
for i = 1:length(rpm)
    ## global frame rotation
    w0 = [0, 0, rpm[i]*(2*pi)/60]

    ## perform nonlinear steady state analysis
    system, converged = steady_state_analysis(assembly,
        angular_velocity = w0,
        prescribed_conditions = prescribed_conditions)

    nonlinear_states[i] = AssemblyState(system, assembly;
        prescribed_conditions = prescribed_conditions)

    ## perform linear steady state analysis
    system, converged = steady_state_analysis(assembly,
        angular_velocity = w0,
        prescribed_conditions = prescribed_conditions,
        linear = true)

    linear_states[i] = AssemblyState(system, assembly;
        prescribed_conditions = prescribed_conditions)
end

nothing ##hide

#
# To visualize the solutions we will plot the root moment and tip deflections against the
# angular speed.
#

using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide

#-

#md @suppress_err begin #hide

## root moment
plot(
    xlim = (0, 760),
    xticks = 0:100:750,
    xlabel = "Angular Speed (RPM)",
    yticks = 0.0:2:12,
    ylabel = "\$M_z\$ at the root (lb-in)",
    grid = false,
    overwrite_figure=false
    )
Mz_nl = [-nonlinear_states[i].points[1].M[3] for i = 1:length(rpm)]
Mz_l = [-linear_states[i].points[1].M[3] for i = 1:length(rpm)]
plot!(rpm, Mz_nl, label="Nonlinear")
plot!(rpm, Mz_l, label="Linear")
#!nb plot!(show=true)
#md savefig("../assets/rotating-Mz.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/rotating-Mz.svg)

#-

#md @suppress_err begin #hide

## x tip deflection
plot(
    xlim = (0, 760),
    xticks = 0:100:750,
    xlabel = "Angular Speed (RPM)",
    ylim = (-0.002, 0.074),
    yticks = 0.0:0.01:0.07,
    ylabel = "\$u_x\$ at the tip (in)",
    grid = false,
    overwrite_figure=false
    )
ux_nl = [nonlinear_states[i].points[end].u[1] for i = 1:length(rpm)]
ux_l = [linear_states[i].points[end].u[1] for i = 1:length(rpm)]
plot!(rpm, ux_nl, label="Nonlinear")
plot!(rpm, ux_l, label="Linear")

#!nb plot!(show=true)
#md savefig("../assets/rotating-ux.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/rotating-ux.svg)

#-

#md @suppress_err begin #hide

## y tip deflection
plot(
    xlim = (0, 760),
    xticks = 0:100:750,
    xlabel = "Angular Speed (RPM)",
    ylim = (-0.01, 0.27),
    yticks = 0.0:0.05:0.25,
    ylabel = "\$u_y\$ at the tip (in)",
    grid = false,
    overwrite_figure=false
    )
uy_nl = [nonlinear_states[i].points[end].u[2] for i = 1:length(rpm)]
uy_l = [linear_states[i].points[end].u[2] for i = 1:length(rpm)]
plot!(rpm, uy_nl, label="Nonlinear")
plot!(rpm, uy_l, label="Linear")

#!nb plot!(show=true)
#md savefig("../assets/rotating-uy.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/rotating-uy.svg)

#-

#md @suppress_err begin #hide

## rotation of the tip
plot(
    xlim = (0, 760),
    xticks = 0:100:750,
    xlabel = "Angular Speed (RPM)",
    ylabel = "\$θ_z\$ at the tip",
    grid = false,
    overwrite_figure=false
    )
theta_z_nl = [4*atan(nonlinear_states[i].points[end].theta[3]/4)
    for i = 1:length(rpm)]
theta_z_l = [4*atan(linear_states[i].points[end].theta[3]/4)
    for i = 1:length(rpm)]

plot!(rpm, theta_z_nl, label="Nonlinear")
plot!(rpm, theta_z_l, label="Linear")

#!nb plot!(show=true)
#md savefig("../assets/rotating-theta_z.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/rotating-theta_z.svg)

#-

# We will now compute the eigenvalues of this system for a range of sweep angles and and
# angular speeds.

sweep = (0:2.5:45) * pi/180
rpm = [0, 500, 750]
nev = 30

λ = Matrix{Vector{ComplexF64}}(undef, length(sweep), length(rpm))
U = Matrix{Matrix{ComplexF64}}(undef, length(sweep), length(rpm))
MV = Matrix{Matrix{ComplexF64}}(undef, length(sweep), length(rpm))
state = Matrix{AssemblyState{Float64}}(undef, length(sweep), length(rpm))
eigenstates = Matrix{Vector{AssemblyState{ComplexF64}}}(undef,
    length(sweep), length(rpm))
for i = 1:length(sweep)

    local L_b1, r_b1, nelem_b1, lengths_b1
    local xp_b1, xm_b1, Cab_b1
    local cs, ss
    local L_b2, r_b2, nelem_b2, frame_b2, lengths_b2
    local xp_b2, xm_b2, Cab_b2
    local nelem, points, start, stop
    local lengths, midpoints, Cab, compliance, mass, assembly

    ## straight section of the beam
    L_b1 = 31.5 # inch
    r_b1 = [2.5, 0, 0]
    nelem_b1 = 20
    lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1)

    ## swept section of the beam
    L_b2 = 6 # inch
    r_b2 = [34, 0, 0]
    nelem_b2 = 20
    cs, ss = cos(sweep[i]), sin(sweep[i])
    frame_b2 = [cs ss 0; -ss cs 0; 0 0 1]
    lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2;
        frame = frame_b2)

    ## combine elements and points into one array
    nelem = nelem_b1 + nelem_b2
    points = vcat(xp_b1, xp_b2[2:end])
    start = 1:nelem_b1 + nelem_b2
    stop = 2:nelem_b1 + nelem_b2 + 1
    lengths = vcat(lengths_b1, lengths_b2)
    midpoints = vcat(xm_b1, xm_b2)
    Cab = vcat(Cab_b1, Cab_b2)

    compliance = fill(Diagonal([1/(E*A), 1/(G*Ay), 1/(G*Az), 1/(G*Jx),
        1/(E*Iyy), 1/(E*Izz)]), nelem)

    mass = fill(Diagonal([ρ*A, ρ*A, ρ*A, ρ*J, ρ*Iyy, ρ*Izz]), nelem)

    ## create assembly
    assembly = Assembly(points, start, stop;
        compliance = compliance,
        mass = mass,
        frames = Cab,
        lengths = lengths,
        midpoints = midpoints)

    ## create system
    system = System(assembly, false)

    for j = 1:length(rpm)
        ## global frame rotation
        w0 = [0, 0, rpm[j]*(2*pi)/60]

        ## eigenvalues and (right) eigenvectors
        system, λ[i,j], V, converged = eigenvalue_analysis!(system, assembly;
            angular_velocity = w0,
            prescribed_conditions = prescribed_conditions,
            nev = nev)

        ## corresponding left eigenvectors
        U[i,j] = left_eigenvectors(system, λ[i,j], V)

        ## post-multiply mass matrix with right eigenvector matrix
        ## (we use this later for correlating eigenvalues)
        MV[i,j] = system.M * V

        ## process state and eigenstates
        state[i,j] = AssemblyState(system, assembly;
            prescribed_conditions = prescribed_conditions)
        eigenstates[i,j] = [AssemblyState(system, assembly, V[:,k];
            prescribed_conditions = prescribed_conditions) for k = 1:nev]
    end
end

#!jl nothing #hide

#
# We can correlate each eigenmode by taking advantage of the fact that left and right
# eigenvectors satisfy the following relationships:
#
# ```math
# \begin{aligned}
# u^H M v &= 1 &\text{if \(u\) and \(v\) correspond to the same eigenmode} \\
# u^H M v &= 0 &\text{if \(u\) and \(v\) correspond to different eigenmodes}
# \end{aligned}
# ```

## set previous left eigenvector matrix
U_p = copy(U[1,1])

for j = 1:length(rpm)
    for i = 1:length(sweep)
        ## construct correlation matrix
        C = U_p*MV[i,j]

        ## correlate eigenmodes
        perm, corruption = correlate_eigenmodes(C)

        ## re-arrange eigenvalues and eigenvectors
        λ[i,j] = λ[i,j][perm]
        U[i,j] = U[i,j][perm,:]
        MV[i,j] = MV[i,j][:,perm]
        eigenstates[i,j] = eigenstates[i,j][perm]

        ## update previous eigenvector matrix
        U_p .= U[i,j]
    end
    ## update previous eigenvector matrix
    U_p .= U[1,j]
end

frequency = [[imag(λ[i,j][k])/(2*pi) for i = 1:length(sweep), j=1:length(rpm)]
    for k = 1:2:nev]

#!jl nothing #hide

# In this case these eigenmode correlations work, but remember that large changes in the
# underlying parameters (or just drastic changes in the eigenvectors themselves due to a
# small perturbation) can cause these automatic eigenmode correlations to fail.

# We'll now plot the frequency of the different eigenmodes against those found by Epps and
# Chandra in "The Natural Frequencies of Rotating Composite Beams With Tip Sweep".

names = ["First Bending Mode", "Second Bending Mode", "Third Bending Mode"]
indices = [1, 2, 4]

experiment_rpm = [0, 500, 750]
experiment_sweep = [0, 15, 30, 45]
experiment_frequencies = [
    [1.4 1.8 1.7 1.6;
     10.2 10.1 10.2 10.2;
     14.8 14.4 14.9 14.7],
    [10.3 10.2 10.4 10.4;
     25.2 25.2 23.7 21.6;
     36.1 34.8 30.7 26.1],
    [27.7 27.2 26.6 24.8;
     47.0 44.4 39.3 35.1;
     62.9 55.9 48.6 44.8]
]
#!jl nothing #hide

#-

#md @suppress_err begin #hide
#nb ph = Vector{Any}(undef, 3)

for k = 1:3
#nb     ph[k] = plot(
#!nb     plot(
        title = names[k],
        xticks = 0:15:45,
        xlabel = "Sweep Angle (degrees)",
        ylim = (0, Inf),
        ylabel = "Frequency (Hz)",
        grid = false,
        overwrite_figure=false
        )

    plot!([], [], color=:black, label="GXBeam")
    scatter!([], [], color=:black, label = "Experiment (Epps and Chandra)")

    for j = length(rpm):-1:1
        plot!(sweep*180/pi, frequency[indices[k]][:,j],
            label="$(rpm[j]) RPM", color=j)
        scatter!(experiment_sweep, experiment_frequencies[k][j,:],
            label="", color=j)
    end

#!nb     plot!(show=true)
#md savefig("../assets/rotating-frequencies-$(k).svg") #hide
#md closeall() #hide
end
#md end #hide
#md nothing #hide

#nb ph[1]
#nb #-
#nb ph[2]
#nb #-
#nb ph[3]
#nb #-

#md # ![](../assets/rotating-frequencies-1.svg)
#md # ![](../assets/rotating-frequencies-2.svg)
#md # ![](../assets/rotating-frequencies-3.svg)

#-

#md @suppress_err begin #hide

names = ["1T/5B", "5B/1T", "4B/1T"]
indices = [5, 7, 6]

experiment_frequencies = [
    95.4 87.5 83.7 78.8;
    106.6 120.1 122.6 117.7;
    132.7 147.3 166.2 162.0
]

plot(
    title = "Coupled Torsion-Bending Modes at 750 RPM",
    xticks = 0:15:45,
    xlabel = "Sweep Angle (degrees)",
    ylim = (0, Inf),
    ylabel = "Frequency (Hz)",
    legend = :bottomleft,
    grid = false,
    overwrite_figure=false
    )

plot!([], [], color=:black, label="GXBeam")
scatter!([], [], color=:black,
    label="Experiment (Epps and Chandra)")

for k = 1:length(indices)
    plot!(sweep*180/pi, frequency[indices[k]][:,end], label=names[k], color=k)
    scatter!(experiment_sweep, experiment_frequencies[k,:], label="", color=k)
end
#!nb plot!(show=true)
#md savefig("../assets/rotating-frequencies-4.svg");
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/rotating-frequencies-4.svg)

#-

# As you can see, the frequency results from the eigenmode analysis in this package
# compare well with experimental results.
#
# We can also visualize eigenmodes using ParaView.  Here we will visualize the first
# bending mode for the 45 degree swept tip at a rotational speed of 750 RPM.  This can be
# helpful for identifying different eigenmodes.

## write the response to vtk files for visualization using ParaView
write_vtk("rotating-eigenmode", assembly, state[end,end],
    λ[end,end][1], eigenstates[end,end][1]; mode_scaling = 100.0)

# ![](../assets/rotating-eigenmode.gif)

ext: jl


text:
"""
    basicblock(inplanes, outplanes, downsample = false)

Create a basic residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
basicblock(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = 2, pad = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2], identity; stride = 1, pad = 1, bias = false)...) :
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = 1, pad = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2], identity; stride = 1, pad = 1, bias = false)...)

"""
    bottleneck(inplanes, outplanes, downsample = false)

Create a bottleneck residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
bottleneck(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = 2, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, bias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = 1, bias = false)...) :
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = 1, bias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, bias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = 1, bias = false)...)

"""
    skip_projection(inplanes, outplanes, downsample = false)

Create a skip projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: the number of output feature maps
- `downsample`: set to `true` to downsample the input
"""
skip_projection(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 2, bias = false)...) :
  Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 1, bias = false)...)

# array -> PaddedView(0, array, outplanes) for zero padding arrays
"""
    skip_identity(inplanes, outplanes[, downsample])

Create a identity projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: the number of output feature maps
- `downsample`: this argument is ignored but it is needed for compatibility with [`resnet`](#).
"""
function skip_identity(inplanes, outplanes)
  if outplanes > inplanes
    return Chain(MaxPool((1, 1), stride = 2),
                 y -> cat(y, zeros(eltype(y),
                                   size(y, 1),
                                   size(y, 2),
                                   outplanes - inplanes, size(y, 4)); dims = 3))
  else
    return identity
  end
end
skip_identity(inplanes, outplanes, downsample) = skip_identity(inplanes, outplanes)

"""
    resnet(block, residuals::NTuple{2, Any}, connection = (x, y) -> @. relu(x) + relu(y);
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `residuals`: a 2-tuple of functions with input `(inplanes, outplanes, downsample=false)`,
               each of which will return a function that will be used as a new "skip" path to match a residual block.
              [`Metalhead.skip_identity`](#) and [`Metalhead.skip_projection`](#) can be used here.
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(block, residuals::NTuple{2, Any}, connection = (x, y) -> @. relu(x) + relu(y);
                channel_config, block_config, nclasses = 1000)
  inplanes = 64
  baseplanes = 64
  layers = []
  append!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = (3, 3)))
  push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
  for (i, nrepeats) in enumerate(block_config)
    # output planes within a block
    outplanes = baseplanes .* channel_config
    # push first skip connection on using first residual
    # downsample the residual path if this is the first repetition of a block
    push!(layers, Parallel(connection, block(inplanes, outplanes, i != 1),
                                       residuals[1](inplanes, outplanes[end], i != 1)))
    # push remaining skip connections on using second residual
    inplanes = outplanes[end]
    for _ in 2:nrepeats
      push!(layers, Parallel(connection, block(inplanes, outplanes, false),
                                         residuals[1](inplanes, outplanes[end], false)))
      inplanes = outplanes[end]
    end
    # next set of output plane base is doubled
    baseplanes *= 2
  end

  return Chain(Chain(layers...),
               Chain(AdaptiveMeanPool((1, 1)), flatten, Dense(inplanes, nclasses)))
end

"""
    resnet(block, shortcut_config::Symbol, connection = (x, y) -> @. relu(x) + relu(y);
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
    - `:A`: uses a [`Metalhead.skip_identity`](#) for all residual blocks
    - `:B`: uses a [`Metalhead.skip_projection`](#) for the first residual block
            and [`Metalhead.skip_identity`](@) for the remaining residual blocks
    - `:C`: uses a [`Metalhead.skip_projection`](#) for all residual blocks
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
resnet(block, shortcut_config::Symbol, args...; kwargs...) =
  (shortcut_config == :A) ? resnet(block, (skip_identity, skip_identity), args...; kwargs...) :
  (shortcut_config == :B) ? resnet(block, (skip_projection, skip_identity), args...; kwargs...) :
  (shortcut_config == :C) ? resnet(block, (skip_projection, skip_projection), args...; kwargs...) :
  error("Unrecognized shortcut config == $shortcut_config passed to resnet (use :A, :B, or :C).")

const resnet_config =
  Dict(:resnet18 => ([1, 1], [2, 2, 2, 2], :A),
       :resnet34 => ([1, 1], [3, 4, 6, 3], :A),
       :resnet50 => ([1, 1, 4], [3, 4, 6, 3], :B),
       :resnet101 => ([1, 1, 4], [3, 4, 23, 3], :B),
       :resnet152 => ([1, 1, 4], [3, 8, 36, 3], :B))

"""
    ResNet(channel_config, block_config, shortcut_config; block, nclasses = 1000)

Create a `ResNet` model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`resnet`](#).

# Arguments
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `nclasses`: the number of output classes
"""
struct ResNet
  layers
end

function ResNet(channel_config, block_config, shortcut_config; block, nclasses = 1000)
  layers = resnet(block,
                  shortcut_config;
                  channel_config = channel_config,
                  block_config = block_config,
                  nclasses = nclasses)

  ResNet(layers)
end

@functor ResNet

(m::ResNet)(x) = m.layers(x)

backbone(m::ResNet) = m.layers[1]
classifier(m::ResNet) = m.layers[2]

"""
    ResNet18(; pretrain = false, nclasses = 1000)

Create a ResNet-18 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `nclasses`: the number of output classes

!!! warning
    `ResNet18` does not currently support pretrained weights.
"""
function ResNet18(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet18]...; block = basicblock, nclasses = nclasses)

  pretrain && loadpretrain!(model, "ResNet18")
  return model
end

"""
    ResNet34(; pretrain = false, nclasses = 1000)

Create a ResNet-34 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet34` does not currently support pretrained weights.
"""
function ResNet34(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet34]...; block = basicblock, nclasses = nclasses)

  pretrain && loadpretrain!(model, "ResNet34")
  return model
end

"""
    ResNet50(; pretrain = false, nclasses = 1000)

Create a ResNet-50 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet50` does not currently support pretrained weights.
"""
function ResNet50(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet50]...; block = bottleneck, nclasses = nclasses)

  pretrain && loadpretrain!(model, "ResNet50")
  return model
end

"""
    ResNet101(; pretrain = false, nclasses = 1000)

Create a ResNet-101 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet101` does not currently support pretrained weights.
"""
function ResNet101(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet101]...; block = bottleneck, nclasses = nclasses)

  pretrain && loadpretrain!(model, "ResNet101")
  return model
end

"""
    ResNet152(; pretrain = false, nclasses = 1000)

Create a ResNet-152 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet152` does not currently support pretrained weights.
"""
function ResNet152(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet152]...; block = bottleneck, nclasses = nclasses)

  pretrain && loadpretrain!(model, "ResNet152")
  return model
end


text:
@enum FileFormat AUTO ASCII ASCII_MATLAB BINARY BINARY_FLAT BMP JPG JPEG MATRIX_MARKET PNG PPM XBM XPM FileFormat_MAX

for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))
    for mattype in ("", "Dist")
        mat = Symbol(mattype, "Matrix")
        @eval begin
            function print(A::$mat{$elty}, title::String)
                ElError(ccall(($(string("ElPrint", mattype, "_", ext)), libEl), Cuint,
                    (Ptr{Cvoid}, Ptr{UInt8}),
                    A.obj, title))
            end

            function write(A::$mat{$elty}, basename::String, format::FileFormat, title::String)
                ElError(ccall(($(string("ElWrite", mattype, "_", ext)), libEl), Cuint,
                    (Ptr{Cvoid}, Ptr{UInt8}, FileFormat, Ptr{UInt8}),
                    A.obj, basename, format, title))
            end
        end
    end
end


text:
using Pkg
Pkg.instantiate()

using PackageCompiler

create_sysimage(:MusicAlbums;
    sysimage_path = "MusicAlbums.so",
    precompile_execution_file = "deploy/pkgprecompile.jl")


text:
@userplot PlotMode

function get_mode_title(pseudo, approx, z, s)
    prefix = pseudo ? "Pseudomode: " : "Eigenmode: "
    rel = approx ? "\\approx" : "="
    formatted = @sprintf("%15.2e", s)
    infix = pseudo ? ("\$\\epsilon = " * formatted * "\$") : ("\$\\kappa(\\lambda) " * rel * formatted * "\$")
    λ_str = @sprintf("%12g%+12gi", real(z), imag(z))
    return prefix * infix * "\n  \$\\lambda = " * λ_str * "\$"
end

"""
    plotmode(A, z, U=I; ...)

Plot the (pseudo)eigenmode of `A` associated with (pseudo)eigenvalue `z`.

# Arguments
* `A` an arbitrary square matrix
* `z` the real or complex (pseudo)eigenvalue whose eigenmode is to be plotted.
* `U` optional Schur factor. The eigenmode is premultiplied by `U` before plotting.
# Keyword arguments
* `pseudo::Bool = false` Whether to plot the eigenmode or pseudomode.
* `approx::Bool = false`
* `verbosity::Int = 0` Verbosity parameter passed when computing the eigenmode.

# Examples
```jldoctest
julia> A = Pseudospectra.grcar(40)
julia> plotmode(A, 0.5+2.0im)
```
"""
function plotmode end

RecipesBase.@recipe function f(p::PlotMode, U = I; pseudo=false, approx=false, verbosity=0)
    if length(p.args) != 2
        error("plotmode must be given 2 arguments")
    end
    A = p.args[1]
    z = p.args[2]

    SS, q = get_mode(A, z, pseudo, verbosity)

    q = U * q

    #showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()

    # NB. The title is broken
    #title_str = get_mode_title(pseudo, approx, z, SS[end])
    #title := title_str

    @series begin
        seriestype := :path
        linecolor := (pseudo ? :magenta : :cyan)
        label := "Real part"
        real.(q)
    end

    seriestype := :path
    linecolor := :black
    linestyle := :dash

    @series begin
        label := "Modulus"
        abs.(q)
    end

    @series begin
        label := ""
        -abs.(q)
    end
end


text:
"""
# module Directives

# Examples

```jldoctest
julia>
```
"""
module Directives

using Parameters

export Directive

@with_kw struct Directive{A <: AbstractString,B <: Union{Nothing,AbstractChar},T}
    long::A
    short::B = nothing
    value::T
end

function Base.string(d::Directive; long::Bool = true, with_value::Bool = true)
    key = if long
        "--$(d.long)"
    else
        isnothing(d.short) ? throw(ArgumentError("Directive `$(d.long)` does not have a short form!")) : "-$(d.short)"
    end
    value = if with_value
        long ? "==$(d.value)" : " $(d.value)"
    else
        ""
    end
    return key * value
end

end


text:
# Vertex
vertex(n::AbstractNode) = n.vertex[]
set_vertex!(n::AbstractNode, v::Int) = n.vertex[] = v

# Ports
freeports(n::AbstractNode) = n.freeports
freeports(::Junction) = [true]
numports(n::AbstractNode) = length(n.freeports)
numports(::Junction) = Inf
updateport!(n::AbstractNode, idx::Int) = freeports(n)[idx] = !freeports(n)[idx]
nextfreeport(n::AbstractNode) = findfirst(freeports(n))

# Nodes in Bonds
srcnode(b::Bond) = b.src.node
dstnode(b::Bond) = b.dst.node
in(n::AbstractNode, b::Bond) = n === srcnode(b) || n === dstnode(b)

# Searching
getnodes(bg::BondGraph, t::Symbol) = filter(x -> x.type == t, bg.nodes)
getnodes(bg::BondGraph, n::AbstractString) = filter(x -> x.name == n, bg.nodes)

getbonds(bg::BondGraph, t::Tuple) = getbonds(bg, t[1], t[2])
getbonds(bg::BondGraph, node1::AbstractNode, node2::AbstractNode) = filter(b -> node1 in b && node2 in b, bg.bonds)

# I/O
show(io::IO, node::AbstractNode) = print(io, "$(node.type):$(node.name)")
show(io::IO, node::Junction) = print(io, "$(node.type)")
show(io::IO, port::Port) = print(io, "Port $(port.node) ($(port.index))")
show(io::IO, b::Bond) = print(io, "Bond $(srcnode(b)) ⇀ $(dstnode(b))")
show(io::IO, bg::BondGraph) = print(io, "BondGraph $(bg.name) ($(lg.nv(bg)) Nodes, $(lg.ne(bg)) Bonds)")

# Comparisons
# These definitions will need to expand when equations etc. are added
==(n1::AbstractNode, n2::AbstractNode) = n1.type == n2.type && n1.name == n2.name
# bondgraph comparison

# Mapping
# TODO - function which creates map of components and bonds to indices

# Easier referencing systems using a.b notation
# TODO create tests
function getproperty(bg::BondGraph, sym::Symbol)
    # Calling getfield explicitly avoids using "a.b"
    # and causing a StackOverflowError
    allnodes = getfield(bg, :nodes)
    names = [getfield(n, :name) for n in allnodes]
    symnodes = allnodes[names .== string(sym)]
    if isempty(symnodes)
        return getfield(bg, sym)
    elseif length(symnodes) == 1
        return symnodes[1]
    else
        return symnodes
    end
end

function getproperty(bgn::BondGraphNode, sym::Symbol)
    bg = getfield(bgn, :bondgraph)
    try
        return getproperty(bg, sym)
    catch
        return getfield(bgn, sym)
    end
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
function add_unique_store!(ls::LoopSet, op::Operation)
  add_vptr!(ls, op)
  pushop!(ls, op, name(op.ref))
end
function cse_store!(ls::LoopSet, op::Operation)
  id = identifier(op)
  ls.operations[id] = op
  ls.opdict[op.variable] = op
end
function add_store!(ls::LoopSet, op::Operation, add_pvar::Bool = !any(r -> r == op.ref, ls.refs_aliasing_syms))
  @assert isstore(op)
  if add_pvar
    push!(ls.syms_aliasing_refs, name(first(parents(op))))
    push!(ls.refs_aliasing_syms, op.ref)
  end
  id = op.identifier
  id == length(operations(ls)) ? add_unique_store!(ls, op) : cse_store!(ls, op)
end
function add_copystore!(
    ls::LoopSet, parent::Operation, mpref::ArrayReferenceMetaPosition, elementbytes::Int
)
    op = add_compute!(ls, gensym!(ls, "identity"), :identity, [parent], elementbytes)
    # pushfirst!(mpref.parents, parent)
    add_store!(ls, mpref, elementbytes, op)
end

function add_store!(
    ls::LoopSet, mpref::ArrayReferenceMetaPosition, elementbytes::Int, parent = getop(ls, varname(mpref), mpref.loopdependencies, elementbytes)
)
    isload(parent) && return add_copystore!(ls, parent, mpref, elementbytes)
    vparents = mpref.parents
    ldref = mpref.loopdependencies
    reduceddeps = mpref.reduceddeps
    pvar = name(parent)
    id = length(ls.operations)
    # try to cse store, by replacing the previous one
    mref = mpref.mref
    add_pvar = true
    for opp ∈ operations(ls)
        if mref == opp.ref
            isstore(opp) && (id = opp.identifier)
            add_pvar = false
            break
        end
        # add_pvar &= (name(first(parents(opp))) != pvar)
    end
    pushfirst!(vparents, parent)
    update_deps!(ldref, reduceddeps, parent)
    op = Operation( id, name(mpref), elementbytes, :setindex!, memstore, mpref )
    add_store!(ls, op, add_pvar)
end

function add_store!(
    ls::LoopSet, var::Symbol, array::Symbol, rawindices, elementbytes::Int
)
    mpref = array_reference_meta!(ls, array, rawindices, elementbytes, var)
    add_store!(ls, mpref, elementbytes)
end
function add_simple_store!(ls::LoopSet, parent::Operation, mref::ArrayReferenceMeta, elementbytes::Int)
    op = Operation( ls, name(mref), elementbytes, :setindex!, memstore, getindices(mref.ref), NODEPENDENCY, [parent], mref )
    add_unique_store!(ls, op)
end
function add_simple_store!(ls::LoopSet, var::Union{Symbol,Operation}, ref::Union{ArrayReference,ArrayReferenceMeta}, elementbytes::Int)
    parent = isa(var, Symbol) ? getop(ls, var, elementbytes) : var
    mref = isa(ref, ArrayReference) ? ArrayReferenceMeta(ref, fill(true, length(getindices(ref)))) : ref
    add_simple_store!(ls, parent, mref, elementbytes)
end
function add_store_ref!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int)
    array, raw_indices = ref_from_ref!(ls, ex)
    add_store!(ls, var, array, raw_indices, elementbytes)
end
function add_store_ref!(ls::LoopSet, var, ex::Expr, elementbytes::Int)
    array, raw_indices = ref_from_ref!(ls, ex)
    mpref = array_reference_meta!(ls, array, raw_indices, elementbytes)
    c = add_constant!(ls, var, loopdependencies(mpref), gensym(:storeconst), elementbytes)
    add_store!(ls, mpref, elementbytes, c)
end

# For now, it is illegal to load from a conditional store.
# if you want that sort of behavior, do a conditional reassignment, and store that result unconditionally.
function add_conditional_store!(ls::LoopSet, LHS, condop::Operation, storeop::Operation, elementbytes::Int)
    array, rawindices = ref_from_ref!(ls, LHS)
    mpref = array_reference_meta!(ls, array, rawindices, elementbytes)
    mref = mpref.mref
    ldref = mpref.loopdependencies

    pvar = name(storeop)
    id = length(ls.operations)
    @assert pvar ∉ ls.syms_aliasing_refs
    # if pvar ∉ ls.syms_aliasing_refs
  # FIXME properly handle CSE of conditional stores.
    push!(ls.syms_aliasing_refs, pvar)
    push!(ls.refs_aliasing_syms, mref)
    storeparents = [storeop, condop]
    # else
    #     # for now, we don't try to cse the store
    #     # later, as an optimization, we could:
    #     # 1. cse the store
    #     # 2. use the mask to combine the vector we're trying to store here with the vector that would have been stored in the now cse-ed 1.
    #     # 3. use a regular (non-masked) store on that vector.
    #     ref = mpref.mref.ref
    #     for opp ∈ operations(ls)
    #         isstore(opp) || continue
    #         if ref == opp.ref.ref# && return cse_store!(ls, identifier(opp), mref, parents, ldref, reduceddeps, elementbytes)
    #             id = opp.identifier
    #             break
    #         end
    #     end
    #     if id != length(ls.operations) # then there was a previous store
    #         prevstore = getop(ls, id + 1)
    #         # @show prevstore prevstore.node_type, loopdependencies(prevstore)
    #         # @show operations(ls)
    #         storeop = add_compute!(ls, gensym(:combinedstoreop), Instruction(:ifelse), [condop, storeop, first(parents(prevstore))], elementbytes)
    #         storeparents = [storeop]
    #         storeinstr = if prevstore.instruction.instr === :conditionalstore!
    #             push!(storeparents, add_compute!(ls, gensym(:combinedmask), Instruction(:|), [condop, last(parents(prevstore))], elementbytes))
    #             :conditionalstore!
    #         else
    #             :setindex!
    #         end
    #         op = Operation( id, name(mref), elementbytes, storeinstr, memstore, ldref, NODEPENDENCY, storeparents, mref )
    #         return cse_store!(ls, op)
    #     end
    # end
    op = Operation( id, name(mref), elementbytes, :conditionalstore!, memstore, ldref, reduceddependencies(storeop), storeparents, mref )
    add_unique_store!(ls, op)
end


text:
@doc """
    MPSState <: QState
# Structure
- `sites::Vector{ITensor}`: Array of ITensors storing MPS sites.
- `map::QubitSiteMap`: A binary map between qubit and site(position).
- `llim::Int`: Sites to the left of llim is guaranteed to be orthogonal.
- `rlim::Int`: Sites to the right of rlim is guaranteed to be orthogonal.
""" ->
MPSState

@doc """
	MPSState(N::Int, init::Vector)

Construct a MPSState with `N` qubits all initialized to given `init` vector.

# Arguments
- `N::Int`: the number of qubits
- `init::Vector{T} where T`: the unit vector representing the initial state of each qubit

""" ->
MPSState(N::Int, init::Vector)

@doc """
	MPSState(N::Int)

Construct |0> MPSState with `N` qubits.

Same as `MPSState(N, [1,0])`
""" ->
MPSState(N::Int)

@doc """
	getindex(state::MPSState, n::Int)

Return the ITensor at `n`th site in the MPSState.

Note: `n`th site is not the same as `n`th qubit.
""" ->
getindex(state::MPSState, n::Int)

@doc """
	getinds(state::MPSState, inds::Vector{Int})

Return ITensors at sites with indices provided by `inds` vector.

""" ->
getinds(state::MPSState, inds::Vector{Int})

@doc """
	setindex!(state::MPSState, T::ITensor, n::Integer)

Set the `n`th site of `state` to be ITensor `T`.

#Note:
No error checking. Should be careful when calling this function.
""" ->
setindex!(state::MPSState, T::ITensor, n::Integer)

@doc """
	length(state::MPSState)

Return number of qubits.

""" ->
length(state::MPSState)

@doc """
	size(state::MPSState)
Return number of qubits.
""" ->
size(state::MPSState)

@doc """
	iterate(state::MPSState, itstate::Int=1)

Iterate through sites.
""" ->
iterate(state::MPSState, itstate::Int=1)

@doc """
	copy(state::MPSState)
Make a copy of `state`.
""" ->
copy(state::MPSState)

@doc """
	isapprox(state1::MPSState, state2::MPSState)
Inexact comparison between two MPSStates `state1` and `state2`.
""" ->
isapprox(state1::MPSState, state2::MPSState)

@doc """
	siteForQubit(state::MPSState, i::Int)
Return site index for qubit `i`.
""" ->
siteForQubit(state::MPSState, i::Int)

@doc """
	sitesForQubits(state::MPSState, inds::Vector{Int})

Return an array of site indices for given qubits.
""" ->
sitesForQubits(state::MPSState, inds::Vector{Int})

@doc """
	qubitAtSite(state::MPSState, i::Int)
Get which qubit is currently stored at site `n` in `state`.
""" ->
qubitAtSite(state::MPSState, i::Int)

@doc """
	qubitsAtSites(state::MPSState, inds::Vector{Int})
Get which qubits are currently stored at sites `inds` in `state`.
""" ->
qubitsAtSites(state::MPSState, inds::Vector{Int})

@doc """
	updateMap!(state::MPSState; kwargs...)
Update state.map according to qubit-site pair proveided in keyword arguments.

# Arguments
- `state`: MPSState
- `q`: [_keyword argument_] Qubit index. Default value: -1.
- `s`: [_keyword argument_] Site index. Default value: -1.

# Note:
This function is intended for internal use.
""" ->
updateMap!(state::MPSState; kwargs...)

@doc """
	function updateLims!(state::MPSState, leftEnd::Int, rightEnd::Int)
Update `state.llim` and `state.rlim`.

# Arguments
- `state` : MPSState
- `leftEnd`: _Site_ index of the most left modified site.
- `rightEnd`: _Site_ index of the most right modified site.

# Note:
This function is intended for internal use.
""" ->
updateLims!(state::MPSState, leftEnd::Int, rightEnd::Int)

@doc """
	updateLims!(state::MPSState, touchedSite::Int)
Update `state.llim` and `state.rlim` according to modified site.

# Note:
This function is intended for internal use.
""" ->
updateLims!(state::MPSState, touchedSite::Int)

@doc """
	getQubit(state::MPSState, i::Int)
Return ITensor for qubit with index `i`.
""" ->
getQubit(state::MPSState, i::Int)

@doc """
	getQubits(state::MPSState, inds::Vector{Int})
Return ITensors for qubits specified in `inds`.
""" ->
getQubits(state::MPSState, inds::Vector{Int})

@doc """
	siteIndex(state::MPSState, i::Int)
Return site index in the ITensor at site `i`.
# Arguments:
- `state`: MPSState
- `i`: Site index

# Return Value:
Index object.

# Note:
This function is intended for internal use.
""" ->
siteIndex(state::MPSState, i::Int)

@doc """
	siteInds(state::MPSState, inds::Vector{Int})
Return site indices in the ITensor at sites `inds`.

# Arguments:
- `state`: MPSState
- `inds`: Site indices

# Return Value:
IndexSet object.

# Note:
This function is intended for internal use.
""" ->
siteInds(state::MPSState, inds::Vector{Int})

@doc """
	toExactState(state::MPSState)
Turn a MPSState to an ExactState.

# Note:
This is a costly operation and is intended for MPSState testing.
""" ->
toExactState(state::MPSState)

@doc """
	normalize!(state::MPSState; kwargs...)
Normalize the MPSState by getting it into canonical form and normalize at the center.

# Arguments
- `state`: MPSState
- `center`: [_keyword argument_] Site index for gauge center. Default value: state.llim+1.
""" ->
normalize!(state::MPSState; kwargs...)

@doc """
	measure!(state::MPSState, qubits::Vector{Int}, shots::Int; kwargs...)
Perform measurement to specified `qubits` in computatioal basis for `shots` times.

# Arguments
- `state`: MPSState
- `qubits`: vector of qubits to be measured
- `shots`: number of measurements to be performed
- `binary`: [_keyword argument_] If `binary == true`, the function returns results in binary form in a 2D array. Otherwise, the function returns measurement results in decimal form. Default value: `false`.
- `probability`: [_keyword argument_] whether or not to keep probability as a by-product.
- `logprob`: [_keyword argument_] switch for getting log(prob) or direct prob (handy when probability is smal
l for all configurations). Default value: `false`.

# Example
``` julia
julia> state = MPSState(3) # generate 3-qubit |0> state
3-qubit MPSState

julia> apply!(state, H, 1) # apply Hadamard gate at qubit 1
3-qubit MPSState

julia> result = measure!(state, [1:3;], 5)
5×3 Array{Int64,2}:
 0  0  0
 1  0  0
 0  0  0
 1  0  0
 0  0  0

julia> result[4,1] # outcome of 1st qubit at 4th measurement
1

julia> measure!(state, [1:3;], 5; binary = false)
5-element Array{Int64,1}:
 0
 0
 1
 1
 1

julia> measure!(state, [1:3;], 3; binary = true, probability = true)
([0 0 0; 0 0 0; 1 0 0], [0.5, 0.5, 0.5])
```
 # Note
Calling this function will not physically change state, but the internal data structure (e.g. the position of gauge center) of the MPS may be changed.
See also: **collapse!**
""" ->
measure!(state::MPSState, qubits::Vector{Int}, shots::Int; kwargs...)

@doc """
	measure!(state::MPSState, qubit::Int, shots::Int; kwargs...)
Perform measuremt to specified qubit in computational basis.

# Arguments
- `state`: MPSState
- `qubit`: qubit to be measured
- `shots`: number of measurements to be performed
- `binary`: [_keyword argument_] If `binary == true`, the function returns results in binary form in a 2D array. Otherwise, the function returns measurement results in decimal form. Default value: `false`.
""" ->
measure!(state::MPSState, qubit::Int, shots::Int; kwargs...)

@doc """
	orderQubits!(state::MPSState; kwargs...)
Reorder internal storage of MPSState; make site index match qubit index.
""" ->
orderQubits!(state::MPSState; kwargs...)

@doc """
	moveQubit!(state::MPSState, q::Int, s::Int; kwargs...)
Move qubit `q` to site `s`. `kwargs` are used for svd truncation (See **svd**.)
""" ->
moveQubit!(state::MPSState, q::Int, s::Int; kwargs...)

@doc """
	localizeQubitsInOrder!(state::MPSState, qubits::Vector{Int}; kwargs...)
Localize `qubits` respecting the given ordering.

The specified qubits are moved to the front of MPS by calling **moveQubit!**.
""" ->
localizeQubitsInOrder!(state::MPSState, qubits::Vector{Int}; kwargs...)

@doc """
	centerAtSite!(state::MPSState, s::Int)
Canonicalize MPS and put the canonical center at site `s`.

# Note:
MPS will be auto-normalized after canonicalization.
""" ->
centerAtSite!(state::MPSState, s::Int)

@doc """
	centerAtQubit!(state::MPSState, q::Int)
Canonicalize MPS and put the canonical center at qubit `q`.
""" ->
centerAtQubit!(state::MPSState, q::Int)

@doc """
	swapSites!(state::MPSState, s1::Int, s2::Int, decomp= "svd"; kwargs...)
Swapping two _neighboring_ sites.
""" ->
swapSites!(state::MPSState, s1::Int, s2::Int, decomp= "svd"; kwargs...)

@doc """
	oneShot(state::MPSState, sites::Vector{Int}; kwargs...)
Doing one measurement on state.

This is a helper function called by **measure!**. The state needs to be pretreated (localize qubits to be measured & canonicalize) before calling this funciton.

# Note: Internal function.
""" ->
oneShot(state::MPSState, sites::Vector{Int}; kwargs...)

@doc """
	collapse!(state::MPSState, qubit::Int; kwargs...)
Collapse the specified qubit and return the measurement result.

# Argument
- `state` : MPSState
- `qubit`: qubit to be collapsed
- `reset`: [_keyword argument_] Whether or not to reset measured qubit to |0>. Default value: false.

# Example
``` julia
julia> state = MPSState(1)
1-qubit MPSState

julia> apply!(state, H, 1)
1-qubit MPSState

julia> measure!(state, 1, 4)
4-element Array{Int64,1}:
 0
 0
 1
 1

julia> collapse!(state, 1)
0

julia> measure!(state, 1, 4)
4-element Array{Int64,1}:
 0
 0
 0
 0

```
""" ->
collapse!(state::MPSState, qubit::Int; kwargs...)

@doc """
	collapse!(state::MPSState, qubits::Vector{Int}; kwargs...)
Collapse the specified qubits and return the measurement result.
""" ->
collapse!(state::MPSState, qubits::Vector{Int}; kwargs...)

@doc """
	dag(state::MPSState)
Return complex conjugate of the MPSState.
""" ->
dag(state::MPSState)

@doc """
	logprobability(state::MPSState, config::Vector{Int})
Return log of the probability to get given configuration.
""" ->
logprobability(state::MPSState, config::Vector{Int})

@doc """
	showStructure(io::IO, state::MPSState)
Print structure of the MPSState.

# Example
```
julia> state = MPSState(3)
3-qubit MPSState

julia> showStructure(stdout, state)
MPSState
Site 1:IndexSet(Index[(1|id=414|Link,l=1), (2|id=485|Site,q=1)])
Site 2:IndexSet(Index[(1|id=414|Link,l=1), (1|id=969|Link,l=2), (2|id=707|Site,q=2)])
Site 3:IndexSet(Index[(1|id=969|Link,l=2), (2|id=363|Site,q=3)])
siteForQubit:[1, 2, 3]
qubitAtSite:[1, 2, 3]
llim = 0, rlim = 4
```
See also: **showData**
""" ->
showStructure(io::IO, state::MPSState)

@doc """
	showData(io::IO,state::MPSState)
Print both structure and stored tensor data.

#Example

```
julia> state = MPSState(3)
3-qubit MPSState

julia> showData(stdout, state)
--- MPS state: ---
Site 1:
ITensor ord=2 (1|id=570|Link,l=1) (2|id=57|Site,q=1)
Main.YQ.ITensors.Dense{Float64}
 1.0  0.0
Site 2:
ITensor ord=3 (1|id=570|Link,l=1) (1|id=598|Link,l=2) (2|id=555|Site,q=2)
Main.YQ.ITensors.Dense{Float64}
[:, :, 1] =
 1.0

[:, :, 2] =
 0.0
Site 3:
ITensor ord=2 (1|id=598|Link,l=2) (2|id=278|Site,q=3)
Main.YQ.ITensors.Dense{Float64}
 1.0  0.0
```

See also: **showStructure**
""" ->
showData(io::IO,state::MPSState)

@doc """
	ρ(state::MPSState, config::Vector{Int})
Return probability density of the given configuration.

# Example
```julia
julia> state = MPSState(5)
5-qubit MPSState

julia> apply!(state, H(1))
5-qubit MPSState

julia> ρ(state, [0,0,0,0,0])
0.7071067811865475 + 0.0im

julia> ρ(state, [1,0,0,0,0])
0.7071067811865475 + 0.0im
```
# Note:
Length of `config` must match number of qubits in the state.
""" ->
ρ(state::MPSState, config::Vector{Int})


text:
################################################################################
#
#   pglCosets.jl : Generating the cosets
#
################################################################################

export pglCosets
"""
    pglCosets(g::RingElem)

Compute the coset representatives of PGL(F_q²)/PGL(F_q).

This algorithm is due to Jincheng Zhuang and Qi Cheng **[3]**, see References in
the module documentation for more details.
"""
function pglCosets(g::RingElem)
    Fq2 = parent(g)
    d = div(degree(Fq2), 2)
    z = string("z", d)
    charac::Int = characteristic(Fq2)
    Fq = FiniteField(charac, d, z)[1]
    Fqx = nonZeroElements(Fq)
    img = findImg(Fq2, Fq)
    MM, piv = projectFindInv(Fq2, Fq)
    q::Int = length(Fq)
    R = zeros(Fq2, 5, q)
    MS = MatrixSpace(Fq2, 2, 2)
    M = MS()
    A = Array{typeof(M), 1}()

    # We compute a table with containing at most q+1-th root
    # for each element in F_q, these roots are taken in F_q²
    for y in Fq2
        t = y^(q+1)
        i = indexin(Fq, projectLinAlg(Fq, t, MM, piv))
        #        i = coeff(t, 0) + q*coeff(t, 1) + 1 # Is coeff(t, 1) always
        #        equal to zero ? (probably)

        for j in 1:5
            if R[j, i] == y
                break
            elseif R[j, i] == 0
                R[j, i] = y
                break
            end
        end
    end

    # STEP 1: We add cosets representatives of the form
    # [1, x]
    # [y, 1]
    # with x, y in F_q²\F_q and xy ≠ 1 and where x and y were chosen in a way
    # that we have to do only O(q³) steps instead of O(q⁴).
    M[1, 1], M[2, 2] = Fq2(1), Fq2(1)
    for v in Fqx
        for w in Fq2
            ι = inv(embed(Fq2, v, img))
            t = (ι*w^q)^(q+1)-ι
            i = indexin(Fq, projectLinAlg(Fq, t, MM, piv))
#            i = coeff(t, 0) + q*coeff(t, 1) + 1
            if i != 1
                for j in 1:5
                    r = R[j, i]
                    if r != 0
                        y = ι*w^q + r
                        x = -embed(Fq2, v, img)*y + w + w^q
                        if x*y != 1 && x^q != x && y^q != y
                            M[1, 2], M[2, 1] = x, y
                            push!(A, deepcopy(M))
                            break
                        end
                    else
                        break
                    end
                end
            else
                y = ι*w^q
                x = -embed(Fq2, v, img)*y + w + w^q
                if x*y != 1 && x^q != x && y^q != y
                    M[1, 2], M[2, 1] = x, y
                    push!(A, deepcopy(M))
                end
            end
        end
    end

    # STEP 2: We add cosets representatives of the form
    # [1,  a]
    # [g, bg]
    # with a, b in F_q*, g in F_q²\F_q, and a ≠ b
    M[2, 1] = g
    for a in Fqx
        for b in Fqx
            if a != b
                M[1, 2] = embed(Fq2, a, img)
                M[2, 2] = embed(Fq2, b, img)*g
                push!(A, deepcopy(M))
            end
        end
    end

    # STEP 3: We add cosets representatives of the form
    # [0,  1]
    # [g, ag]
    # with a in F_q
    M[1, 1] = Fq2(0)
    M[1, 2] = Fq2(1)
    for a in Fq
        M[2, 2] = embed(Fq2, a, img)*g
        push!(A, deepcopy(M))
    end

    # STEP 4: We add cosets representatives of the form
    # [     0,  1]
    # [1 + ag, bg]
    # with a, b in F_q
    for a in Fq
        for b in Fq
            M[2, 1] = 1 + embed(Fq2, a, img)*g
            M[2, 2] = embed(Fq2, b, img)*g
            push!(A, deepcopy(M))
        end
    end

    # STEP 5: We add cosets representatives of the form
    # [  1, 0]
    # [ ag, g]
    # with a in F_q
    M[1, 1] = Fq2(1)
    M[1, 2] = Fq2(0)
    M[2, 2] = g
    for a in Fq
        M[2, 1] = embed(Fq2, a, img)*g
        push!(A, deepcopy(M))
    end

    # STEP 6: We add cosets representatives of the form
    # [  1,      0]
    # [ ag, 1 + bg]
    # with a, b in F_q
    for a in Fq
        for b in Fq
            M[2, 1] = embed(Fq2, a, img)*g
            M[2, 2] = 1 + embed(Fq2, b, img)*g
            push!(A, deepcopy(M))
        end
    end

    return A
end


text:
fatal = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
any_errors = false

push!(LOAD_PATH, "../src")

using FsBert
using Base.Test

tests = [ "encode_decode.jl" ]

println("Running tests ...")

for test in tests
    try
        include(test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
    catch e
        any_errors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
        if fatal
            rethrow(e)
        elseif !quiet
            showerror(STDOUT, e, backtrace())
            println()
        end
    end
end

if any_errors
    throw("Tests failed")
end

stdin = joinpath(dirname(@__FILE__), "stdin.sh")
ENV2 = copy(ENV)
ENV2["JULIA_HOME"] = JULIA_HOME
run(setenv(`bash $stdin`, ENV2))

max_stars_repo_name: esl/FsBert

max_issues_repo_name: esl/FsBert

max_forks_repo_name: esl/FsBert


text:
"""
Bernoulli scheme is a generalization of Bernoulli proceses to dimension N where the random
variable may take on values with probability p₁,...,p_N that sum to 1.

In our case, we assume p₁,...,p_{N-1} are probabilities of different types of success,
p_N is the probability

The constructor takes probabilities p_1, ... p_N-1 and calculates p_N = 1-sum_i=1^N-1(pᵢ)
"""
struct NonhomogenousBernoulliScheme<:DiscreteCountingProcess
    p_matrix::Array{Float64, 2}

    function NonhomogenousBernoulliScheme(p_matrix)
        # @assert 0<sum(p_suc)<=1
        new(p_matrix)
    end
end

"""
Get parameters of the process
"""
Distributions.params(bs::NonhomogenousBernoulliScheme) = (bs.p_matrix)

"""
Get outcome distribution for given index.
Support for the outcome distribution is 1...N+1. N+1 means failure.
Index can be in the range 1, ..., n (n is the number of random variables in the scheme, not the number of outcomes).
"""
Base.getindex(bs::NonhomogenousBernoulliScheme, i::Integer) = 1<=i<=size(bs.p_matrix)[2] ? Categorical([bs.p_matrix[:, i]..., 1-sum(bs.p_matrix[:, i])]) : throw(BoundsError(bs, i))


text:
@marginalrule typeof(+)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, vout = mean_var(m_out)
    return (in1 = prod(ProdAnalytical(), NormalMeanVariance(mout - mean(m_in2), vout), m_in1), in2 = m_in2)
end

@marginalrule typeof(+)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily) = begin
    mout, vout = mean_var(m_out)
    return (in1 = m_in1, in2 = prod(ProdAnalytical(), NormalMeanVariance(mout - mean(m_in1), vout), m_in2))
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalMeanCovariance, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, Vout = mean_cov(m_out)
    return (in1 = prod(ProdAnalytical(), MvNormalMeanCovariance(mout - mean(m_in2), Vout), m_in1), in2 = m_in2)
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    mout, Wout = mean_precision(m_out)
    return (in1 = prod(ProdAnalytical(), MvNormalMeanPrecision(mout - mean(m_in2), Wout), m_in1), in2 = m_in2)
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalWeightedMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    xiout, Wout = weightedmean_precision(m_out)
    tmp = -Wout*mean(m_in2)
    tmp .+= xiout
    return (in1 = prod(ProdAnalytical(), MvNormalWeightedMeanPrecision(tmp, Wout), m_in1), in2 = m_in2)
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalMeanCovariance, m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    mout, Vout = mean_cov(m_out)
    return (in1 = m_in1, in2 = prod(ProdAnalytical(), MvNormalMeanCovariance(mout - mean(m_in1), Vout), m_in2))
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalMeanPrecision, m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    mout, Wout = mean_precision(m_out)
    return (in1 = m_in1, in2 = prod(ProdAnalytical(), MvNormalMeanPrecision(mout - mean(m_in1), Wout), m_in2))
end

@marginalrule typeof(+)(:in1_in2) (m_out::MvNormalWeightedMeanPrecision, m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    xiout, Wout = weightedmean_precision(m_out)
    tmp = -Wout * mean(m_in1)
    tmp .+= xiout
    return (in1 = m_in1, in2 = prod(ProdAnalytical(), MvNormalWeightedMeanPrecision(tmp, Wout), m_in2))
end

@marginalrule typeof(+)(:in1_in2) (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_in1, W_in1 = weightedmean_precision(m_in1)
    xi_in2, W_in2 = weightedmean_precision(m_in2)

    return MvNormalWeightedMeanPrecision([ xi_in1 + xi_out; xi_in2 + xi_out ], [ W_in1 + W_out W_out; W_out W_in2 + W_out])
end

@marginalrule typeof(+)(:in1_in2) (m_out::MultivariateNormalDistributionsFamily, m_in1::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_in1, W_in1 = weightedmean_precision(m_in1)
    xi_in2, W_in2 = weightedmean_precision(m_in2)

    T = promote_type(eltype(W_out), eltype(W_in1), eltype(W_in2))
    d = length(xi_out)
    Λ = Matrix{T}(undef, (2*d, 2*d))
    @inbounds for k2 in 1:d
        @inbounds for k1 in 1:d
            tmp = W_out[k1,k2]
            k1d = k1+d
            k2d = k2+d
            Λ[k1,k2] = tmp + W_in1[k1,k2]
            Λ[k1d,k2] = tmp
            Λ[k1,k2d] = tmp
            Λ[k1d,k2d] = tmp + W_in2[k1,k2]
        end
    end

    xi = Vector{T}(undef, 2*d)
    @inbounds for k = 1:d
        tmp = xi_out[k]
        xi[k]   = tmp + xi_in1[k]
        xi[k+d] = tmp + xi_in2[k]
    end

    # naive: return MvNormalWeightedMeanPrecision([ xi_in1 + xi_out; xi_in2 + xi_out ], [ W_in1 + W_out W_out; W_out W_in2 + W_out])
    return MvNormalWeightedMeanPrecision(xi, Λ)
end


text:
module TestQuery

import PublicAPI
using Test
import ..PublicAPITests

function test_of()
    apis = PublicAPI.of(PublicAPI)
    @test sort!(nameof.(apis)) == [Symbol("@public"), Symbol("@strict"), :of]
    @test Module.(apis) == fill(PublicAPI, 3)
end

function test_recursive()
    apis = PublicAPI.of(PublicAPITests)
    @test sort!(nameof.(apis)) == [:InplaceSamples, :ONE, :TWO]
end

end  # module


text:
function graph(T, m::Model)
    g = T()

    mvars = variables(m)
    for v in mvars
        add!(g, v)
    end

    for (v, expr) in pairs(m.vals)
        for p in variables(expr) ∩ variables(m)
            add!(g, p, v)
        end
    end

    for (v, expr) in pairs(m.dists)
        for p in variables(expr) ∩ variables(m)
            add!(g, p, v)
        end
    end

    g
end

export digraph
digraph(m::Model) = graph(SimpleDigraph{Symbol}, m)

export poset
poset(m::Model) = graph(SimplePoset{Symbol}, m)


text:
"""In-place version of `signed_exponent(::Array)`."""
function signed_exponent!(A::Array{T}) where {T<:Union{Float16,Float32,Float64}}

    # sign&fraction mask
    sfmask = Base.sign_mask(T) | Base.significand_mask(T)
    emask = Base.exponent_mask(T)

    sbits = Base.significand_bits(T)
    bias  = Base.exponent_bias(T)
    ebits = Base.exponent_bits(T)-1

    for i in eachindex(A)
        ui = reinterpret(Unsigned,A[i])
        sf = ui & sfmask                    # sign & fraction bits
        e = ((ui & emask) >> sbits) - bias  # de-biased exponent
        eabs = e == -bias ? 0 : abs(e)      # for iszero(A[i]) e == -bias, set to 0
        esign = (e < 0 ? 1 : 0) << ebits    # determine sign of exponent
        esigned = ((esign | eabs) % typeof(ui)) << sbits    # concatentate exponent

        A[i] = reinterpret(T,sf | esigned)  # concatenate everything back together
    end
end

"""
```julia
B = signed_exponent(A::Array{T}) where {T<:Union{Float16,Float32,Float64}}
```
Converts the exponent bits of Float16,Float32 or Float64-arrays from its
conventional biased-form into a sign&magnitude representation.

# Example

```julia
julia> bitstring(10f0,:split)
"0 10000010 01000000000000000000000"

julia> bitstring.(signed_exponent([10f0]),:split)[1]
"0 00000011 01000000000000000000000"
```

In the former the exponent 3 is interpret from 0b10000010=130 via subtraction of
the exponent bias of Float32 = 127. In the latter the exponent is inferred from
sign bit (0) and a magnitude represetation 2^1 + 2^1 = 3."""
function signed_exponent(A::Array{T}) where {T<:Union{Float16,Float32,Float64}}
    B = copy(A)
    signed_exponent!(B)
    return B
end


text:
"""
    Berlage(; <keyword arguments>)

Create a Berlage wavelet.

# Arguments
- `dt=0.002`: sampling interval in secs.
- `f0=20.0`: central frequency in Hz.
- `m::Real=2`: exponential parameter of Berlage wavelet.
- `alpha::Real=180.0`: alpha parameter of Berlage wavelet in rad/secs.
- `phi0::Real`: phase rotation in radians.

# Example
```julia
julia> w = Berlage(); plot(w);
```
**Reference**
* Aldridge, David F., 1990, The berlage wavelet: GEOPHYSICS, 55, 1508--1511.
"""
function Berlage(; dt=0.002, f0=20.0, m::Real=2, alpha::Real=180.0,
                 phi0::Real=0.0)

    nw = floor(Int, 2.2/(f0*dt))
    t = dt*collect(0:1:nw-1)
    w = (t .^ m).* exp.(-alpha*t) .* cos.(2*pi*f0*t .+ phi0);
    w = w/maximum(w)

end


text:
###############################################################################
#
# Convert energy to photons
#
###############################################################################
_FAC(FT) = FT(1e-9) / (H_PLANCK(FT) * LIGHT_SPEED(FT) * AVOGADRO(FT))

"""
    e2phot(λ::Array{FT}, E::Array{FT}) where {FT<:AbstractFloat}

Calculates the number of moles of photons, given
- `λ` An array of wave length in `[nm]`, converted to `[m]` by _FAC
- `E` Joules of energy
"""
function e2phot(
            λ::Array{FT},
            E::Array{FT}
) where {FT<:AbstractFloat}
    return (E .* λ) .* _FAC(FT)
end

"""
    e2phot!(
                λ::Array{FT,1},
                E::Array{FT,1},
                cache::Array{FT,1}
    ) where {FT<:AbstractFloat}

Calculates the number of moles of photons, given
- `λ` An array of wave length in `[nm]`, converted to `[m]` by _FAC
- `E` Joules of energy
- `cache` Cache to avoid memory allocations
"""
function e2phot!(
            λ::Array{FT,1},
            E::Array{FT,1},
            cache::Array{FT,1}
) where {FT<:AbstractFloat}
    cache .= (E .* λ) .* _FAC(FT);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using SFML

# Create the window with a width of 800 and a height of 600
window = RenderWindow("Camera/Views", 800, 600)
# Always set the framerate limit
set_framerate_limit(window, 60)
# This is for listening for events
event = Event()

# Make a circle shape
circle = CircleShape()
# Do all the setup
set_radius(circle, 40)
set_position(circle, Vector2f(500, 200))
set_fillcolor(circle, SFML.red)

# Make a square
square = RectangleShape()
set_size(square, Vector2f(60, 60))
set_position(square, Vector2f(200, 400))
set_fillcolor(square, SFML.cyan)

# Cameras are views in SFML
# Pass it center of the view, and size of the view (same as window in this case)
view = View(Vector2f(400, 300), Vector2f(800, 600))

while isopen(window)
    # Check for any events
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            # Close the window if a closed event is detected
            close(window)
        end
    end

    # Check keypresses
    if is_key_pressed(KeyCode.LEFT)
        # Move left
        move(view, Vector2f(-0.5, 0))
    end
    if is_key_pressed(KeyCode.RIGHT)
        # Move right
        move(view, Vector2f(0.5, 0))
    end
    if is_key_pressed(KeyCode.UP)
        # Move up
        move(view, Vector2f(0, -0.5))
    end
    if is_key_pressed(KeyCode.DOWN)
        # Move down
        move(view, Vector2f(0, 0.5))
    end
    # Zoom out
    if is_key_pressed(KeyCode.Z)
        zoom(view, 0.5)
        set_size(view, Vector2f(1600, 1200))
    end
    # Zoom in
    if is_key_pressed(KeyCode.X)
        zoom(view, 1.5)
        set_size(view, Vector2f(400, 300))
    end

    # Attach the view to the window
    set_view(window, view)

    # Clear and redraw the screen
    clear(window, SFML.white)
    draw(window, square)
    draw(window, circle)
    display(window)
end


text:
# Test the mut_evolve_exact() function from src/Evolve.jl

#include("../src/CGP.jl")
using Main.CGP
const MyInt = UInt8

goal1 = Main.CGP.MyInt[ 0xb, 0xe ]
goal2 = Main.CGP.MyInt[ 0x6, 0x9 ]
goal3 = Main.CGP.MyInt[ 0xa, 0xe ]
#goallist = [goal1, goal2,goal3]

goallist = [[0x0,0x1],[0x07, 0x01],[0x0d, 0x0e],[0x07, 0x0d],[0xb,0xe],[0x6,0x9],[0xa,0xe]]

numinputs = 2
funcs = default_funcs(numinputs)
p = Parameters( numinputs=numinputs, numoutputs=2, numinteriors=6, numlevelsback=6 )
context = construct_context(numinputs)
c = random_chromosome( p, funcs )
max_steps = 20
result = mut_evolve_exact( random_chromosome( p, funcs ), context, goallist, max_steps )

Random.seed!(5);
c = random_chromosome( p, funcs )
ec = eval_chromosome(ch5,goallist)
println("ec: ",ec)
@assert ec[1] == [0, 0, 2, 0, 1, 0, 1]
@assert chrome_check(c,goallist,0) == 1
@assert chrome_check(c,goallist,1) == 1
@assert chrome_check(c,goallist,2) == 0
@assert chrome_check(c,goallist,3) == -1


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Box(min, max)

An axis-aligned box with `min` and `max` corners.
See https://en.wikipedia.org/wiki/Hyperrectangle.

## Example

```julia
Box(Point(0,0,0), Point(1,1,1)) # unit cube
```
"""
struct Box{Dim,T} <: Primitive{Dim,T}
  min::Point{Dim,T}
  max::Point{Dim,T}
end

Box(min::Tuple, max::Tuple) = Box(Point(min), Point(max))

paramdim(::Type{<:Box{Dim}}) where {Dim} = Dim

Base.minimum(b::Box) = b.min
Base.maximum(b::Box) = b.max
Base.extrema(b::Box) = b.min, b.max
center(b::Box) = Point((coordinates(b.max) + coordinates(b.min)) / 2)
measure(b::Box) = prod(b.max - b.min)
diagonal(b::Box) = norm(b.max - b.min)
sides(b::Box) = b.max - b.min

function Base.in(p::Point{Dim}, b::Box{Dim}) where {Dim}
  l, u = coordinates.((b.min, b.max))
  x = coordinates(p)
  for i in 1:Dim
    l[i] ≤ x[i] && x[i] ≤ u[i] || return false
  end
  true
end


text:
# Autogenerated wrapper script for TetGen_jll for aarch64-apple-darwin
export libtet

JLLWrappers.@generate_wrapper_header("TetGen")
JLLWrappers.@declare_library_product(libtet, "@rpath/libtet.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtet,
        "lib/libtet.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ExcelReaders

function getdice2013excelparameters(filename)
    p = Dict{Symbol,Any}()

    T = 60

    #Open RICE_2010 Excel File to Read Parameters
    f = openxl(filename)

    p[:a1]          = getparams(f, "B25:B25", :single, "Base", 1)       #Damage coefficient on temperature
    p[:a2]          = getparams(f, "B26:B26", :single, "Base", 1)       #Damage quadratic term
    p[:a3]          = getparams(f, "B27:B27", :single, "Base", 1)       #Damage exponent
    p[:al]          = getparams(f, "B21:BI21", :all, "Base", T)         #Level of total factor productivity
    p[:b11]         = getparams(f, "B65:B65", :single, "Base", 1)/100   #Carbon cycle transition matrix atmosphere to atmosphere
    p[:b12]         = getparams(f, "B67:B67", :single, "Base", 1)/100   #Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b21]         = getparams(f, "B66:B66", :single, "Base", 1)/100   #Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
    p[:b22]         = getparams(f, "B68:B68", :single, "Base", 1)/100   #Carbon cycle transition matrix shallow ocean to shallow oceans
    p[:b23]         = getparams(f, "B70:B70", :single, "Base", 1)/100   #Carbon cycle transition matrix shallow to deep ocean
    p[:b32]         = getparams(f, "B69:B69", :single, "Base", 1)/100   #Carbon cycle transition matrix deep ocean to shallow ocean
    p[:b33]         = getparams(f, "B71:B71", :single, "Base", 1)/100   #Carbon cycle transition matrix deep ocean to deep oceans
    p[:c1]          = getparams(f, "B82:B82", :single, "Base", 1)       #Speed of adjustment parameter for atmospheric temperature (per 5 years)
    p[:c3]          = getparams(f, "B83:B83", :single, "Base", 1)       #Coefficient of heat loss from atmosphere to oceans
    p[:c4]          = getparams(f, "B84:B84", :single, "Base", 1)       #Coefficient of heat gain by deep oceans
    p[:cca0]        = getparams(f, "B92:B92", :single, "Base", T)       #Initial cumulative industrial emissions
    p[:cost1]       = getparams(f, "B32:BI32", :all, "Base", T)         #Abatement cost function coefficient
    p[:damadj]      = getparams(f, "B56:B56", :single, "Parameters", 1) #Adjustment exponent in damage function
    p[:dk]          = getparams(f, "B6:B6", :single, "Base", 1)         #Depreciation rate on capital (per year)
    p[:elasmu]      = getparams(f, "B19:B19", :single, "Base", 1)       #Elasticity of MU of consumption
    p[:eqmat]       = getparams(f, "B71:B71", :single, "Parameters", 1) #Equilibirum concentration of CO2 in atmosphere (GTC)
    p[:etree]       = getparams(f, "B44:BI44", :all, "Base", T)         #Carbon emissions form land use change (GTCO2 per year)
    p[:expcost2]    = getparams(f, "B39:B39", :single, "Base", 1)       #Exponent of control cost function
    p[:fco22x]      = getparams(f, "B80:B80", :single, "Base", 1)       #Forcings of equilibrium CO2 doubling (Wm-2)
    p[:forcoth]     = getparams(f, "B73:BI73", :all, "Base", T)         #Exogenous forcing for other greenhouse gases
    p[:fosslim]     = getparams(f, "B57:B57", :single, "Base", 1)       #Maximum carbon resources (Gtc)
    p[:gama]        = getparams(f, "B5:B5", :single, "Base", 1)         #Capitail Share
    p[:k0]          = getparams(f, "B12:B12", :single, "Base", 1)       #Initial capital
    p[:l]           = getparams(f, "B53:BI53", :all, "Base", T)         #Level of population and labor (millions)
    p[:mat0]        = getparams(f, "B60:B60", :single, "Base", 1)       #Initial Concentration in atmosphere at end 2007 (GtC)
    p[:MIU]         = getparams(f, "B135:BI135", :all, "Base", T)       #Optimized emission control rate results from DICE2013 (base case)
    p[:ml0]         = getparams(f, "B63:B63", :single, "Base", 1)       #Initial Concentration in deep oceans 2008 (GtC)
    p[:mu0]         = getparams(f, "B62:B62", :single, "Base", 1)       #Initial Concentration in biosphere/shallow oceans 2008 (GtC)
    p[:partfract]   = getparams(f, "B47:BI47", :all, "Base", T)         #Fraction of emissions in control regime
    p[:pbacktime]   = getparams(f, "B34:BI34", :all, "Base", T)         #Backstop price ($1000 per ton CO2)
    p[:rr]          = getparams(f, "B18:BI18", :all, "Base", T)         #Social Time Preference Factor
    p[:S]           = getparams(f, "B131:BI131", :all, "Base", T)       #Optimized savings rate (fraction of gross output) results from DICE2013 (base case)
    p[:scale1]      = getparams(f, "B49:B49", :single, "Base", 1)       #Multiplicative scaling coefficient
    p[:scale2]      = getparams(f, "B50:B50", :single, "Base", 1)       #Additive scaling coefficient
    p[:sigma]       = getparams(f, "B41:BI41", :all, "Base", T)         #(industrial, MTCO2/$1000    2000 US$)
    p[:t2xco2]      = getparams(f, "B79:B79", :single, "Base", 1)       #Equilibrium temp impact (oC per doubling CO2)
    p[:tatm0]       = getparams(f, "B76:B76", :single, "Base", 1)       #Initial atmospheric temp change 2008-2011 (C from 1900)
    p[:tocean0]     = getparams(f, "B77:B77", :single, "Base", 1)       #Initial temperature of deep oceans (deg C aboce 1900)
    p[:usedamadj]   = true

    return p
end

function getdice2013gamsparameters(filename)
    p = Dict{Symbol,Any}()

    T = 60

    #Open RICE_2010 Excel File to Read Parameters
    f = openxl(filename)
    sheet = "DICE2013_Base"

    p[:a1]          = getparams(f, "B41:B41", :single, sheet, 1)       #Damage coefficient on temperature
    p[:a2]          = getparams(f, "B42:B42", :single, sheet, 1)       #Damage quadratic term
    p[:a3]          = getparams(f, "B43:B43", :single, sheet, 1)       #Damage exponent
    p[:al]          = getparams(f, "B5:BI5", :all, sheet, T)           #Level of total factor productivity
    p[:b11]         = getparams(f, "B22:B22", :single, sheet, 1)       #Carbon cycle transition matrix atmosphere to atmosphere
    p[:b12]         = getparams(f, "B20:B20", :single, sheet, 1)       #Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b21]         = getparams(f, "B23:B23", :single, sheet, 1)       #Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
    p[:b22]         = getparams(f, "B24:B24", :single, sheet, 1)       #Carbon cycle transition matrix shallow ocean to shallow oceans
    p[:b23]         = getparams(f, "B21:B21", :single, sheet, 1)       #Carbon cycle transition matrix shallow to deep ocean
    p[:b32]         = getparams(f, "B25:B25", :single, sheet, 1)       #Carbon cycle transition matrix deep ocean to shallow ocean
    p[:b33]         = getparams(f, "B26:B26", :single, sheet, 1)       #Carbon cycle transition matrix deep ocean to deep oceans
    p[:c1]          = getparams(f, "B36:B36", :single, sheet, 1)       #Speed of adjustment parameter for atmospheric temperature (per 5 years)
    p[:c3]          = getparams(f, "B37:B37", :single, sheet, 1)       #Coefficient of heat loss from atmosphere to oceans
    p[:c4]          = getparams(f, "B38:B38", :single, sheet, 1)       #Coefficient of heat gain by deep oceans
    p[:cca0]        = getparams(f, "B14:B14", :single, sheet, 1)       #Initial cumulative industrial emissions
    p[:cost1]       = getparams(f, "B46:BI46", :all, sheet, T)         #Abatement cost function coefficient
    p[:damadj]      = getparams(f, )                                   #Adjustment exponent in damage function
    p[:dk]          = getparams(f, "B8:B8", :single, sheet, 1)         #Depreciation rate on capital (per year)
    p[:elasmu]      = getparams(f, "B54:BI54", :single, sheet, 1)      #Elasticity of MU of consumption
    p[:eqmat]       = getparams(f, 588.)                               #Equilibirum concentration of CO2 in atmosphere (GTC)
    p[:etree]       = getparams(f, "B13:BI13", :all, sheet, T)         #Carbon emissions form land use change (GTCO2 per year)
    p[:expcost2]    = getparams(f, "B48:B48", :single, sheet, 1)       #Exponent of control cost function
    p[:fco22x]      = getparams(f, "B30:B30", :single, sheet, 1)       #Forcings of equilibrium CO2 doubling (Wm-2)
    p[:forcoth]     = getparams(f, "B29:BI29", :all, sheet, T)         #Exogenous forcing for other greenhouse gases
    p[:gama]        = getparams(f, "B7:B7", :single, sheet, 1)         #Capitail Share
    p[:k0]          = getparams(f, "B9:B9", :single, sheet, 1)         #Initial capital
    p[:l]           = getparams(f, "B6:BI6", :all, sheet, T)           #Level of population and labor (millions)
    p[:mat0]        = getparams(f, "B17:B17", :single, sheet, 1)       #Initial Concentration in atmosphere at end 2007 (GtC)
    p[:MIU]         = getparams(f, "B47:BI47", :all, sheet, T)         #Optimized emission control rate results from DICE2013 (base case)
    p[:ml0]         = getparams(f, "B19:B19", :single, sheet, 1)       #Initial Concentration in deep oceans 2008 (GtC)
    p[:mu0]         = getparams(f, "B18:B18", :single, sheet, 1)       #Initial Concentration in biosphere/shallow oceans 2008 (GtC)
    p[:partfract]   = getparams(f, "B49:BI49", :all, sheet, T)         #Fraction of emissions in control regime
    p[:pbacktime]   = getparams(f, "B50:BI50", :all, sheet, T)         #Backstop price ($1000 per ton CO2)
    p[:rr]          = getparams(f, "B55:BI55", :all, sheet, T)         #Social Time Preference Factor
    p[:S]           = getparams(f, "B51:BI51", :all, sheet, T)         #Optimized savings rate (fraction of gross output) results from DICE2013 (base case)
    p[:scale1]      = getparams(f, "B56:B56", :single, sheet, 1)       #Multiplicative scaling coefficient
    p[:scale2]      = getparams(f, "B57:B57", :single, sheet, 1)       #Additive scaling coefficient
    p[:sigma]       = getparams(f, "B12:BI12", :all, sheet, T)         #(industrial, MTCO2/$1000    2000 US$)
    p[:t2xco2]      = getparams(f, "B33:B33", :single, sheet, 1)       #Equilibrium temp impact (oC per doubling CO2)
    p[:tatm0]       = getparams(f, "B34:B34", :single, sheet, 1)       #Initial atmospheric temp change 2008-2011 (C from 1900)
    p[:tocean0]     = getparams(f, "B35:B35", :single, sheet, 1)       #Initial temperature of deep oceans (deg C aboce 1900)
    p[:usedamadj]   = false

    return p
end


text:
"""
# Available lower level routines
* [`besseljs3d`](@ref) spherical Bessel function eval
"""
function lower_level_routs() end

"""
    function fj, fjder = besseljs3d(nterms,z;scale=1.0,ifder=0)
This subroutine evaluates the first `nterms` spherical Bessel
functions, and if requested, their derivatives.
It incorporates a scaling parameter `scale` so that

      	fjs_n(z)=j_n(z)/SCALE^n
      	fjder_n(z)=\\frac{\\partial fjs_n(z)}{\\partial z}
# Input
* `nterms::Integer` order of expansion of output array `fjs`
* `z::ComplexF64` argument of the spherical Bessel functions
* `scale::Float64` scaling factor
* `ifder::Integer1` flag indicating whether to calculate `fjder`
      	          0	NO
      	          1	YES
OUTPUT:
* `fjs::Array{ComplexF64}` array of length `nterms+1` of scaled Bessel functions.
* `fjder::Array{ComplexF64}` array of derivatives of scaled Bessel functions, if requested.
"""
function besseljs3d(nterms::Integer,z::ComplexF64;scale::Float64=1.0,
                    ifder::Integer=0)

    @assert (nterms >= 0)
    if (ifder !=0 && ifder != 1)
        @warn "unexpected value in ifder, no ders computed"
    end

    fjs = Array{ComplexF64}(undef,nterms+1)
    fjder = ComplexF64(0)

    if ifder == 1
        fjder = Array{ComplexF64}(undef,nterms+1)
    end

    # fortran interface
    # subroutine besseljs3d(nterms,z,scale,fjs,ifder,fjder)

    ccall((:besseljs3d_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,Fc),
          nterms,z,scale,fjs,ifder,fjder)

    if ifder != 1
        fjder = nothing
    end

    return fjs, fjder

end


text:
module GeneralWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractAnalyticalPropagator = @jimport org.orekit.propagation.analytical.AbstractAnalyticalPropagator
const AbstractPropagator = @jimport org.orekit.propagation.AbstractPropagator
const AdditionalStateProvider = @jimport org.orekit.propagation.AdditionalStateProvider
const Appendable = @jimport java.lang.Appendable
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BoundedPropagator = @jimport org.orekit.propagation.BoundedPropagator
const BufferedReader = @jimport java.io.BufferedReader
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const EphemerisFile = @jimport org.orekit.files.general.EphemerisFile
const EphemerisFileParser = @jimport org.orekit.files.general.EphemerisFileParser
const EphemerisFileWriter = @jimport org.orekit.files.general.EphemerisFileWriter
const EphemerisSegmentPropagator = @jimport org.orekit.files.general.EphemerisSegmentPropagator
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const Frame = @jimport org.orekit.frames.Frame
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const Object = @jimport java.lang.Object
const OrekitEphemerisFile = @jimport org.orekit.files.general.OrekitEphemerisFile
const OrekitEphemerisFile_OrekitSatelliteEphemeris = @jimport org.orekit.files.general.OrekitEphemerisFile$OrekitSatelliteEphemeris
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates

export EphemerisFile
export EphemerisFileParser
export EphemerisFileWriter
export EphemerisSegmentPropagator
export OrekitEphemerisFile
export add_additional_state_provider
export add_event_detector
export add_satellite
export clear_events_detectors
export equals
export get_additional_state_providers
export get_attitude_provider
export get_class
export get_events_detectors
export get_frame
export get_generated_ephemeris
export get_initial_state
export get_managed_additional_states
export get_max_date
export get_min_date
export get_mode
export get_pv_coordinates
export get_pv_provider
export get_satellites
export hash_code
export is_additional_state_managed
export notify
export notify_all
export parse
export propagate
export reset_initial_state
export set_attitude_provider
export set_ephemeris_mode
export set_master_mode
export set_slave_mode
export to_string
export wait
export write

include("ephemeris_file.jl")
include("ephemeris_file_parser.jl")
include("ephemeris_file_writer.jl")
include("ephemeris_segment_propagator.jl")
include("orekit_ephemeris_file.jl")

end


text:
abstract Type1FuzzySet <: AbstractFuzzySet
export Type1FuzzySet

#####################
abstract Membershipfunction <: Type1FuzzySet
export Membershipfunction

import Base: getindex

getindex(mf::Membershipfunction, range) = [mf[value]::Float64 for value in range]
getindex(mf::Membershipfunction, values...) = [mf[value]::Float64 for value in values]

export getindex

#####################
# Membership Functions
"Triangle Membership Function"
immutable TriangleMF <: Membershipfunction
	a::Float64
	b::Float64
	c::Float64

	TriangleMF(a,b,c) = a <= b <= c ? new(a,b,c) : error("Indices out of order")
end
export TriangleMF

import Base: getindex
function getindex(trimf::TriangleMF, value::Real)
	(value <= trimf.a || trimf.c <= value) && return zero(Float64)
	(value == trimf.b) && return one(Float64)

	if trimf.a < value < trimf.b
		slope = (1.0 - 0.0)/(trimf.b - trimf.a)
		return slope * (value - trimf.a)
	else
		slope = (0.0 - 1.0)/(trimf.c - trimf.b)
		return (slope * (value - trimf.b) + 1.0)
	end
end
export getindex

"Trapezoid Membership Function"
immutable TrapezoidMF <: Membershipfunction
	a::Float64
	b::Float64
	c::Float64
	d::Float64
end
export TrapezoidMF

import Base: getindex
function getindex(trapmf::TrapezoidMF, value::Real)
	(value <= trapmf.a || trapmf.d <= value) && return zero(Float64)
	(trapmf.b <= value <= trapmf.c) && return one(Float64)

	if trapmf.a < value < trapmf.b
		slope = (1 - 0)/(trapmf.b - trapmf.a)
		return slope * (value - trapmf.a)
	else
		slope = (0 - 1)/(trapmf.d - trapmf.c)
		return slope * (value - trapmf.c) + 1
	end
end
export getindex

"Generic Membership Function"
immutable GenericMF <: Membershipfunction
	mf::Function
end
export GenericMF

import Base: getindex
function getindex(gmf::GenericMF, value::Real)
	gmf.mf(value)
end
export getindex

immutable GaussMF <: Membershipfunction
	mu::Float64
	sigma::Float64
end
export GaussMF
immutable GaussMF <: Membershipfunction
	mu::Float64
	sigma::Float64
end
export GaussMF

import Base: getindex
function getindex(gmf::GaussMF, value::Real)
	GenericMF(x -> exp(-(x-gmf.mu).^2/(2.*gmf.sigma.^2)))[value]
end
export getindex

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


text:
## --------------  Debugging

const dbgLow = true
const dbgMedium = true
const dbgHigh = true

## ------------  Symbols

# Name of the constant term in regressions
const RegrIntercept = :cons

# ------------