text:
"""
    Household

A module to solve a standard household life-cycle problem stated in:
[page](https://clpr.github.io/pages/blogs/190622_The_household_life_cycle_problem_with_a_continuous_medical_expenditure_shock.html)

depends on our custom module `src`
"""
module Household
    import src # custom module, also requires Distributions package; used in dynamic programming

    export LifeCycleDatPkg, solve, get_a, solve_dp

# --------------- CRRA UTILITY
"""
    u_CRRA(c::Real ; γ::Real = 1.5, ψ::Real = 1E-12)

CRRA utility `` u(c) = [ (c+\\psi)^{1-\\gamma} - 1 ] / (1 - \\gamma) ``,
where `ψ` is a minor amount to avoid absolute zero, `c` must be non-negative and finite,
and `γ` is relative risk aversion. It must be greater than 0.
When `γ` is 1, CRRA degenerates to logarithm utility ``u(c) = \\ln(c+\\psi)``.
This function returns a `Real` value of utility.

Timing information:
```
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.201 ns (0.00% GC)
  median time:      1.202 ns (0.00% GC)
  mean time:        1.478 ns (0.00% GC)
  maximum time:     50.175 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
```
"""
function u_CRRA(c::Real ; γ::Real = 1.5, ψ::Real = 1E-6)
    # assertion
    @assert( 0 <= c < Inf, "c must be non-negative and finite" )
    @assert( 0 <= γ < Inf, "γ must be non-negative and finite" )
    @assert( ψ >= 0, "ψ must be non-negative" )
    # branch
    γ == 1 ? (return log(c + ψ)::Real) : nothing
    return ( ((c+ψ)^(1-γ) - 1) / (1-γ) )::Real
end # u_CRRA

# ---------------- Types
# In this section, we define a structure to define a standard household problem.
# It helps standard I/O stream
struct LifeCycleDatPkg
    S::Int # length of life
    a1::Real # initial asset
    β::Real # utility discounting factor
    γ::Real # risk aversion of CRRA utility
    # budget: ``A_{s}a_{s+1} = B_{s}a_{s} - E_{s}c_{s} + F_{s},s=1,\\dots,S-1 ``
    A::Vector{T} where T <: Real # coef on a_{s+1}
    B::Vector{T} where T <: Real # coef on a_{s}
    E::Vector{T} where T <: Real # coef on c_{s} (expectation, if there are random components)
    F::Vector{T} where T <: Real # constant term
    # --------------------
    function LifeCycleDatPkg(A::Vector, B::Vector, E::Vector, F::Vector ; a1::Real = 0, β::Real = 0.99, γ::Real = 1.5, S::Int = length(A))
        # domain check
        @assert(0<S<Inf, "S must be positive and finite integer")
        @assert(0<= a1 <Inf, "a1 must be non-negative and finite")
        @assert(0<β<Inf, "beta must be positive and finite")
        @assert(0<γ<Inf, "gamma must be positive and finite")
        @assert(all(A .> 0.0), "A_{s} must be greater than 0")
        @assert(all(B .> 0.0), "B_{s} must be greater than 0")
        @assert(all(E .> 0.0), "The expectations of E_{s} must be greater than 0")
        @assert(all(isfinite.(F)), "F_{s} must be finite")
        # length check
        @assert(S==length(A)==length(B)==length(F)==length(E), "uncompatible vector length(s)")
        return new(S,a1,β,γ, A,B,E,F)
    end # LifeCycleDatPkg
end

# ---------------- testing dataset
demodat = LifeCycleDatPkg(
    cat(fill(0.99,39),[1.0,],dims=1), # As
    fill(1.05,40), # Bs
    fill(1.19,40), # E
    cat(fill(0.7,20),zeros(20),dims=1), # Fs
    a1 = 0.0, β = 0.99, γ = 1.5, S = 40
)

# ------------------------------------------------------------------------------
# NOTE: the following section is for ANALYTICAL solution WITHOUT borrowing constraints
# ---------------- Analytical Solver (no borrowing constraint)
"""
    solve(m::LifeCycleDatPkg)

solves a standard deterministic household life-cycle problem, using given data package.
returns a `NamedTuple` of two elements (in order): `a`, `c`;
where `a` is the expecation path of every period's asset,
and `c` is the deterministic decision path (controller) of consumption.
"""
function solve(m::LifeCycleDatPkg)
    # case 1: S == 1 (only one period)
    # NOTE: in this case, household just spend all their asset
    #       they follow the budget: `` 0 = A * a1 - E c1 + F ``
    if m.S == 1
        return ( a = Float64[m.a1,], c = Float64[ (m.A[1] * m.a1 + m.F[1]) / m.E ,] )::NamedTuple
    end # if

    # case 2: S >= 2 (at least two periods)
    # NOTE: in this case, we use general solution
    # 0. define \\tilde{\\prod} B_i/A_i (len= S)
    local tmpB2A = m.B ./ m.A
    local tildeBA = Float64[]
    for x in 2:m.S
        push!(tildeBA, prod(tmpB2A[x:m.S]) )
    end # for x
    push!(tildeBA,1.0)

    # 1. define d𝔼G/dc, the derivatives of compressed budgeting G on c_{s} (len= S)
    local d𝔼Gdc::Vector = tildeBA .* m.E ./ m.A
    # 2. define H̄, the Euler equation multiplier (len= S-1), H̄[s]: c[s] -> c[s+1]
    local H̄::Vector = d𝔼Gdc[1:(m.S-1)] ./ d𝔼Gdc[2:m.S]; H̄ .*= m.β; H̄ .^= -m.γ
    # 3. define M̄, the "cumulative" Euler equation multiplier (len= S), M̄[s]: c[1] -> c[s]
    local M̄::Vector = cumprod(H̄); insert!(M̄, 1, 1.0)
    println(M̄)
    # 4. define Ȳ, the denominator
    local Ȳ::Real = sum( m.E .* M̄ ./ m.A .* tildeBA )
    # 5. define X̄, the numerator
    local X::Real = m.a1 .* prod(m.B ./ m.A) .+ sum( m.F ./ m.A .* tildeBA )

    # check if X is valid
    X >= 0 ? nothing : throw(DomainError("negative consumption found"))

    # if alright, go on to get c[1], and extend the path
    local c::Vector{Float64} = (X/Ȳ) .* M̄
    # now, get the path/distribution of asset expectation in each period
    local a::Vector = Float64[m.a1,]
    for x in 1:(m.S-1); push!(a, (m.B[x]*a[x] - m.E[x]*c[x] + m.F[x])/m.A[x] ); end
    # finally, return a tuple
    return (a = a, c = c)::NamedTuple
end # solve

# ----------------- compute a[s], s=2 to S+1 when given c[s], s=1 to S
"""
    get_a(d::LifeCycleDatPkg, C::Vector)

when `s` is 1, `d.a1` will be directly returned;
otherwise, `length(C) >= s-1` requried
"""
function get_a(d::LifeCycleDatPkg, s::Int, C::Vector)
    if s==1; return (d.a1)::Real; end
    # if not a[1]
    @assert(length(C) >= (s-1), "length(C) >= s-1 requried")
    # computing
    local as::Real = d.a1
    for x in 1:(s-1)
        as = ( d.B[x] * as - d.E[x] * C[x] + d.F[x] ) / d.A[x]
    end # for x
    return as::Real
end # get_a

# =============== STD DP SOLVER (DETERMINISTIC)
# NOTE: using the same datastructure LifeCycleDatPkg

# --------------- small type to pass parameters
"""
    BellmanParSet

a simple structure to pass in the parameters required by `bellman_val()`
"""
struct BellmanParSet <: Any
    Athis::Real; Bthis::Real; Ethis::Real; Fthis::Real # pars in budget, used to compute c[s]
    β::Real # utility discounting factor
    γ::Real # risk aversion
    # ------- define a constructor to allow any-order parameter list
    function BellmanParSet(; Athis::Real = 0.0, Bthis::Real = 0.0, Ethis::Real = 0.0, Fthis::Real = 0.0, β::Real = 1.0, γ::Real = 2.0)
        return new(Athis,Bthis,Ethis,Fthis,β,γ)
    end
end # BellmanParSet
# --------------- a convenient reloaded external constructor
BellmanParSet(m::LifeCycleDatPkg, s::Int) = begin
    return BellmanParSet(Athis = m.A[s], Bthis = m.B[s], Ethis = m.E[s], Fthis = m.F[s], β = m.β, γ = m.γ)
end

# -------------- minor function: c[s] = [ B[s]a[s] + F[s] - A[s]a[s+1] ] / E[s]
"""
    get_cthis( athis::Real, anext::Real, pars::BellmanParSet )

computes c[s] when given a[s], a[s+1] and other parameters in the inter-temporal budget constraint.
returns a `Float64`
"""
get_cthis( athis::Real, anext::Real, pars::BellmanParSet ) = begin
    return ( (pars.Bthis * athis + pars.Fthis - pars.Athis * anext) / pars.Ethis )::Float64
end # get_cthis
# -------------- function of the value of value functions
"""
    bellman_val( athis::Real, anext::Real, anext_vec::Vector, vnext_vec::Vector, pars::BellmanParSet )

computes `u(c[s]) + βV_{s+1}(a[s+1])` (no `max` yet!) when given a[s] (`athis`),
a[s+1] (`anext`), the INCREASINGLY-SORTED discrete function series ``(a_{s+1},V_{s+1})`` (`anext_vec` & `vnext_vec`),
and other parameters (passed in `pars`) used in inter-temporal budgeting.
returns a `Float64`. If computed `c[s]` is infeasible (<0), this function returns `-9.99E99`

Please note that, `next_vec` & `vnext_vec` actually together define the function ``V_{s+1}(a_{s+1})``,
as independent variable and function value. They are in the form of discrete 2-dim points.
"""
function bellman_val( athis::Real, anext::Real, anext_vec::Vector, vnext_vec::Vector, pars::BellmanParSet )
    # get consumption c[s]
    local cthis::Real = get_cthis(athis,anext,pars)
    # c[s] is valid? (>=0 required by definition)
    if cthis < 0.0
        return (-9.99E99)::Float64
    else
        return ( u_CRRA(cthis, γ = pars.γ) + pars.β * src.linear_function_interpolation(anext, anext_vec, vnext_vec) )::Float64
    end # if
end # bellman_val

# ------------------------- standard solver
"""
    solve_dp(m::LifeCycleDatPkg ; arange::Tuple{Float64,Float64} = (0.0,sum(m.F) + a.a1), gridnum::Int = 200 )

standard solver for a deterministic Bellman equation which uses CRRA utility (with parameter γ),
where households do not intend to leave bequests.
The bellman equation follows the following inter-temporal budget:
``
A_s k_{s+1} = B_s k_{s} - E_{s} c_{s} + F_{s}, s=1,\\dots,S
``
where A,B,E,F are given parameters (deterministic);
for E, it is defined as member `E` in `LifeCycleDatPkg`.

Note: we use `sum(m.F)+a.a1` as the guess of ``a_s`` upper bound because someone
cannot consume more than his/her life-time total income when meeting borrowing constraint (`arange[1]=0.0`).
"""
function solve_dp(m::LifeCycleDatPkg ; arange::Tuple{Float64,Float64} = (0.0,sum(m.F) + m.a1), gridnum::Int = 200 )
    # validation
    @assert(-Inf < arange[1] < arange[2] < Inf, "arange must be finite and the lower bound must be less than the upper bound")
    @assert(1 < gridnum < Inf, "the number of grid points must be positive and finite integer")
    # case: m.S == 1 (only one period, no need to compute)
    if m.S == 1
        tmppars = BellmanParSet(m,1)
        return ( a = Float64[m.a1,], c = Float64[ get_cthis(m.a1,0.0,tmppars) ,] )::NamedTuple
    end # if
    # if not, go on
    # malloc
    local agrid = Array(LinRange(arange[1],arange[2],gridnum)) # griding vector of a_{s}
    local optv = zeros(gridnum,m.S) # policy function(s) component matrix, for value functions' values
    local opta = zeros(gridnum,m.S) # policy function(s) component matrix, for optimal asset levels
    local optc = zeros(gridnum,m.S) # policy function(s) component matrix, for optimal consumption
    local apath = zeros(m.S) # solved path for asset, the 1st element is m.a1
    local cpath = zeros(m.S) # solved path for consumption
    # ---------- let`s go
    # period: s = m.S
    idxS = m.S
        for x in 1:gridnum
            # first, fill the last column (policy function of period S) with griding asset(s)
            opta[x, m.S] = agrid[x]
            # then, compute consumption for every possible a_{S}; NOTE: no intend to leave bequest
            tmppars = BellmanParSet(m, idxS)
            optc[x, m.S] = get_cthis(agrid[x], 0.0, tmppars)
            # finally, use u(c[S]) as V_{S}
            optv[x, m.S] = u_CRRA(optc[x,m.S], γ = tmppars.γ)
        end # for x
    # period: 1 <= s < m.S (NOTE: this block compatible with the case that m.S == 2)
    for idxS in (m.S-1):-1:1  # NOTE: on this layer, we compute policy functions for every period
        for x in 1:gridnum  # NOTE: on this layer, we compute policy function points for every possible a[s]
            # first, define the parameter set used in this loop of x
            tmppars = BellmanParSet(m, idxS)
            # then, define a one-parameter function for Bisection search later
            # NOTE: we have to define it here because it depends on `x` and `idxS`
            local afunc(anext::Real) = bellman_val(agrid[x],  anext,  agrid,optv[:,idxS+1],tmppars)
            # next, compute a temporary value function series for griding a[s]
            # NOTE: this temp series is exactly the values of ``[u(c) + βV_{s+1}(a_{s+1})]`` for every possible a_{s+1};
            #       we will find its maxmimum as the value of ``V_{s}`` later, i.e. the Bellman equation's definition
            local tmpv = Float64[]
            for j in 1:gridnum; push!( tmpv, afunc(agrid[j]) ); end # for j
            # now, let's find the maximum of `tmpv`
            local tmploc = findmax(tmpv)[2]  # locate the maximum
            if (tmploc == 1) | (tmploc == gridnum)
                isfinite(tmpv[tmploc]) ? nothing : throw(DomainError("infeasible DP problem found when evaluating Bellman equation in period $(idxS)"))
                opta[x,idxS] = agrid[tmploc]
            else
                # if not the 1st or the last element of griding asset series, we use Golden Section search to get a more precise a[s] in a possible range
                # NOTE: assuming we find the maximum at location j, then we search in range (agrid[j-1], agrid[j+1])
                # NOTE: the golden_section() function checks if gold_lb or gold_ub is infinite or NaN; if infinite or NaN, an error will be thrown by the function
                local gold_lb = agrid[tmploc-1]; local gold_ub = agrid[tmploc+1]
                local tmpres = src.golden_section(afunc, gold_lb, gold_ub, atol = 1E-8, maxiter = 5000)
                # check if the solution is feasible
                if isfinite(tmpres[1]) & isfinite(tmpres[2]) # check optimal a[s+1] and V[s]
                    opta[x,idxS] = tmpres[1] # optimal a[s+1] (please note where I save it in `opta`)
                else
                    throw(ErrorException("no solution, stopped at evaluating the policy function of period $(idxS) at asset level $(agrid[x]); the stopping optimal V is $(tmpres[2])"))
                end # if
            end # if tmploc

            # finally, we evalute the (interpolated/precise/accurate/continuous/real) ``V_{s}`` and ``c^*_{s}`` and save
            # NOTE: we need this step because we (might) have get a non-grid asset point (opta[x,idxS]) as our optimum
            optv[x,idxS] = afunc(opta[x,idxS])
            optc[x,idxS] = get_cthis(agrid[x], opta[x,idxS], tmppars)

        end # for x
    end # for idxS

    # now, let's extract optimal paths conditional on a[1] (as input parameter in `m`)
    apath[1] = m.a1
    for x in 1:(m.S-1)
        # query & interpolate a[s+1] with src.linear_function_interpolation
        # NOTE: dont forget the policy function is defined as such a mapping: agrid (a[s]) --> opta[:,x] (optimal a[s+1])
        #       this helps to better understand `linear_function_interpolation()`
        apath[x+1] = src.linear_function_interpolation(apath[x], agrid, opta[:,x])
        # using inter-temporal budget to get c[s+1]
        tmppars = BellmanParSet(m, x)
        cpath[x] = get_cthis(apath[x],apath[x+1],tmppars)
    end # for x
    # finally, get c[S] (because we have the boundary condition that a[S+1]=0, the computing of c[S] needs to be separated from the above loop)
    tmppars = BellmanParSet(m, m.S)
    cpath[m.S] = get_cthis(apath[m.S],0.0,tmppars)

    # finally, returns the results in a NamedTuple
    return ( a = apath, c = cpath )::NamedTuple
end # solve_dp

end # Household

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
@testset "Getter functions for types" begin
  n = 5
  g = complete_graph(n)
  ctqw = CTQW(g)
  qws = QWSearch(ctqw, [1])

  @test QuantumWalk.graph(ctqw) == ctqw.graph
  @test QuantumWalk.graph(qws) == QuantumWalk.graph(qws.model)
  @test QuantumWalk.model(qws) == qws.model
  @test QuantumWalk.parameters(qws) == qws.parameters
end

ext: jl


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

@doc raw"""RuntimeClassSpec is a specification of a RuntimeClass. It contains parameters that are required to describe the RuntimeClass to the Container Runtime Interface (CRI) implementation, as well as any other components that need to understand how the pod will be run. The RuntimeClassSpec is immutable.

    IoK8sApiNodeV1alpha1RuntimeClassSpec(;
        overhead=nothing,
        runtimeHandler=nothing,
        scheduling=nothing,
    )

    - overhead::IoK8sApiNodeV1alpha1Overhead : Overhead represents the resource overhead associated with running a pod for a given RuntimeClass. For more details, see https://git.k8s.io/enhancements/keps/sig-node/20190226-pod-overhead.md This field is alpha-level as of Kubernetes v1.15, and is only honored by servers that enable the PodOverhead feature.
    - runtimeHandler::String : RuntimeHandler specifies the underlying runtime and configuration that the CRI implementation will use to handle pods of this class. The possible values are specific to the node & CRI configuration.  It is assumed that all handlers are available on every node, and handlers of the same name are equivalent on every node. For example, a handler called \"runc\" might specify that the runc OCI runtime (using native Linux containers) will be used to run the containers in a pod. The RuntimeHandler must conform to the DNS Label (RFC 1123) requirements and is immutable.
    - scheduling::IoK8sApiNodeV1alpha1Scheduling : Scheduling holds the scheduling constraints to ensure that pods running with this RuntimeClass are scheduled to nodes that support it. If scheduling is nil, this RuntimeClass is assumed to be supported by all nodes.
"""
mutable struct IoK8sApiNodeV1alpha1RuntimeClassSpec <: SwaggerModel
    overhead::Any # spec type: Union{ Nothing, IoK8sApiNodeV1alpha1Overhead } # spec name: overhead
    runtimeHandler::Any # spec type: Union{ Nothing, String } # spec name: runtimeHandler
    scheduling::Any # spec type: Union{ Nothing, IoK8sApiNodeV1alpha1Scheduling } # spec name: scheduling

    function IoK8sApiNodeV1alpha1RuntimeClassSpec(;overhead=nothing, runtimeHandler=nothing, scheduling=nothing)
        o = new()
        validate_property(IoK8sApiNodeV1alpha1RuntimeClassSpec, Symbol("overhead"), overhead)
        setfield!(o, Symbol("overhead"), overhead)
        validate_property(IoK8sApiNodeV1alpha1RuntimeClassSpec, Symbol("runtimeHandler"), runtimeHandler)
        setfield!(o, Symbol("runtimeHandler"), runtimeHandler)
        validate_property(IoK8sApiNodeV1alpha1RuntimeClassSpec, Symbol("scheduling"), scheduling)
        setfield!(o, Symbol("scheduling"), scheduling)
        o
    end
end # type IoK8sApiNodeV1alpha1RuntimeClassSpec

const _property_map_IoK8sApiNodeV1alpha1RuntimeClassSpec = Dict{Symbol,Symbol}(Symbol("overhead")=>Symbol("overhead"), Symbol("runtimeHandler")=>Symbol("runtimeHandler"), Symbol("scheduling")=>Symbol("scheduling"))
const _property_types_IoK8sApiNodeV1alpha1RuntimeClassSpec = Dict{Symbol,String}(Symbol("overhead")=>"IoK8sApiNodeV1alpha1Overhead", Symbol("runtimeHandler")=>"String", Symbol("scheduling")=>"IoK8sApiNodeV1alpha1Scheduling")
Base.propertynames(::Type{ IoK8sApiNodeV1alpha1RuntimeClassSpec }) = collect(keys(_property_map_IoK8sApiNodeV1alpha1RuntimeClassSpec))
Swagger.property_type(::Type{ IoK8sApiNodeV1alpha1RuntimeClassSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNodeV1alpha1RuntimeClassSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiNodeV1alpha1RuntimeClassSpec }, property_name::Symbol) =  _property_map_IoK8sApiNodeV1alpha1RuntimeClassSpec[property_name]

function check_required(o::IoK8sApiNodeV1alpha1RuntimeClassSpec)
    (getproperty(o, Symbol("runtimeHandler")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiNodeV1alpha1RuntimeClassSpec }, name::Symbol, val)
end


text:
using Gadfly, RDatasets, Compose

set_default_plot_size(6inch, 9inch)

vstack(
    plot(sin, 0, 25, xintercept=[0, pi, 2pi, 3pi], yintercept=[0, -1, 1],
        Geom.hline(style=[:dot,[1mm,1mm],:solid]), Geom.vline(style=:dashdot)),
    plot(dataset("datasets", "iris"), x="SepalLength", y="SepalWidth", Geom.point,
        yintercept=[2.5, 4.0], Geom.hline(color=["orange","red"], size=[2mm,3mm])),
    plot(dataset("ggplot2", "mpg"), x="Cty", y="Hwy", Geom.point,
        intercept=[0], slope=[1], Geom.abline(color="red", style=:dash),
        Guide.annotation(compose(context(), text(0.1w, 0.9h, "y=x", hright, vbottom),
            fill(colorant"red"), svgclass("marker")))) )


text:
include("joins.jl")
include("data.jl")
include("measurement.jl")
include("classification.jl")
include("misc.jl")


text:
module ExamplePackage

using ForwardDiff
include("extra_file.jl")

export my_f, derivative_of_my_f

end # module


text:
using .Delegates

a = Delegator(nothing)
b = Delegator("string")

d = Delegate()
c = Delegator(d)

@show Delegates.operation(a)
@show Delegates.operation(b)
@show Delegates.operation(c)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
SLS2 = RheoModelClass(
        # Model name
        name="SLS2",
        # Model parameters,
        p = (:G₀, :G₁, :η₁, :G₂, :η₂),
        # Relaxation modulus
        G = quote
              G₀ + G₁*exp(-t*G₁/η₁) + G₂*exp(-t*G₂/η₂)
            end,
        # Creep modulus
        J = quote
              tau1 = G₁/η₁
              tau2 = G₂/η₂
              Jbar(s) = (1/s^2)*(G₀/s+G₁*1/(s+tau1)+G₂*1/(s+tau2))
              invLaplace(s -> Jbar(s), t)
            end,
        # Network
        info= "
                                       ___
                 ___╱╲  ╱╲  ╱╲  ________| |_____
                 |    ╲╱  ╲╱  ╲╱       _|_|     |
                 |      G₂              η₂      |
                 |                     ___      |
              ___|__╱╲  ╱╲  ╱╲  ________| |_____|___
                 |    ╲╱  ╲╱  ╲╱       _|_|     |
                 |      G₁              η₁      |
                 |                              |
                 |__________╱╲  ╱╲  ╱╲  ________|
                              ╲╱  ╲╱  ╲╱
                                G₀
               "
        )

#  Experimental / work in progress
#  not functional yet
#
#  Generate array of parameter symbols and expression of the relaxation modulus
#  for the generalised maxwell model
#

function G_GSLS(n)

    params=(:k0,)
    G_terms=(:(k0),)

    for i=1:n
        k=Symbol("k_"*string(i))
        η=Symbol("η_"*string(i))
        e=:($k * exp.( (- $k/$η) * t))

        params=merge(params, (k,η))
        G_terms=merge(G_terms,(e,))
    end
    G_Expr=Expr(:call,:+,G_terms...)

    return (params, G_Expr)
end


text:
using Mimi
using MimiFUND
using Plots
using Statistics
using Distributions

# get the model
m = MimiFUND.get_model()
run(m)
explore(m)

sd = @defsim begin
    # Define random variables
    climatedynamics.climatesensitivity = truncated(Normal(2, 1), 0.1, 5) # truncated to be stricly positive
    save(climatedynamics.temp)
end

# Run 100 trials
si = run(sd, m, 100)

test = getdataframe(si, :climatedynamics, :temp)

test = filter("time" => ==(2100), test)

histogram(test[:,:temp], label = nothing, xlabel = "Temperature Increase", ylabel = "Nber of trials")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
using StackOverflow

"""
    makeanswersarray(json_obj::Dict)

Returns a Vector of `Answer` type extracted from a json object from stackexchange API.
"""
function makeanswersarray(json_obj::Dict)
    answerholder = Vector()
    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                creation_date = 0
                owner = Dict()
                last_activity_date = 0
                score = 0
                question_id = 0
                answer_id = 0
                locked_date = 0
                community_owned_date = 0
                is_accepted = 0
                last_edit_date = 0

                for (key, value) in item
                    if occursin("locked_date", key)
                        locked_date = value
                    elseif occursin("community_owned_date", key)
                        community_owned_date = value
                    elseif occursin("creation_date", key)
                        creation_date = value
                    elseif occursin("last_edit_date", key)
                        last_edit_date = value
                    elseif occursin("owner", key)
                        owner = value
                    elseif occursin("last_activity_date", key)
                        last_activity_date = value
                    elseif occursin("score", key)
                        score = value
                    elseif occursin("answer_id", key)
                        answer_id = value
                    elseif occursin("question_id", key)
                        question_id = value
                    end
                end
                answer = Answer(owner,
                 is_accepted,
                 score,
                 last_activity_date,
                 last_edit_date,
                 community_owned_date,
                 locked_date,
                 creation_date,
                 answer_id,
                 question_id
                 )
                push!(answerholder, answer)

            end
        end
    end
    answerholder
end

"""
    getquestions(sort::string = "creation", order::String = "desc")

Returns 30 most recent answers given a certain sort parameter.
"""
function getanswers(sort::String = "creation", order::String = "desc")
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/answers?order=$(order)&sort=$(sort)&site=stackoverflow")

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)
    answers = makeanswersarray(json_obj)

end

"""
    getanswerforquestion(q::Question; sort::String = "activity", order::String = "desc")

Returns a `Vector` of `Answer`s belonging to a `Question` `q`.
"""
function getanswerfromquestion(q::Question; sort::String = "activity", order::String = "desc")
    q_id = q.question_id
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions/$(q_id)/answers?order=$(order)&sort=$(sort)&site=stackoverflow")
    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)
    answers = makeanswersarray(json_obj)
end


text:
using Random
using LinearAlgebra
using PyPlot
using Seaborn
#PyPlot.svg(true)

function iterate_CG(xk, rk, pk, A)
    """
    Basic iteration of the conjugate gradient algorithm
    Parameters:
    xk: current iterate
    rk: current residual
    pk: current direction
    A: matrix of interest
    """

    # form products
    rkk = rk' * rk
    Apk = A * pk

    # construct step size
    αk = rkk / (pk' * Apk)

    # take a step in current conjugate direction
    xk_new = xk + αk * pk

    # construct new residual
    rk_new = rk + αk * Apk

    # construct new linear combination
    betak_new = (rk_new' * rk_new) / rkk

    # generate new conjugate vector
    pk_new = -rk_new + betak_new * pk

    return xk_new, rk_new, pk_new

end

function run_conjugate_gradient(x0, A, b; max_iter = 2000)
    """
    Conjugate gradient algorithm

    Parameters:
    x0: initial point
    A: matrix of interest
    b: vector in linear system (Ax = b)
    max_iter: max number of iterations to run CG
    """

    # initial iteration
    xk = x0
    rk = A * xk - b
    pk = -rk

    err = sum(abs.(rk))
    errors = [err]
    D = pk

    #while sum(abs.(rk)) > 10e-6
    for i in 1:max_iter
        xk, rk, pk = iterate_CG(xk, rk, pk, A)

        D = hcat(D,pk)

        err = sum(abs.(rk))
        push!(errors, err)

        # print iteration
        i % 500 == 0 && println("Iteration $i: error: $err")

        # break if we reach desired error
        if err < 10e-6
            println("Terminated in $i iterations")
            break
        end
    end

    return xk, errors, D
end

Random.seed!(123)
b = rand(5)

A = zeros(5,5)
A[1,1]=10; A[2,2]=10.1; A[3,3]=10.2; A[4,4]=2; A[5,5] = 1

x = A \ b

x0 = [0.3,0.1,0.1,0.1,0.1]
xsol, errors, D = run_conjugate_gradient(x0, A, b);
xsol

Seaborn.semilogy(0:length(errors)-1, errors)
xlabel("iteration"); ylabel("total absolute error")

# confirm these are 0
@assert D[:,1]' * A * D[:,2] < 10e-6
@assert D[:,1]' * A * D[:,3] < 10e-6
@assert D[:,1]' * A * D[:,4] < 10e-6

# confirm these aren't 0
@assert D[:,1]' * A * D[:,1] < 10e-6

# Construct a random symmetric positive definite matrix
Random.seed!(111)
A = rand(5,5); A = 0.5*(A+A'); A = A + 5*I;

eigvals(A)

x2 = A \ b

x0 = [0.,0.,0.,0.,0.]
xsol, errors, D = run_conjugate_gradient(x0, A, b);
xsol

Seaborn.semilogy(0:length(errors)-1, errors)
xlabel("iteration"); ylabel("total absolute error")

N = 1000
Random.seed!(123)
b = rand(N)
A = rand(N,N); A = 0.5*(A+A'); A = A + 50*I;
println("The condition number is: ", cond(A))

# distribution of the eigenvalues
plt[:hist](eigvals(A), bins=500); xlim(0,100)

x0 = zeros(N)
xsol, errors, D = run_conjugate_gradient(x0, A, b);
xsol;

Seaborn.semilogy(0:length(errors)-1, errors)
xlabel("iteration"); ylabel("total absolute error")

D[:,1]' * A * D[:,2]

D[:,1]' * A * D[:,10]

D[:,1]' * A * D[:,33]

Random.seed!(123)

# Generate random s.p.d. matrix
A = rand(N,N); A = 0.5*(A+A'); A = A + 13*I;
println("The condition number is: ", cond(A))

# distribution of the eigenvalues
plt[:hist](eigvals(A), bins=500); xlim(0,100)

x0 = zeros(N)
xsol, errors, D = run_conjugate_gradient(x0, A, b; max_iter=5000);

Seaborn.semilogy(0:length(errors)-1, errors)
xlabel("iteration"); ylabel("total absolute error")

D[:,1]' * A * D[:,2]

D[:,1]' * A * D[:,7]

D[:,1]' * A * D[:,75]


text:
"""
    exf(ex)

Return the name of a called function in a `:call` `Expr`.
"""
exf(ex) = ex.args[1]

"""
    exfargs(ex)

Return the arguments to a function call in a `:call` `Expr`.
"""
exfargs(ex) = ex.args[2:end]

name(x::Symbol) = string(x)
# name(x::DataType) = name(x.name)
# name(x::TypeName) = string(x.name)
name(x::DataType) = name(x.name, x.parameters)
name(x::TypeName, parameters) =
    length(parameters) > 0 ? string(x.name, "{:$(parameters[1])}") : string(x.name)

function as end

"""
    hasnulls(itr)

Return `true` if any element in `itr` is null; otherwise return `false`.
"""
@inline function hasnulls(itr::Any)::Bool
    res = false
    for el in itr
        if isa(el, Nullable)
            res |= isnull(el)
        end
    end
    return res
end


text:
using CairoMakie
using DifferentialEquations

noto_sans = "./resources/NotoSans-Regular.ttf"

tickfontsize    = 28
labelfontsize   = 32

function vanderpol!(du,u,μ,t)
    x,y = u
    du[1] = μ * (x - x^3/3 - y)
    du[2] = x/μ
end

function vanderpol(u,μ,t)
    du = zeros(2)
    vanderpol!(du,u,μ,t)
    return du
end

#%%
f = Figure(font=noto_sans, resolution=(800,400))

ax = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Time", ylabel = "State",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize)

curve_count = 100
u0 = 2*randn(2,curve_count)

solution_curves = zeros(401,2,curve_count)
timings = 0:0.05:20

for i in 1:curve_count
    μ = 2.0
    tspan = (0.0,20.0)

    prob = ODEProblem(vanderpol!,u0[:,i],tspan,μ)

    sol = solve(prob,saveat=0.05)
    solution_curves[:,1,i] = sol[1,:]
    solution_curves[:,2,i] = sol[2,:]

end

# Find the last peak
peaks = zeros(Integer,curve_count)
for i in 1:curve_count
    ~, peaks[i] = findmax(solution_curves[250:400,1,i])
end

# Now actually plot finally
for i in 1:curve_count
    lines!(ax,timings[1:100],solution_curves[peaks[i]:peaks[i]+100-1,1,i], color=:black)
    lines!(ax,timings[1:100],solution_curves[peaks[i]:peaks[i]+100-1,2,i], color=:blue)
end
xlims!(ax, (0,5.0))
ylims!(ax, (-4,4))
f

save("converge.pdf",f)


text:
module ImagesTests

using FactCheck, Base.Test, Images, Colors, FixedPointNumbers
using Graphics

testing_units = Int == Int64
if testing_units
    using SIUnits, SIUnits.ShortUnits
end

include("core.jl")
include("map.jl")
include("overlays.jl")
include("algorithms.jl")
include("edge.jl")
include("writemime.jl")
include("corner.jl")
include("distances.jl")

isinteractive() || FactCheck.exitstatus()

end

include("parallel.jl")


text:
n = 5 # rand(2:10)

A = matrixdepot("wilkinson", n)

@test A == matrixdepot("wilkinson", Float64, n)

@test issymmetric(Matrix(A))

@test diag(Matrix(A),1) ≈ ones(Float64, n - 1)

# symmetric along antidiagonal
@test issymmetric(Matrix(A)[n:-1:1,:])

θ = matrixdepot("wilkinson", 1)

println("'wilkinson' passed test...")


text:
using GXBeam
using LinearAlgebra
using DifferentialEquations
using Test
import Elliptic
using ForwardDiff

@testset "Linear Analysis of a Cantilever Partially Under a Uniform Distributed Load" begin

    nelem = 12

    # create points
    a = 0.3
    b = 0.7
    L = 1.0
    n1 = n3 = div(nelem, 3)
    n2 = nelem - n1 - n3
    x1 = range(0, a, length=n1+1)
    x2 = range(a, b, length=n2+1)
    x3 = range(b, L, length=n3+1)
    x = vcat(x1, x2[2:end], x3[2:end])
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints for each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # create compliance matrix for each beam element
    EI = 1e9
    stiffness = fill(Diagonal([0, 0, 0, 0, EI, 0]), nelem)

    # create the assembly
    assembly = Assembly(points, start, stop, stiffness=stiffness)

    # set prescribed conditions (fixed right endpoint)
    prescribed_conditions = Dict(
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0)
    )

    # create distributed load
    q = 1000
    distributed_loads = Dict()
    for ielem in n1+1:n1+n2
        distributed_loads[ielem] = DistributedLoads(assembly, ielem; fz = (s) -> q)
    end

    system, converged = static_analysis(assembly, prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads, linear=true)

    state = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)

    # analytical solution obtained using superposition
    initial_slope = -q/(6*EI)*((L-a)^3 - (L-b)^3)
    initial_deflection = q/(24*EI)*((L-a)^3*(3*L + a) - (L-b)^3*(3*L + b))
    analytical_M = function(x)
        if 0 < x <= a
            M = 0.0
        elseif a < x <= b
            M = q/2*(x-a)^2
        else
            M = q/2*((x-a)^2 - (x-b)^2)
        end
        return M
    end
    analytical_slope = function(x)
        slope = initial_slope
        if 0 < x <= a
            slope += 0.0
        elseif a < x <= b
            slope += q/(6*EI)*(x-a)^3
        else
            slope += q/(6*EI)*((x-a)^3 - (x-b)^3)
        end
        return slope
    end
    analytical_deflection = function(x)
        deflection = initial_deflection + initial_slope*x
        if 0 < x <= a
            deflection += 0.0
        elseif a < x <= b
            deflection += q/(24*EI)*(x-a)^4
        else
            deflection += q/(24*EI)*((x-a)^4 - (x-b)^4)
        end
        return deflection
    end

    # test element properties
    for i = 1:length(assembly.elements)
        xi = assembly.elements[i].x[1]
        @test isapprox(state.elements[i].u[3], analytical_deflection(xi), atol=1e-9)
        @test isapprox(state.elements[i].theta[2], -4*analytical_slope(xi)/4, atol=1e-9)
        @test isapprox(state.elements[i].M[2], -analytical_M(xi), atol=2)
    end

    # test point properties
    for i = 1:length(assembly.points)
        xi = assembly.points[i][1]
        @test isapprox(state.points[i].u[3], analytical_deflection(xi), atol=1e-8)
        @test isapprox(state.points[i].theta[2], -4*analytical_slope(xi)/4, atol=1e-7)
    end
end

@testset "Linear Analysis of a Beam Under a Linear Distributed Load" begin

    nelem = 16

    # create points
    L = 1
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints for each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # create compliance matrix for each beam element
    EI = 1e7
    compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

    # create assembly
    assembly = Assembly(points, start, stop, compliance=compliance)

    # set prescribed conditions
    prescribed_conditions = Dict(
        # simply supported left endpoint
        1 => PrescribedConditions(uz=0),
        # clamped right endpoint
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0)
    )

    # create distributed load
    qmax = 1000
    distributed_loads = Dict()
    for i = 1:nelem
        distributed_loads[i] = DistributedLoads(assembly, i; s1=x[i],
            s2=x[i+1], fz = (s) -> qmax*s)
    end

    # solve system
    system, converged = static_analysis(assembly, prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads, linear=true)

    # post-process the results
    state = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)

    # construct analytical solution
    analytical_deflection = (x) -> qmax*(1-x)^2/(120*EI)*(4 - 8*(1-x) + 5*(1-x)^2 - (1-x)^3)
    analytical_slope = (x) -> -qmax*(1-x)/(120*EI)*(8 - 24*(1-x) + 20*(1-x)^2 - 5*(1-x)^3)
    analytical_M = (x) -> qmax/120*(8 - 48*(1-x) + 60*(1-x)^2 - 20*(1-x)^3)

    # test element properties
    for i = 1:length(assembly.elements)
        xi = assembly.elements[i].x[1]
        @test isapprox(state.elements[i].u[3], analytical_deflection(xi), atol=1e-8)
        @test isapprox(state.elements[i].theta[2], -4*analytical_slope(xi)/4, atol=1e-7)
        @test isapprox(state.elements[i].M[2], -analytical_M(xi), atol=1)
    end

    # test point properties
    for i = 1:length(assembly.points)
        xi = assembly.points[i][1]
        @test isapprox(state.points[i].u[3], analytical_deflection(xi), atol=1e-8)
        @test isapprox(state.points[i].theta[2], -4*analytical_slope(xi)/4, atol=1e-8)
    end
end

@testset "Nonlinear Analysis of a Cantilever Subjected to a Constant Tip Load" begin

    L = 1
    EI = 1e6

    # shear force (applied at end)
    λ = 0:0.5:16
    p = EI/L^2
    P = λ*p

    # create points
    nelem = 16
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # compliance matrix for each beam element
    compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop, compliance=compliance)

    # pre-initialize system storage
    system = System(assembly, true)

    # run an analysis for each prescribed tip load
    states = Vector{AssemblyState{Float64}}(undef, length(P))
    for i = 1:length(P)

        # create dictionary of prescribed conditions
        prescribed_conditions = Dict(
            # fixed left side
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # shear force on right tip
            nelem+1 => PrescribedConditions(Fz = P[i])
        )

        # perform a static analysis
        static_analysis!(system, assembly, prescribed_conditions=prescribed_conditions)

        # post-process the results
        states[i] = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)

    end

    # construct analytical solution
    δ = range(pi/4, pi/2, length=10^5)[2:end-1]

    k = @. cos(pi/4)/sin(δ)
    λ_a = @. (Elliptic.F(pi/2, k^2) - Elliptic.F(δ,  k^2))^2

    θ_a = @. 2*(pi/4 - acos(k))

    ξ_a = @. sqrt(2*sin(θ_a)/λ_a) .- 1

    η_a = @. 1-2/sqrt(λ_a)*(Elliptic.E(pi/2, k^2) - Elliptic.E(δ, k^2))

    # test tip displacements
    for i = 1:length(P)
        i_a = argmin(abs.(λ[i] .- λ_a))
        @test isapprox(states[i].points[end].u[1]/L, ξ_a[i_a], atol=1e-3)
        @test isapprox(states[i].points[end].u[3]/L, η_a[i_a], atol=1e-3)
        @test isapprox(states[i].points[end].theta[2], -4*tan(θ_a[i_a]/4), atol=1e-2)
    end
end

@testset "Nonlinear Analysis of a Cantilever Subjected to a Constant Moment" begin

    L = 12 # inches
    h = w = 1 # inches
    E = 30e6 # lb/in^4 Young's Modulus

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

    # bending moment (applied at end)
    # note that solutions for λ > 1.8 do not converge
    λ = [0.0, 0.4, 0.8, 1.2, 1.6, 1.8, 2.0]
    m = pi*E*Iyy/L
    M = λ*m

    # create points
    nelem = 16
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints for each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # compliance matrix for each beam element
    compliance = fill(Diagonal([1/(E*A), 0, 0, 0, 1/(E*Iyy), 1/(E*Izz)]), nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop, compliance=compliance)

    # pre-initialize system storage
    system = System(assembly, true)

    # run an analysis for each prescribed bending moment
    states = Vector{AssemblyState{Float64}}(undef, length(M))
    for i = 1:length(M)

        # create dictionary of prescribed conditions
        prescribed_conditions = Dict(
            # fixed left side
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # moment on right side
            nelem+1 => PrescribedConditions(Mz = M[i])
        )

        # perform a static analysis
        static_analysis!(system, assembly, prescribed_conditions=prescribed_conditions)

        # post-process the results
        states[i] = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)

    end

    # analytical solution (ρ = E*I/M)
    analytical(x, ρ) = ifelse(ρ == Inf, zeros(3), [ρ*sin(x/ρ)-x, ρ*(1-cos(x/ρ)), 0])

    # test element properties
    for i = 1:length(M)
        for ielem = 1:length(assembly.elements)
            xi = assembly.elements[ielem].x[1]
            u_a, v_a, w_a = analytical(xi, E*Iyy/M[i])
            @test isapprox(states[i].elements[ielem].u[1], u_a, atol=5e-2)
            @test isapprox(states[i].elements[ielem].u[2], v_a, atol=5e-2)
        end

        # test point properties
        for ipoint = 1:length(assembly.points)
            xi = assembly.points[ipoint][1]
            u_a, v_a, w_a = analytical(xi, E*Iyy/M[i])
            @test isapprox(states[i].points[ipoint].u[1], u_a, atol=5e-2)
            @test isapprox(states[i].points[ipoint].u[2], v_a, atol=5e-2)
        end
    end
end

@testset "Nonlinear Analysis of the Bending of a Curved Beam in 3D Space" begin

    # problem constants
    R = 100
    L = R*pi/4 # inches
    h = w = 1 # inches
    E = 1e7 # psi Young's Modulus
    ν = 0.0
    G = E/(2*(1+ν))

    # beam starting point, frame, and curvature
    r = [0, 0, 0]
    frame = [0 -1 0; 1 0 0; 0 0 1]
    curvature = [0, 0, -1/R]

    # cross section properties
    A = h*w
    Ay = A
    Az = A
    Iyy = w*h^3/12
    Izz = w^3*h/12
    J = Iyy + Izz

    # discretize the beam
    nelem = 16
    ΔL, xp, xm, Cab = discretize_beam(L, r, nelem; frame=frame, curvature = curvature)

    # force
    P = 600 # lbs

    # index of left and right endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # compliance matrix for each beam element
    compliance = fill(Diagonal([1/(E*A), 1/(G*Ay), 1/(G*Az), 1/(G*J), 1/(E*Iyy), 1/(E*Izz)]), nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(xp, start, stop, compliance=compliance, frames=Cab,
        lengths=ΔL, midpoints=xm)

    # create dictionary of prescribed conditions
    prescribed_conditions = Dict(
        # fixed left endpoint
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        # force on right endpoint
        nelem+1 => PrescribedConditions(Fz=P)
    )

    # perform static analysis
    system, converged = static_analysis(assembly, prescribed_conditions=prescribed_conditions)

    # post-process results
    state = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)

    # Results from "Large Displacement Analysis of Three-Dimensional Beam
    # Structures" by Bathe and Bolourch:
    # - Tip Displacement: [-13.4, -23.5, 53.4]

    # Note that these results are comparing computational solutions, rather than
    # the computational to the analytical solution, so some variation is expected.

    @test isapprox(state.points[end].u[1], -13.4, atol=0.2) # -13.577383726758564
    @test isapprox(state.points[end].u[2], -23.5, atol=0.1) # -23.545303336988038
    @test isapprox(state.points[end].u[3],  53.4, atol=0.1) #  53.45800757548929
end

@testset "Rotating Beam with a Swept Tip" begin
    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

    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 #hide
        local xp_b1, xm_b1, Cab_b1 #hide
        local cs, ss #hide
        local L_b2, r_b2, nelem_b2, frame_b2, lengths_b2 #hide
        local xp_b2, xm_b2, Cab_b2 #hide
        local nelem, points, start, stop #hide
        local lengths, midpoints, Cab, compliance, mass, assembly #hide

        # 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

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

    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]
    ]

    for k = 1:length(experiment_frequencies)
        for j = 1:length(experiment_sweep)
            for i = 1:length(experiment_rpm)
                ii = argmin(abs.(rpm .- experiment_rpm[i]))
                jj = argmin(abs.(sweep*180/pi .- experiment_sweep[j]))
                kk = indices[k]
                @test isapprox(frequency[kk][jj,ii], experiment_frequencies[k][i,j], atol=1, rtol=0.1)
            end
        end
    end

    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
    ]

    for k = 1:size(experiment_frequencies, 1)
        for j = 1:length(experiment_sweep)
            ii = argmin(abs.(rpm .- 750))
            jj = argmin(abs.(sweep*180/pi .- experiment_sweep[j]))
            kk = indices[k]
            @test isapprox(frequency[kk][jj,ii], experiment_frequencies[k,j], rtol=0.1)
        end
    end
end

@testset "Nonlinear Dynamic Analysis of a Wind Turbine Blade" begin

    L = 60 # m

    # create points
    nelem = 10
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # stiffness matrix for each beam element
    stiffness = fill(
        [2.389e9  1.524e6  6.734e6 -3.382e7 -2.627e7 -4.736e8
         1.524e6  4.334e8 -3.741e6 -2.935e5  1.527e7  3.835e5
         6.734e6 -3.741e6  2.743e7 -4.592e5 -6.869e5 -4.742e6
        -3.382e7 -2.935e5 -4.592e5  2.167e7 -6.279e5  1.430e6
        -2.627e7  1.527e7 -6.869e5 -6.279e5  1.970e7  1.209e7
        -4.736e8  3.835e5 -4.742e6  1.430e6  1.209e7  4.406e8],
        nelem)

    # mass matrix for each beam element
    mass = fill(
        [258.053      0.0        0.0      0.0      7.07839  -71.6871
           0.0      258.053      0.0     -7.07839  0.0        0.0
           0.0        0.0      258.053   71.6871   0.0        0.0
           0.0       -7.07839   71.6871  48.59     0.0        0.0
           7.07839    0.0        0.0      0.0      2.172      0.0
         -71.6871     0.0        0.0      0.0      0.0       46.418],
         nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop; stiffness=stiffness, mass=mass)

    # simulation time
    tvec = 0:0.001:2.0

    # prescribed conditions
    prescribed_conditions = (t) -> begin
        Dict(
            # fixed left side
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # force on right side
            nelem+1 => PrescribedConditions(Fz=1e5*sin(20*t))
            )
    end

    system, history, converged = time_domain_analysis(assembly, tvec; prescribed_conditions=prescribed_conditions)

    @test converged
end

@testset "Nonlinear Static Analysis of a Joined-Wing" begin

    # Set endpoints of each beam
    p1 = [-7.1726, -12, -3.21539]
    p2 = [-5.37945, -9, -2.41154]
    p3 = [-3.5863, -6, -1.6077]
    p4 = [-1.79315, -3, -0.803848]
    p5 = [0, 0, 0]
    p6 = [7.1726, -12, 3.21539]

    # get transformation matrix for left beams

    # transformation from intermediate to global frame
    tmp1 = sqrt(p1[1]^2 + p1[2]^2)
    c1, s1 = -p1[1]/tmp1, -p1[2]/tmp1
    rot1 = [c1 -s1 0; s1 c1 0; 0 0 1]

    # transformation from local to intermediate frame
    tmp2 = sqrt(p1[1]^2 + p1[2]^2 + p1[3]^2)
    c2, s2 = tmp1/tmp2, -p1[3]/tmp2
    rot2 = [c2 0 -s2; 0 1 0; s2 0 c2]

    Cab_1 = rot1*rot2

    # get transformation matrix for right beam

    # transformation from intermediate frame to global frame
    tmp1 = sqrt(p6[1]^2 + p6[2]^2)
    c1, s1 = p6[1]/tmp1, p6[2]/tmp1
    rot1 = [c1 -s1 0; s1 c1 0; 0 0 1]

    # transformation from local beam frame to intermediate frame
    tmp2 = sqrt(p6[1]^2 + p6[2]^2 + p6[3]^2)
    c2, s2 = tmp1/tmp2, p6[3]/tmp2
    rot2 = [c2 0 -s2; 0 1 0; s2 0 c2]

    Cab_2 = rot1*rot2

    # beam 1
    L_b1 = norm(p2-p1)
    r_b1 = p1
    nelem_b1 = 5
    lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1, frame=Cab_1)
    compliance_b1 = fill(Diagonal([1.05204e-9, 3.19659e-9, 2.13106e-8, 1.15475e-7, 1.52885e-7, 7.1672e-9]), nelem_b1)

    # beam 2
    L_b2 = norm(p3-p2)
    r_b2 = p2
    nelem_b2 = 5
    lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2, frame=Cab_1)
    compliance_b2 = fill(Diagonal([1.24467e-9, 3.77682e-9, 2.51788e-8, 1.90461e-7, 2.55034e-7, 1.18646e-8]), nelem_b2)

    # beam 3
    L_b3 = norm(p4-p3)
    r_b3 = p3
    nelem_b3 = 5
    lengths_b3, xp_b3, xm_b3, Cab_b3 = discretize_beam(L_b3, r_b3, nelem_b3, frame=Cab_1)
    compliance_b3 = fill(Diagonal([1.60806e-9, 4.86724e-9, 3.24482e-8, 4.07637e-7, 5.57611e-7, 2.55684e-8]), nelem_b3)

    # beam 4
    L_b4 = norm(p5-p4)
    r_b4 = p4
    nelem_b4 = 5
    lengths_b4, xp_b4, xm_b4, Cab_b4 = discretize_beam(L_b4, r_b4, nelem_b4, frame=Cab_1)
    compliance_b4 = fill(Diagonal([2.56482e-9, 7.60456e-9, 5.67609e-8, 1.92171e-6, 2.8757e-6, 1.02718e-7]), nelem_b4)

    # beam 5
    L_b5 = norm(p6-p5)
    r_b5 = p5
    nelem_b5 = 20
    lengths_b5, xp_b5, xm_b5, Cab_b5 = discretize_beam(L_b5, r_b5, nelem_b5, frame=Cab_2)
    compliance_b5 = fill(Diagonal([2.77393e-9, 7.60456e-9, 1.52091e-7, 1.27757e-5, 2.7835e-5, 1.26026e-7]), nelem_b5)

    # combine elements and points into one array
    nelem = nelem_b1 + nelem_b2 + nelem_b3 + nelem_b4 + nelem_b5
    points = vcat(xp_b1, xp_b2[2:end], xp_b3[2:end], xp_b4[2:end], xp_b5[2:end])
    start = 1:nelem
    stop = 2:nelem + 1
    lengths = vcat(lengths_b1, lengths_b2, lengths_b3, lengths_b4, lengths_b5)
    midpoints = vcat(xm_b1, xm_b2, xm_b3, xm_b4, xm_b5)
    Cab = vcat(Cab_b1, Cab_b2, Cab_b3, Cab_b4, Cab_b5)
    compliance = vcat(compliance_b1, compliance_b2, compliance_b3, compliance_b4, compliance_b5)

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

    Fz = range(0, 70e3, length=141)

    # pre-allocate memory to reduce run-time
    system = System(assembly, true)

    linear_states = Vector{AssemblyState{Float64}}(undef, length(Fz))
    for i = 1:length(Fz)

        # create dictionary of prescribed conditions
        prescribed_conditions = Dict(
            # fixed endpoint on beam 1
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # force applied on point 4
            nelem_b1 + nelem_b2 + nelem_b3 + nelem_b4 + 1 => PrescribedConditions(Fz = Fz[i]),
            # fixed endpoint on last beam
            nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        )

        _, converged = static_analysis!(system, assembly, prescribed_conditions=prescribed_conditions, linear=true)

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

        @test converged
    end

    reset_state!(system)
    nonlinear_states = Vector{AssemblyState{Float64}}(undef, length(Fz))
    for i = 1:length(Fz)

        # create dictionary of prescribed conditions
        prescribed_conditions = Dict(
            # fixed endpoint on beam 1
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # force applied on point 4
            nelem_b1 + nelem_b2 + nelem_b3 + nelem_b4 + 1 => PrescribedConditions(Fz = Fz[i]),
            # fixed endpoint on last beam
            nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        )

        _, converged = static_analysis!(system, assembly, prescribed_conditions=prescribed_conditions,
            reset_state = false)

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

        @test converged
    end

    reset_state!(system)
    nonlinear_follower_states = Vector{AssemblyState{Float64}}(undef, length(Fz))
    for i = 1:length(Fz)
        # create dictionary of prescribed conditions
        prescribed_conditions = Dict(
            # fixed endpoint on beam 1
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # force applied on point 4
            nelem_b1 + nelem_b2 + nelem_b3 + nelem_b4 + 1 => PrescribedConditions(Fz_follower = Fz[i]),
            # fixed endpoint on last beam
            nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        )

        _, converged = static_analysis!(system, assembly, prescribed_conditions=prescribed_conditions,
            reset_state = false)

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

        @test converged
    end
end

@testset "Nonlinear Dynamic Analysis of a Joined-Wing" begin

    # Set endpoints of each beam
    p1 = [0, 0, 0]
    p2 = [-7.1726, -12, -3.21539]
    p3 = [7.1726, -12,  3.21539]

    Cab_1 = [
    0.5         0.866025  0.0
    0.836516    -0.482963  0.258819
    0.224144     -0.12941   -0.965926
    ]

    Cab_2 = [
    0.5         0.866025  0.0
    -0.836516    0.482963 0.258819
    0.224144    -0.12941   0.965926
    ]

    # beam 1
    L_b1 = norm(p1-p2)
    r_b1 = p2
    nelem_b1 = 8
    lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1, frame=Cab_1)

    # beam 2
    L_b2 = norm(p3-p1)
    r_b2 = p1
    nelem_b2 = 8
    lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2, frame=Cab_2)

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

    # assign all beams the same compliance and mass matrix
    compliance = fill(Diagonal([2.93944738387698e-10, 8.42991725049126e-10, 3.38313996669689e-08,
        4.69246721094557e-08, 6.79584100559513e-08, 1.37068861370898e-09]), nelem)
    mass = fill(Diagonal([4.86e-2, 4.86e-2, 4.86e-2,
        1.0632465e-2, 2.10195e-4, 1.042227e-2]), nelem)

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

    # time
    tvec = range(0, 0.04, length=1001)

    F_L = (t) -> begin
        if 0.0 <= t < 0.01
            1e6*t
        elseif 0.01 <= t < 0.02
            -1e6*(t-0.02)
        else
            zero(t)
        end
    end

    F_S = (t) -> begin
        if 0.0 <= t < 0.02
            5e3*(1-cos(pi*t/0.02))
        else
            1e4
        end
    end

    # assign boundary conditions and point load
    prescribed_conditions = (t) -> begin
        Dict(
        # fixed endpoint on beam 1
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        # force applied on point 4
        nelem_b1 + 1 => PrescribedConditions(Fx=F_L(t), Fy=F_L(t), Fz=F_S(t)),
        # fixed endpoint on last beam
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        )
    end

    system, history, converged = time_domain_analysis(assembly, tvec;
        prescribed_conditions=prescribed_conditions)

    @test converged
end

@testset "DifferentialEquations" begin

    L = 60 # m

    # create points
    nelem = 10
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # stiffness matrix for each beam element
    stiffness = fill(
        [2.389e9  1.524e6  6.734e6 -3.382e7 -2.627e7 -4.736e8
         1.524e6  4.334e8 -3.741e6 -2.935e5  1.527e7  3.835e5
         6.734e6 -3.741e6  2.743e7 -4.592e5 -6.869e5 -4.742e6
        -3.382e7 -2.935e5 -4.592e5  2.167e7 -6.279e5  1.430e6
        -2.627e7  1.527e7 -6.869e5 -6.279e5  1.970e7  1.209e7
        -4.736e8  3.835e5 -4.742e6  1.430e6  1.209e7  4.406e8],
        nelem)

    # mass matrix for each beam element
    mass = fill(
        [258.053      0.0        0.0      0.0      7.07839  -71.6871
           0.0      258.053      0.0     -7.07839  0.0        0.0
           0.0        0.0      258.053   71.6871   0.0        0.0
           0.0       -7.07839   71.6871  48.59     0.0        0.0
           7.07839    0.0        0.0      0.0      2.172      0.0
         -71.6871     0.0        0.0      0.0      0.0       46.418],
         nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop; stiffness=stiffness, mass=mass)

    # prescribed conditions
    prescribed_conditions = (t) -> begin
        Dict(
            # fixed left side
            1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
            # force on right side
            nelem+1 => PrescribedConditions(Fz=1e5*sin(20*t))
            )
    end

    # define simulation time
    tspan = (0.0, 2.0)

    # run initial condition analysis to get consistent set of initial conditions
    system, converged = initial_condition_analysis(assembly, tspan[1]; prescribed_conditions)

    # construct ODEProblem
    prob = ODEProblem(system, assembly, tspan; prescribed_conditions)

    # solve ODEProblem
    sol = solve(prob, Rodas4())

    # test that solution worked
    @test sol.t[end] == 2.0

    # construct DAEProblem
    prob = DAEProblem(system, assembly, tspan; prescribed_conditions)

    # solve DAEProblem
    sol = solve(prob, DABDF2())

    # test that solution worked
    @test sol.t[end] == 2.0
end

@testset "ForwardDiff" begin

    # Linear Analysis of a Beam Under a Linear Distributed Load

    function linear_analysis_test_with_AD(length) # this should affect just about everything

        nelem = 16

        # create points
        L = length[1]
        x = collect(range(0, L, length=nelem+1))
        y = zero(x)
        z = zero(x)

        points = [[x[i],y[i],z[i]] for i = 1:size(x,1)]

        # index of endpoints for each beam element
        start = 1:nelem
        stop = 2:nelem+1

        # create compliance matrix for each beam element
        EI = 1e7
        compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

        # create assembly
        assembly = Assembly(points, start, stop, compliance=compliance)

        # set prescribed conditions
        prescribed_conditions = Dict(
            # simply supported left endpoint
            1 => PrescribedConditions(uz=0),
            # clamped right endpoint
            nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0)
        )

        # create distributed load
        qmax = 1000
        distributed_loads = Dict()
        for i = 1:nelem
            distributed_loads[i] = DistributedLoads(assembly, i; s1=x[i],
                s2=x[i+1], fz = (s) -> qmax*s)
        end

        # solve system
        system, converged = static_analysis(assembly, prescribed_conditions=prescribed_conditions,
            distributed_loads=distributed_loads, linear=true)

        return system.x
    end

    # run FrowardDiff - no specific test, just make sure it runs fine
    J = ForwardDiff.jacobian(linear_analysis_test_with_AD, [1.0]) #length=1
end

@testset "Jacobian and Mass Matrix Calculations" begin

    L = 60 # m

    # create points
    nelem = 1
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    # index of endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    # stiffness matrix for each beam element
    stiffness = fill(
        [2.389e9  1.524e6  6.734e6 -3.382e7 -2.627e7 -4.736e8
         1.524e6  4.334e8 -3.741e6 -2.935e5  1.527e7  3.835e5
         6.734e6 -3.741e6  2.743e7 -4.592e5 -6.869e5 -4.742e6
        -3.382e7 -2.935e5 -4.592e5  2.167e7 -6.279e5  1.430e6
        -2.627e7  1.527e7 -6.869e5 -6.279e5  1.970e7  1.209e7
        -4.736e8  3.835e5 -4.742e6  1.430e6  1.209e7  4.406e8],
        nelem)

    # mass matrix for each beam element
    mass = fill(
        [258.053      0.0        0.0      0.0      7.07839  -71.6871
           0.0      258.053      0.0     -7.07839  0.0        0.0
           0.0        0.0      258.053   71.6871   0.0        0.0
           0.0       -7.07839   71.6871  48.59     0.0        0.0
           7.07839    0.0        0.0      0.0      2.172      0.0
         -71.6871     0.0        0.0      0.0      0.0       46.418],
         nelem)

    # create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop; stiffness=stiffness, mass=mass)

    # prescribed conditions
    pcond = Dict(
        # fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        )

    # distributed loads
    dload = Dict()

    # --- Static System --- #
    static_system = System(assembly, true)

    force_scaling = static_system.force_scaling
    mass_scaling = static_system.mass_scaling
    irow_point = static_system.irow_point
    irow_elem = static_system.irow_elem
    irow_elem1 = static_system.irow_elem1
    irow_elem2 = static_system.irow_elem2
    icol_point = static_system.icol_point
    icol_elem = static_system.icol_elem

    x = rand(length(static_system.x))
    J = similar(x, length(x), length(x))

    f = (x) -> GXBeam.system_residual!(similar(x), x, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem)

    GXBeam.system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    @test all(isapprox.(J, ForwardDiff.jacobian(f, x), atol=1e-10))

    # --- Dynamic System --- #
    system = System(assembly, false)

    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem
    x0 = rand(3)
    v0 = rand(3)
    ω0 = rand(3)

    x = rand(length(system.x))
    J = similar(x, length(x), length(x))

    f = (x) -> GXBeam.system_residual!(similar(x), x, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0)

    GXBeam.system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point,
        icol_elem, x0, v0, ω0)

    @test all(isapprox.(J, ForwardDiff.jacobian(f, x), atol=1e-10))

    # initial condition residual

    u0 = [rand(3) for ielem = 1:length(assembly.elements)]
    theta0 = [rand(3) for ielem = 1:length(assembly.elements)]
    udot0 = [rand(3) for ielem = 1:length(assembly.elements)]
    thetadot0 = [rand(3) for ielem = 1:length(assembly.elements)]

    x = rand(length(system.x))
    J = similar(x, length(x), length(x))

    f = (x) -> GXBeam.system_residual!(similar(x), x, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, u0, theta0, udot0, thetadot0)

    GXBeam.system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point,
        icol_elem, x0, v0, ω0, u0, theta0, udot0, thetadot0)

    @test all(isapprox.(J, ForwardDiff.jacobian(f, x), atol=1e-10))

    # time domain residual

    udot = [rand(3) for ielem = 1:length(assembly.elements)]
    θdot = [rand(3) for ielem = 1:length(assembly.elements)]
    Pdot = [rand(3) for ielem = 1:length(assembly.elements)]
    Hdot = [rand(3) for ielem = 1:length(assembly.elements)]
    dt = rand()

    x = rand(length(system.x))
    J = similar(x, length(x), length(x))

    f = (x) -> GXBeam.system_residual!(similar(x), x, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, udot, θdot, Pdot, Hdot, dt)

    GXBeam.system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point,
        icol_elem, x0, v0, ω0, udot, θdot, Pdot, Hdot, dt)

    @test all(isapprox.(J, ForwardDiff.jacobian(f, x), atol=1e-10))

    # system jacobian and mass matrix

    x = rand(length(system.x))
    dx = rand(length(x))
    J = similar(x, length(x), length(x))
    M = similar(x, length(x), length(x))

    fx = (x) -> GXBeam.dynamic_system_residual!(similar(x), x, dx, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0)

    fdx = (dx) -> GXBeam.dynamic_system_residual!(similar(dx), x, dx, assembly, pcond, dload,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0)

    GXBeam.dynamic_system_jacobian!(J, x, dx, assembly,
        pcond, dload, force_scaling, mass_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0)

    GXBeam.system_mass_matrix!(M, x, assembly, force_scaling, mass_scaling, irow_point,
        irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    @test all(isapprox.(J, ForwardDiff.jacobian(fx, x), atol=1e-10))

    @test all(isapprox.(M, ForwardDiff.jacobian(fdx, dx), atol=1e-10))

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -8.186501529651974 <= q <= 9.622483160840265)

                     add_NL_constraint(m, :(sigmoid(-0.9342689023247241 + -0.06643065035242302*sigmoid(-0.6967397920955185 + 0.4612305646899797*sigmoid(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.5917923558177942*sigmoid(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))) + 0.13022424784864484*sigmoid(-0.8092681495291991 + -0.6577545463201151*sigmoid(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.816846389391447*sigmoid(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))) + sigmoid(0.9502905692637253 + 0.2731932353707758*sigmoid(-0.6967397920955185 + 0.4612305646899797*sigmoid(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.5917923558177942*sigmoid(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))) + 0.7525923835832526*sigmoid(-0.8092681495291991 + -0.6577545463201151*sigmoid(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])) + 0.816846389391447*sigmoid(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#Assumes cluster started elsewhere
#register global variables so that host and ip are available elsewhere
function init(; ip::AbstractString = "localhost", port::Int = 54321)

	global H2Oip = ip
	global H2Oport = port
	return nothing

end

#Shut down a cluster
#Prompt set to false because prompt doesn't show up in Jupyter Notebook
#Using PyCall here to punt on re-writing the logic for time being
#shutdown(;conn = nothing, prompt::Bool = false) = h2o.shutdown(conn = conn, prompt = prompt)

#Return cluster info (after you've started the cluster)
function cluster_info()

	results = h2o_get("Cloud", "skip_ticks=true")
	result_json = JSON.parse(readall(results))

	return result_json

end


text:
@testset "Bridges" begin
    include("Constraint/Constraint.jl")
end


text:
# BAI
# ===
#
# Index for BAM files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# An index type for the BAM file format.
struct BAI
    # BGZF file index
    index::BGZFIndex

    # number of unmapped reads
    n_no_coors::Nullable{Int}
end

"""
    BAI(input::AbstractString)

Load a BAI index from `input`.
"""
function BAI(input::AbstractString)
    return open(x->BAI(x), input)
end

"""
    BAI(input::IO)

Load a BAI index from `input`.
"""
function BAI(input::IO)
    # check magic bytes
    B = read(input, UInt8)
    A = read(input, UInt8)
    I = read(input, UInt8)
    x = read(input, UInt8)
    if B != UInt8('B') || A != UInt8('A') || I != UInt8('I') || x != 0x01
        error("input is not a valid BAI file")
    end

    # read contents
    n_refs = read(input, Int32)
    index = read_bgzfindex(input, n_refs)
    if !eof(input)
        n_no_coors = Nullable{Int}(read(input, UInt64))
    else
        n_no_coors = Nullable{Int}()
    end

    return BAI(index, n_no_coors)
end

function overlapchunks(bai::BAI, seqid::Integer, interval::UnitRange{<:Integer})
    return overlapchunks(bai.index, seqid, first(interval):last(interval))
end


text:
using CuArrays, LinearAlgebra, Images, FileIO

MINFLOAT32 = nextfloat( Float32( -Inf ) );

l22d( x1, y1, x2, y2 ) = ( x1 - x2 )^2 + ( y1 - y2 )^2;

cost2dToMatrix( cost, n, m ) = [ cost( i÷n + 1, i%n + 1, j÷n + 1, j%n + 1) for i in 0:(n*m - 1), j in 0:(n*m - 1) ];

zrs( n... ) = CuArrays.fill( 0.0f0, n... )

function entropy( dist )
    dist /= sum( dist ) # we want it to sum to one
    return -(dist' * max.( log.( dist ), MINFLOAT32 ))
end

function wassBarycenters( dists, weights, precision, costMatrix, n_iter, sharpen=false )
	n, k = size( dists )
	nbarys = size( weights, 2 )
	# only accept probability distributions
	totalMass = sum( dists )
	@assert isapprox( totalMass, k ) && totalMass == sum( abs.( dists ) )

	logmus = max.( log.( dists ), MINFLOAT32 )
	# our barycenter should not have higher entropy than the inputs if we sharpen
	entropicBound = -minimum( diag( dists' * logmus ) )

	# renormalizing to avoid underflows in division
	costMatrix = costMatrix ./ max( costMatrix... )
	K = CuArray( exp.( -costMatrix/precision ) )
	weights = CuArray( weights )

	vs, ws, ds = zrs( n, k, nbarys ), zrs( n, k, nbarys ), zrs( n, k, nbarys )
	logbarycenters = zrs( n, nbarys )

	for j in 1:n_iter
		for i in 1:nbarys # CUDA.jl is smart and does this in parallel
			# ws[i] = dists[i] ./ K*vs[i]
			ws[:,:,i] = logmus - log.( K*exp.( vs[:,:,i] ) )
			# ds[i] = vs[i]  .* K*ws[i]
			ds[:,:,i] = vs[:,:,i] + log.( K*exp.( ws[:,:,i] ) )

			logbarycenters[:,i] = ds[:,:,i] * weights[:,i]

			if sharpen && entropy( exp.( logbarycenters[:,i] ) ) > entropicBound
				logbarycenters[:,i] *= 2
			end

			vs[:,:,i] = vs[:,:,i] - ds[:,:,i] .+ logbarycenters[:,i]
		end
	end
	return Array( exp.( logbarycenters ) )
end

function wassBarycenters2d( dists, weights, precision, cost, n_iter, sharpen=false )
	n, m = size( dists[1] )
	nbarys = size( weights, 2)
	C = cost2dToMatrix( cost, n, m )
	dists = CuArray( hcat( ([dist...] for dist in dists)... ) )
	bar = wassBarycenters( dists, weights, precision, C, n_iter, sharpen )
	return [ reshape( bar[:,i], (n, m) ) for i in 1:nbarys ]
end

function imageBarycenter( images, t, sharpen=false )
	inverted = [ 1 .- channelview( img ) for img in images ]
	dists    = [ img/sum( img ) for img in inverted ]
	bar = wassBarycenters2d( dists, [1-t, t], 0.0005, l22d, 100, sharpen )
	return Gray.(1 .- bar[1]./maximum( bar[1] ))
end

function imageBarycenters( images, step=0.05, sharpen=false )
	weights = [ j==1 ? 1 - i : i for j in 1:2, i in 0.0:step:1.0 ]
	inverted = [ 1 .- channelview( Gray.( img ) ) for img in images ]
	dists    = [ img/sum( img ) for img in inverted ]
	barys = wassBarycenters2d( dists, weights, 0.0005, l22d, 100, sharpen )
	imgs = [ Gray.(1 .- bary./maximum( bary )) for bary in barys ]
	return imgs
end

function saveImages( folderName, images )
	mkdir( folderName )
	for (i, img) in enumerate( images )
		save( folderName*"/$i.png", img )
	end
end

max_stars_repo_name: undefdev/barycuda

max_issues_repo_name: undefdev/barycuda

max_forks_repo_name: undefdev/barycuda