text:
include("Continuum.jl")
export FiniteStrainContinuum

"""
    FiniteStrainContinuum

Constructs a finite strain element.

- `eledim`: spatial dimension of the element (default = 2).
- `mat`: constitutive law, a length `#elem` vector of materials such as [`PlaneStress`](@ref)
- `elnodes`: the node indices in this finite element, an integer array
- `props`: property dictionary
- `coords`: coordinates of the vertices of the element
- `dhdx`, `weights`, `hs`: data for integral
- `stress`: stress at each quadrature points

# Example
```julia
#   Local degrees of freedom
#   4 ---- 3
#
#   1 ---- 2

nx = 10
ny = 5
h = 0.1
element = FiniteStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
for j = 1:ny
    for i = 1:nx
        n = (nx+1)*(j-1) + (i-1)+1
        elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        ngp = 3 # 3 x 3 Gauss points per element
        coords = [(i-1)*h (j-1)*h
                    i*h (j-1)*h
                    i*h j*h
                    (i-1)*h j*h]
        push!(element, FiniteStrainContinuum(coords,elnodes, prop, ngp))
    end
end
```
"""
mutable struct FiniteStrainContinuum <: Continuum
    eledim::Int64
    mat  # constitutive law
    elnodes::Array{Int64}   # the node indices in this finite element
    props::Dict{String, Any}
    coords::Array{Float64}
    dhdx::Array{Array{Float64}}  # 4nPointsx2 matrix
    weights::Array{Float64}
    hs::Array{Array{Float64}}
    stress::Array{Array{Float64}}
end

Base.show(io::IO, z::FiniteStrainContinuum) = print(io, "FiniteStrainContinuum with $(length(z.weights)) Gauss quadrature points and $(z.props["name"]) material.")

"""
    FiniteStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
"""
function FiniteStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
    eledim = 2
    dhdx, weights, hs = get2DElemShapeData( coords, ngp )
    nGauss = length(weights)
    name = props["name"]
    if name=="PlaneStrain"
        mat = [PlaneStrain(props) for i = 1:nGauss]
    elseif name=="Scalar1D"
        mat = [Scalar1D(props) for i = 1:nGauss]
    elseif name=="PlaneStress"
        mat = [PlaneStress(props) for i = 1:nGauss]
    elseif name=="PlaneStressPlasticity"
        mat = [PlaneStressPlasticity(props) for i = 1:nGauss]
    elseif name=="PlaneStrainViscoelasticityProny"
        mat = [PlaneStrainViscoelasticityProny(props) for i = 1:nGauss]
    elseif name=="PlaneStressViscoelasticityProny"
        mat = [PlaneStressViscoelasticityProny(props) for i = 1:nGauss]
    elseif name=="PlaneStressIncompressibleRivlinSaunders"
        mat = [PlaneStressIncompressibleRivlinSaunders(props) for i = 1:nGauss]
    elseif name=="NeuralNetwork2D"
        mat = [NeuralNetwork2D(props) for i = 1:nGauss]
    else
        error("Not implemented yet: $name")
    end
    strain = Array{Array{Float64}}(undef, length(weights))
    FiniteStrainContinuum(eledim, mat, elnodes, props, coords, dhdx, weights, hs, strain)
end

function getStiffAndForce(self::FiniteStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]
    Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:2*nnodes]
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        E = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy]
        DE = [Dux+0.5*(Dux*Dux+Dvx*Dvx); Dvy+0.5*(Duy*Duy+Dvy*Dvy); Duy+Dvx+Dux*Duy+Dvx*Dvy]

        # #@show "+++",E, DE
        S, dS_dE = getStress(self.mat[k], E, DE, Δt)
        # error()
        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8

        S∂∂E∂∂u = [g1*g1'*S[1]+g2*g2'*S[2]+g1*g2'*S[3]+g2*g1'*S[3] zeros(nnodes, nnodes);
                  zeros(nnodes, nnodes)  g1*g1'*S[1]+g2*g2'*S[2]+g1*g2'*S[3]+g2*g1'*S[3]]

        stiff += (∂E∂u * dS_dE * ∂E∂u' + S∂∂E∂∂u)*self.weights[k] # 8x8
    end
    return fint, stiff
end

function getInternalForce(self::FiniteStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    n = dofCount(self)
    fint = zeros(Float64,n)
    out = Array{Float64}[]
    nnodes = length(self.elnodes)
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]
    Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:2*nnodes]
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]
        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        E = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy]
        DE = [Dux+0.5*(Dux*Dux+Dvx*Dvx); Dvy+0.5*(Duy*Duy+Dvy*Dvy); Duy+Dvx+Dux*Duy+Dvx*Dvy]
        S, _ = getStress(self.mat[k],E, DE, Δt)

        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8
    end
    return fint
end

function getStrain(self::FiniteStrainContinuum, state::Array{Float64})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    nGauss = length(self.weights)
    E = zeros(nGauss, 3)
    w∂E∂u = zeros(nGauss, 2nnodes, 3)
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]  # 8x3
        E[k,:] = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy] # 3x1
        w∂E∂u[k,:,:] = ∂E∂u*self.weights[k] # 8x3
    end

    return E, w∂E∂u
end

function  getForceAndDforceDstress(self::FiniteStrainContinuum, state::Array{Float64},
    stress::Array{Float64,2})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    nStrain = 3
    nGauss = length(self.weights)
    fint = zeros(Float64, ndofs)
    dfint_dstress = zeros(Float64,  ndofs, nGauss * nStrain)

    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        # compute  ∂E∂u, 3 by 2nnodes array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        S = stress[k, :]

        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8

        dfint_dstress[:, (k-1)*nStrain+1:k*nStrain] = ∂E∂u * self.weights[k]

    end
    return fint , dfint_dstress
end


text:
# Something like `DataFrames.combine`
# https://juliadata.github.io/DataFrames.jl/stable/man/split_apply_combine/
# https://juliadata.github.io/DataFrames.jl/stable/lib/functions/#DataFrames.select

"""
    oncol(iname₁ => spec₁, ..., inameₙ => specₙ) -> f::Function
    oncol(; \$iname₁ = spec₁, ..., \$inameₙ = specₙ) -> f::Function

Combine functions that work on a column and create a function that
work on an entire row.

It constructs a reducing step function acting on a table row where
`specᵢ` is either a reducing step function or a `Pair` of a reducing
step function and an output column name.

It also defines a unary function when `specᵢ` is either a unary
function or a `Pair` of a unary function and an output column name.

This function is inspired by the "`Pair` notation" in DataFrames.jl
(see also [Split-apply-combine ·
DataFrames.jl](https://juliadata.github.io/DataFrames.jl/stable/man/split_apply_combine/)
and
[`DataFrames.select`](https://juliadata.github.io/DataFrames.jl/stable/lib/functions/#DataFrames.select)).

# Examples
```jldoctest oncol
julia> using DataTools
       using Transducers

julia> rf = oncol(a = +, b = *);

julia> foldl(rf, Map(identity), [(a = 1, b = 2), (a = 3, b = 4)])
(a = 4, b = 8)

julia> rf((a = 1, b = 2), (a = 3, b = 4))
(a = 4, b = 8)

julia> rf = oncol(:a => (+) => :sum, :a => max => :max);

julia> foldl(rf, Map(identity), [(a = 1,), (a = 2,)])
(sum = 3, max = 2)

julia> rf((sum = 1, max = 1), (a = 2,))
(sum = 3, max = 2)

julia> rf = oncol(:a => min, :a => max);

julia> foldl(rf, Map(identity), [(a = 2,), (a = 1,)])
(a_min = 1, a_max = 2)

julia> rf((a_min = 2, a_max = 2), (a = 1,))
(a_min = 1, a_max = 2)

julia> foldl(rf, Map(x -> (a = x,)), [5, 2, 6, 8, 3])
(a_min = 2, a_max = 8)
```

`oncol` also defines a unary function

```jldoctest oncol
julia> f = oncol(a = string);

julia> f((a = 1, b = 2))
(a = "1",)
```

Note that `oncol` does not verify the arity of input functions.  If
the input functions have unary and binary methods, `oncol` is callable
with both arities:

```jldoctest oncol
julia> f((a = 1, b = 2), (a = 3, b = 4))
(a = "13",)
```
"""
oncol

struct Property{name} end
Property(p::Property) = p
Property(name::Symbol) = Property{name}()
Property(::Val{name}) where {name} = Property{name}()

@inline getprop(x, ::Property{name}) where {name} = getproperty(x, name)
@inline Base.Symbol(::Property{name}) where {name} = name::Symbol

const PropertyLike = Union{Symbol,Val,Property}

struct OnRowFunction{FS} <: Function
    functions::FS
end

# :x => f                    (x_f = f(a.x, b.x),)
# :x => f => :y              (y = f(a.x, b.x),)
# (:x, :y) => f => :z        (z = f((a.x, a.y), (b.x, b.y)),)   ???

@inline (f::OnRowFunction)(x) =
    mapfoldl(merge, f.functions; init = NamedTuple()) do (iname, g, oname)
        Base.@_inline_meta
        (; Symbol(oname) => g(getprop(x, iname)))
    end

@inline (rf::OnRowFunction)(acc, x) = next(rf, acc, x)
@inline Transducers.next(rf::OnRowFunction, acc, x) =
    mapfoldl(merge, rf.functions; init = NamedTuple()) do (iname, op, oname)
        Base.@_inline_meta
        (; Symbol(oname) => next(op, getprop(acc, oname), getprop(x, iname)))
    end

Transducers.start(rf::OnRowFunction, init) =
    mapfoldl(merge, rf.functions; init = NamedTuple()) do (_, op, oname)
        (; Symbol(oname) => start(op, init))
    end

# TODO: dispatch on "Initializer" type instead
Transducers.start(rf::OnRowFunction, init::RowLike) =
    mapfoldl(merge, rf.functions; init = NamedTuple()) do (_, op, oname)
        (; Symbol(oname) => start(op, getprop(init, oname)))
    end

Transducers.complete(rf::OnRowFunction, acc) =
    mapfoldl(merge, rf.functions; init = NamedTuple()) do (_, op, oname)
        (; Symbol(oname) => complete(op, getprop(acc, oname)))
    end

Transducers.combine(rf::OnRowFunction, a, b) =
    mapfoldl(merge, rf.functions; init = NamedTuple()) do (_, op, oname)
        (; Symbol(oname) => combine(op, getprop(a, oname), getprop(b, oname)))
    end

# TODO: define better API:
Transducers._asmonoid(rf::OnRowFunction) =
    OnRowFunction(map(modifying(@optic(_[2]) => Transducers._asmonoid), rf.functions))
Transducers.Completing(rf::OnRowFunction) =
    OnRowFunction(map(modifying(@optic(_[2]) => Transducers.Completing), rf.functions))

@inline process_spec_kwargs((iname, spec)::Pair) =
    process_spec(Property(iname), spec, Property(iname))

@inline process_spec(name::PropertyLike) = (Property(name), identity, Property(name))
@inline process_spec((iname, spec)::Pair) = process_spec(iname, spec)
@inline process_spec(iname, (op, oname)::Pair) = process_spec(iname, op, oname)
@inline process_spec(iname, op) = process_spec(iname, op, outname(iname, op))
@inline process_spec(iname, op, oname) = (Property(iname), op, Property(oname))

outname(iname, op) = Symbol(Symbol(iname), '_', funname(op))

funname(op::ComposedFunction) = Symbol(funname(op.f), :_, funname(op.g))
function funname(op)
    name = string(op)
    return startswith(name, "#") ? :function : Symbol(op)
end

oncol(; specs...) =
    if any(((_, spec),) -> spec isa Pair, specs)
        fs = map(process_spec, Tuple(specs))
        @assert allunique(map(last, fs))  # TODO: uniquify function names
        OnRowFunction(fs)
    else
        OnRowFunction(map(process_spec_kwargs, Tuple(specs)))
    end

oncol(specs::Union{Pair{<:PropertyLike},PropertyLike}...) =
    OnRowFunction(map(process_spec, specs))

ext: jl


text:
"""
	$(SIGNATURES)

Legends are 1xn string matrices. Make that from a vector.
`nothing` is passed through.
"""
make_legend(v) = string.(reshape(v, 1, length(v)));
make_legend(::Nothing) = nothing;

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


text:
module Domains

import ..Geometry: Geometry, float_type
using IntervalSets
export RectangleDomain

"""
    AbstractDomain

A domain represents a region of space.
"""
abstract type AbstractDomain end

const BCTagType = Union{Nothing, Tuple{Symbol, Symbol}}

float_type(domain::AbstractDomain) = float_type(coordinate_type(domain))

"""
    boundary_names(obj::Union{AbstractDomain, AbstractMesh, AbstractTopology})

A tuple or vector of unique boundary names of a spatial domain.
"""
function boundary_names end

struct IntervalDomain{CT, B} <: AbstractDomain where {
    CT <: Geometry.Abstract1DPoint{FT},
    B <: BCTagType,
} where {FT}
    coord_min::CT
    coord_max::CT
    boundary_names::B
end

isperiodic(domain::IntervalDomain) = isnothing(domain.boundary_names)
boundary_names(domain::IntervalDomain) =
    isperiodic(domain) ? () : unique(domain.boundary_names)

"""
    IntervalDomain(coord⁻, coord⁺; periodic=true)
    IntervalDomain(coord⁻, coord⁺; boundary_names::Tuple{Symbol,Symbol})

Construct a `IntervalDomain`, the closed interval is given by `coord⁻`, `coord⁺` coordinate arguments.

Either a `periodic` or `boundary_names` keyword argument is required.
"""
function IntervalDomain(
    coord_min::Geometry.Abstract1DPoint,
    coord_max::Geometry.Abstract1DPoint;
    periodic = false,
    boundary_tags = nothing, # TODO: deprecate this
    boundary_names::BCTagType = boundary_tags,
)
    if !periodic && isnothing(boundary_names)
        throw(
            ArgumentError(
                "if `periodic=false` then an `boundary_names::Tuple{Symbol,Symbol}` keyword argument is required.",
            ),
        )
    end
    IntervalDomain(promote(coord_min, coord_max)..., boundary_names)
end
IntervalDomain(coords::ClosedInterval; kwargs...) =
    IntervalDomain(coords.left, coords.right; kwargs...)

coordinate_type(::IntervalDomain{CT}) where {CT} = CT
Base.eltype(domain::IntervalDomain) = coordinate_type(domain)

function Base.show(io::IO, domain::IntervalDomain)
    print(
        io,
        nameof(typeof(domain)),
        "($(domain.coord_min), $(domain.coord_max); ",
    )
    if isperiodic(domain)
        print(io, "periodic=true)")
    else
        print(io, "boundary_names = $(domain.boundary_names))")
    end
end

struct RectangleDomain{I1 <: IntervalDomain, I2 <: IntervalDomain} <:
       AbstractDomain
    interval1::I1
    interval2::I2
end
Base.:*(interval1::IntervalDomain, interval2::IntervalDomain) =
    RectangleDomain(interval1, interval2)

boundary_names(domain::RectangleDomain) = unique([
    boundary_names(domain.interval1)...,
    boundary_names(domain.interval2)...,
])

"""
    RectangleDomain(x1::ClosedInterval, x2::ClosedInterval;
        x1boundary::Tuple{Symbol,Symbol},
        x2boundary::Tuple{Symbol,Symbol},
        x1periodic = false,
        x2periodic = false,
    )

Construct a `RectangularDomain` in the horizontal.
If a given x1 or x2 boundary is not periodic, then `x1boundary` or `x2boundary` boundary name keyword arguments must be supplied.
"""
function RectangleDomain(
    x1::ClosedInterval{X1CT},
    x2::ClosedInterval{X2CT};
    x1periodic = false,
    x2periodic = false,
    x1boundary::BCTagType = nothing,
    x2boundary::BCTagType = nothing,
) where {X1CT <: Geometry.Abstract1DPoint, X2CT <: Geometry.Abstract1DPoint}
    interval1 =
        IntervalDomain(x1; periodic = x1periodic, boundary_names = x1boundary)
    interval2 =
        IntervalDomain(x2; periodic = x2periodic, boundary_names = x2boundary)
    return interval1 * interval2
end

function Base.show(io::IO, domain::RectangleDomain)
    indent = get(io, :indent, 0)
    println(io, nameof(typeof(domain)), "(")
    println(io, " "^(indent + 2), domain.interval1)
    println(io, " "^(indent + 2), domain.interval2)
    print(io, " "^indent, ")")
end
coordinate_type(domain::RectangleDomain) = typeof(
    Geometry.product_coordinates(
        domain.interval1.coord_min,
        domain.interval2.coord_min,
    ),
)

"""
    SphereDomain(radius)

A domain representing the surface of a sphere with radius `radius`.
"""
struct SphereDomain{FT} <: AbstractDomain where {FT <: AbstractFloat}
    radius::FT
end
Base.show(io::IO, domain::SphereDomain) =
    print(io, nameof(typeof(domain)), "($(domain.radius))")

boundary_names(::SphereDomain) = ()
coordinate_type(::SphereDomain{FT}) where {FT} = Geometry.Cartesian123Point{FT}

end # module


text:
# Model
include_model("hopper3D")

function surf(x)
	m = 1.0

	if x[1] >= 0.0
		return m * x[1]
	else
		return -m * x[1]
	end
end

function surf_grad(x)
	m = 1.0

	if x[1] >= 0.0
		return [m; 0.0]
	else
		return [-m; 0.0]
	end
end

model = Hopper3D{Discrete, FixedTime}(n, m, d,
			mb, ml, Jb, Jl,
			μ, g,
			qL, qU,
            nq, nu, nc, nf, nb,
            idx_u,
            idx_λ,
            idx_b,
            idx_ψ,
            idx_η,
            idx_s,
			surf,
			surf_grad)

model_ft = free_time_model(model)

function ϕ_func(model::Hopper3D, q)
	k = kinematics(model, q)
    @SVector [k[3] - model.surf(k[1:2])]
end

# Horizon
T = 31

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[end] = 2.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -Inf
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
mrp_init = MRP(UnitQuaternion(RotZ(0.0) * RotY(0.0) * RotX(0.0)))

z_h = 0.0
# q1 = [0.0, 0.0, 0.5, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
# ϕ_func(model, q1)
# x1 = [q1; q1]
_qT = [0.5, 0.0, 0.0, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]

zh = -1.0 * ϕ_func(model_ft, _qT)[1]

q1 = [-0.5, 0.0, zh, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
ϕ_func(model_ft, q1)
x1 = [q1; q1]

qT = [0.5, 0.0, zh, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
ϕ_func(model_ft, qT)
xT = [qT; qT]

xl, xu = state_bounds(model_ft, T,
    [model_ft.qL; model_ft.qL], [model_ft.qU; model_ft.qU])#,
    # x1 = x1, xT = xT)

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])

qp = [1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(0.0 * [qp; qp]) for t = 1:T],
    [Diagonal([1.0e-2, 1.0e-2, 1.0e-2,
		1.0e-3 * ones(model_ft.nc)..., 1.0e-3 * ones(model_ft.nb)...,
		zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb - 1)..., 0.0])
		for t = 1:T-1],
    [[qT; qT] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    1.0)

obj_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)

obj_velocity = velocity_objective(
    [Diagonal(1.0e-1 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3]))

obj = MultiObjective([obj_tracking, obj_penalty, obj_penalty])

# Constraints
include_constraints(["free_time", "contact", "stage", "loop"])
con_free_time = free_time_constraints(T)
con_contact = contact_constraints(model_ft, T)

p1_ref = kinematics(model_ft, q1)
pT_ref = kinematics(model_ft, qT)

function pinned1!(c, x, u, t)
    q = view(x, 1:7)
    c[1:3] = p1_ref - kinematics(model_ft, q)
	nothing
end

function pinnedT!(c, x, u, t)
    q = view(x, 7 .+ (1:7))
	c[1:3] = pT_ref - kinematics(model_ft, q)
	nothing
end

T_fix = 5
n_stage = 3
t_idx1 = vcat([t for t = 1:T_fix]...)
t_idxT = vcat([(T - T_fix + 1):T]...)

con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
con_pinnedT = stage_constraints(pinnedT!, n_stage, (1:0), t_idxT)

con_loop = loop_constraints(model, collect([(2:7)...,(9:14)...]), 1, T)

con = multiple_constraints([con_free_time,
	con_contact,
	con_pinned1,
	con_pinnedT,
	con_loop])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
q_ref = linear_interpolation(q1, qT, T+1)
x0 = configuration_to_state(q_ref) # linear interpolation on state
u0 = [[1.0e-3 * rand(model.m); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time z̄, info = solve(prob, copy(z0),
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
tf, t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
h̄ = mean(h̄)

# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# render(vis)
visualize!(vis, model_ft,
	# state_to_configuration(x̄),
	q,
	Δt = h̄)

# f(x) = x[3] - model_ft.surf(x[1:2])[1]
#
# xlims = [-1.0, 1.0]
# ylims = [-2.0, 2.0]
# mesh = GeometryBasics.Mesh(f,
# 	HyperRectangle(Vec(xlims[1], ylims[1], -2.0), Vec(xlims[2]-xlims[1], ylims[2]-ylims[1], 4.0)),
#     Meshing.MarchingCubes(), samples=(200, 200, Int(floor(200/8))))
# setobject!(vis["surface"], mesh,
# 	MeshPhongMaterial(color=RGBA{Float32}(1.0, 0.0, 0.0, 0.5)))

ref = Diagonal([-1.0; 1.0; 1.0])
mrp_init = MRP(RotZ(0.0) * RotY(-0.25 * π) * RotX(-0.0 * π))
# mrp_mirror = MRP(mrp_init')
mrp_mirror = MRP(mrp_init * ref)
q1 = [0.0, 0.0, 1.0, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
qT = [0.0, 0.0, 1.0, mrp_mirror.x, mrp_mirror.y, mrp_mirror.z, 0.5]

visualize!(vis, model_ft,
	[q1],
	# qa,
	# q_ref,
	Δt = h̄)

qa = linear_interpolation(q1, qT, 10)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using BehaviorTree
using AbstractTrees
using Test

function doSuccess()
    :success
end
function isTrue()
    :success
end
function doSuccess!()
    :success
end

function doFailure()
    :failure
end
function doRunning()
    :running
end
function isPositive(x)
    if x>0
        return :success
    end
    :failure
end

@testset "BehaviorTree.jl" begin
    tree = Sequence([doSuccess])
    @test tick(tree)[1] == :success
    @test tick(tree)[2] == [:success]
    tree = Sequence([doFailure])
    @test tick(tree)[1] == :failure
    tree = Sequence([doRunning])
    @test tick(tree)[1] == :running
    tree = Sequence([doSuccess, doFailure, doSuccess])
    @test tick(tree)[1] == :failure
    @test length([ c for c in children(tree)]) == 3
    tree = Selector([doFailure])
    @test tick(tree)[1] == :failure
    tree = Selector([doSuccess])
    @test tick(tree)[1] == :success
    tree = Selector([doRunning])
    @test tick(tree)[1] == :running
    tree = Selector([doFailure, doSuccess, doFailure])
    @test tick(tree)[1] == :success
    # nesting
    tree = Selector([doFailure, Sequence([doSuccess]), doFailure])
    @test tick(tree)[1] == :success
    @info tick(tree)
    @test tick(tree)[2] == [:failure,[:success]]
    # args
    tree = Sequence([isPositive])
    @test tick(tree, 1)[1] == :success
    @test tick(tree, -1)[1] == :failure
    tree = Selector([isPositive])
    @test tick(tree, 1)[1] == :success
    @test tick(tree, -1)[1] == :failure
end
@testset "dot" begin
    bt = Selector([
        doFailure,
        Sequence([isTrue, doFailure, doSuccess], "choice"),
        doRunning,
        doSuccess!
    ], "head")
    dot_graph = toDot(bt)
    #println(dot_graph)
    r = tick(bt)
    #png_graph =dot2png(dot_graph)

    bt = Selector([
        doFailure,
        Sequence([isTrue, doFailure, doSuccess], "choice"),
        Sequence([isTrue, doSuccess, doRunning], "choice2"),
        doSuccess!
    ], "head")
    dot_graph = toDot(bt)
    r = tick(bt)
    dot_graph_status = toDot(bt, r[2])
    println(dot_graph_status)
    png_graph = dot2png(dot_graph_status)
    filename= "test.png"
    open(filename, "w") do png_file
        write(png_file, png_graph)
    end
end


text:
"""
    A20{N} <: AbstractContinuousPost

Implementation of the reachability method for large linear systems with uncertain
inputs in the Krylov subspace from [[A20]](@ref).

### Fields

- `δ`          -- step-size of the discretization
- `max_order`  -- (optional, default: `5`) maximum zonotope order

### References

See [[A20]](@ref) and references therein.
"""
struct A20{N} <: AbstractContinuousPost
    δ::N
    max_order::Int
end

# convenience constructor using keywords
function A20(; δ::N, max_order::Int=5) where {N}
    return A20(δ, max_order)
end

step_size(alg::A20) = alg.δ
numtype(::A20{N}) where {N} = N

function rsetrep(::A20{N}) where {N}
    VT = Vector{N}
    MT = Matrix{N}
    RT = ReachSet{N, Zonotope{N, VT, MT}}
end

include("post.jl")
include("reach_homog.jl")
include("reach_inhomog.jl")


text:
# Rule mixins

@chain columns @default_kw @flattenable @bounds @description

@mix @columns struct Probabilistic{PT}
    # Field            | Default | Flatn | Limits     | Description
    prob_threshold::PT | 0.1     | true  | (0.0, 1.0) | "A real number between one and zero"
end

@mix @columns struct Timestep{TS,S}
    timestep::TS       | nothing | false | _          | "Timestep converted from sim data. Needs to be separate from rate for DateTime"
    nsteps::S          | 1.0     | false | _          | "The exact nsteps timestep, updated by precalcrule"
end

precalctimestep(rule, data) = precalctimestep(rule.timestep, rule, data)
precalctimestep(ruletimestep::DatePeriod, rule, data) =
    @set rule.nsteps = currenttimestep(data) / Millisecond(ruletimestep)
precalctimestep(ruletimestep::Nothing, rule, data) = @set rule.nsteps = 1
precalctimestep(ruletimestep, rule, data) =
    @set rule.nsteps = timestep(data) / ruletimestep


text:
###
# Example translated from the original necpp package:
# https://github.com/tmolteno/necpp/blob/master/example/test.py

using NECPP

function frequency_response()
    # Scan through frequencies from 1 to 30 MHz
    for f = 1.0:30.0
        nec = nec_create()
        handle_nec(nec_wire(nec, 1, 17, 0.0, 0.0, 2.0, 0.0, 0.0, 11.0, 0.1, 1.0, 1.0))
        handle_nec(nec_geometry_complete(nec, 1))
        handle_nec(nec_gn_card(nec, 1, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
        handle_nec(nec_fr_card(nec, 0, 1, f, 0.0))
        handle_nec(nec_ex_card(nec, 0, 0, 5, 0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0))
        handle_nec(nec_rp_card(nec, 0, 90, 1, 0, 5, 0, 0, 0.0, 90.0, 1.0, 0.0, 0.0, 0.0))
        result_index = 0
        z = Complex(nec_impedance_real(nec,result_index),
            nec_impedance_imag(nec,result_index))
        @printf "f=%0.2fMHz \t(%6.1f,%+6.1fI) Ohms\r\n" f z.re z.im
        nec_delete(nec)
    end
end

#frequency_response()


text:
"""
    struct BuffersFaceIntegral

Stores all temporary buffers needed during face integrals.
Avoids costly re-allocation.
"""
struct BuffersFaceIntegral
    dofsface::Array{Float64,2}
    dofsfaceneigh::Array{Float64,2}
    fluxface::Array{Float64,2}
    fluxfaceneigh::Array{Float64,2}
    numericalflux::Array{Float64,2}
    numericalflux_scaled::Array{Float64,2}

    function BuffersFaceIntegral(basis, ndofs)
        @assert basis.dimensions == 2
        basissize_1d = size(basis, 1)
        basissize_nd = length(basis)
        dofsface = zeros(basissize_1d,ndofs)
        dofsfaceneigh = similar(dofsface)
        fluxface = zeros(basissize_1d * 2,ndofs)
        fluxfaceneigh = similar(fluxface)
        numericalflux = similar(dofsface)
        numericalflux_scaled = similar(numericalflux)

        new(dofsface, dofsfaceneigh, fluxface, fluxfaceneigh, numericalflux, numericalflux_scaled)
    end
end

"""
    rusanov(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)

Computes the Rusanov (or local Lax-Friedrichs) numerical flux for degrees of freedom
`dofs`, degrees of freedoms of the neighbor `dofsneigh`, flux `flux`, flux of neighbor
`fluxneigh`, cellsize `dx`.
All quantities are assumed to be represented by a basis on the reference line.
The face is parametrized by a index `normalidx`, where 1 stands for a face in x-direction
and 2 for a face in y-direction.
The sign of the outer normal of the face is given by `normalsign`.

The numerical flux is stored in `numericalflux`.
Method also returns the maximal eigenvalue.
"""

function rusanov(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)
    maxeigenval_left = max_eigenval(eq, dofs, normalidx)
    maxeigenval_right = max_eigenval(eq, dofsneigh, normalidx)
    maxeigenval = max(maxeigenval_left, maxeigenval_right)

    basissize_1d = size(numericalflux)[1]
    for i=1:size(dofs)[2]
        #check which value of flux need to use and also the direction of the normal
        for j = 1:basissize_1d
            if normalidx == 1 #For faces in the X direction that is left and right
            #@info dofs[i]
                numericalflux[j,i] = ( 1/2 * ( flux[j,i] + fluxneigh[j,i] )   ) +
                                                 ( maxeigenval* normalsign * 1/2 * (dofs[j,i] - dofsneigh[j,i]) )

            else #For faces in the Y direction that is Top and Bottom
                numericalflux[j,i] = ( 1/2 * ( flux[j+basissize_1d,i] + fluxneigh[j+basissize_1d,i] )   ) +
                                                 ( maxeigenval* normalsign * 1/2 * (dofs[j,i] - dofsneigh[j,i]) )
            end

            if eq == Acoustic()
                if i > 3
                    numericalflux[j, i] = 0  # For acoustic, material flux to 0
                end
            end
        end
    end

    #numericalflux =  (1/dx) * ( ( 1/2 * ( flux + fluxneigh ) )  +  ( maxeigenval* normalsign * (dofs - dofsneigh ) * 1/2 ) )
    maxeigenval # maximum eigenvalue
end

"""
    exactRiemann(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)

Computes the exact flux for degrees of freedom `dofs`,
degrees of freedoms of the neighbor `dofsneigh`, flux `flux`,
flux of neighbor `fluxneigh`, cellsize `dx`.
All quantities are assumed to be represented by a basis on the reference line.
The face is parametrized by a index `normalidx`, where 1 stands for a face in x-direction
and 2 for a face in y-direction.
The sign of the outer normal of the face is given by `normalsign`.

The flux is stored in `numericalflux`.
Method also returns the maximal eigenvalue.
"""

function exactRiemann(eq, dofs, dofsneigh, flux, fluxneigh, normalidx, normalsign, numericalflux)
    maxeigenval_left = max_eigenval(eq, dofs, normalidx)
    maxeigenval_right = max_eigenval(eq, dofsneigh, normalidx)
    maxeigenval = max(maxeigenval_left, maxeigenval_right)

    # constants relevant to elastic equation
    λ = 2
    μ = 1
    ρ = 1
    cp = sqrt((λ + 2*μ)/ρ) # speed p-waves
    cs = sqrt(μ/ρ) #speed s-waves

    basissize_1d = size(numericalflux)[1]
    avr_flux = 1/2 * (flux + fluxneigh)
    dif_flux = 1/2 * (flux - fluxneigh)
    for j = 1:basissize_1d
        if normalidx == 1 #For faces in the X direction that is left and right
            numericalflux[j,1] = avr_flux[j,1] - cp*ρ * normalsign * dif_flux[j,4]
            numericalflux[j,2] = (avr_flux[j,2] - normalsign * dif_flux[j,2]
                                    - λ/(cp) * normalsign * dif_flux[j,4]
                                    + λ/(cp^2*ρ) * normalsign * dif_flux[j,1])
            numericalflux[j,3] = avr_flux[j,3] - cs*ρ * normalsign * dif_flux[j,5]
            numericalflux[j,4] = avr_flux[j,4] - 1/(cp*ρ) * normalsign * dif_flux[j,1]
            numericalflux[j,5] = avr_flux[j,5] - 1/(cs*ρ) * normalsign * dif_flux[j,3]
        else #For faces in the Y direction that is Top and Bottom
            numericalflux[j,1] = (avr_flux[j+basissize_1d,1] - normalsign * dif_flux[j+basissize_1d,1]
                                    - λ/cp * normalsign * dif_flux[j+basissize_1d,5]
                                    + λ/(cp^2*ρ) * normalsign * dif_flux[j+basissize_1d,2])
            numericalflux[j,2] = avr_flux[j+basissize_1d,2] - cp*ρ* normalsign * dif_flux[j+basissize_1d,5]
            numericalflux[j,3] = avr_flux[j+basissize_1d,3] - cs*ρ* normalsign * dif_flux[j+basissize_1d,4]
            numericalflux[j,4] = avr_flux[j+basissize_1d,4] - 1/(cs*ρ)* normalsign * dif_flux[j+basissize_1d,3]
            numericalflux[j,5] = avr_flux[j+basissize_1d,5] - 1/(cp*ρ)* normalsign * dif_flux[j+basissize_1d,2]
        end
    end

    maxeigenval # maximum eigenvalue
end

"""
    project_to_faces(globals, dofs, flux, dofsface, fluxface, face)

Projects degrees of freedom `dofs` and `flux` to `face`.
Result is stored in `dofsface` and `fluxface`.
Projection matrices are stored in `globals`.
"""

function project_to_faces(globals, dofs, flux, dofsface, fluxface, face)
    #dofsface .= dofs
    i = globals.basissize_1d
    dofsface .= globals.project_dofs_to_face[face] * dofs
    fluxface[1:i,:] = globals.project_flux_to_face[face] * flux[1:i*i,:]
    fluxface[i+1:end,:] = globals.project_flux_to_face[face] * flux[i*i+1:end,:]
    #dofs = 4x3  and flux = 8x3
    #face proj matrix = 2x4
    #order = 1
    #fluxface .= flux
end

"""
    evaluate_face_integral(eq, globals, buffers, cell, face, celldu)

Computes face integrals for equation `eq`, cell `cell` on face `face`.
Global matrices are passed in `globals`, buffers in `buffers`.
Result is stored in `celldu`. `solver` is specified in config file.
"""

function evaluate_face_integral(eq, globals, buffers, cell, face, celldu, solver)
    normalidx = globals.normalidxs[face]
    normalsign = globals.normalsigns[face]

    # Compute Riemann solver (in normal)
    buffers.numericalflux .= 0

    if solver == "exact"
        maxeigenval = exactRiemann(eq, buffers.dofsface, buffers.dofsfaceneigh, buffers.fluxface, buffers.fluxfaceneigh, normalidx, normalsign, buffers.numericalflux)
    else
        maxeigenval = rusanov(eq, buffers.dofsface, buffers.dofsfaceneigh, buffers.fluxface, buffers.fluxfaceneigh, cell.size[1], normalidx, normalsign, buffers.numericalflux)
    end

    #scale numerical flux with quadweights.
    buffers.numericalflux_scaled .= globals.quadweights_nd.*buffers.numericalflux*area(cell)

    celldu .= celldu .+ -normalsign*globals.project_dofs_from_face[face] * buffers.numericalflux_scaled

    maxeigenval
end


text:
@ghdef mutable struct GitCommit
    sha::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    author::Union{Dict, Nothing}
    commiter::Union{Dict, Nothing}
    message::Union{String, Nothing}
    tree::Union{Dict, Nothing}
    parents::Union{Vector, Nothing}
    verification::Union{Dict, Nothing}
end

namefield(gitcommit::GitCommit) = gitcommit.sha

@api_default function gitcommit(api::GitHubAPI, repo, commit_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/git/commits/$(name(commit_obj))"; options...)
    return GitCommit(result)
end

@api_default function create_gitcommit(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/git/commits"; options...)
    return GitCommit(result)
end


text:
"""
This module defines the following data types:
* `Geom` : acquisition geometry, i.e., positions of supersources, sources and receivers
* `Src` : source related acquisition parameters, e.g., source wavelet
It also provides methods that either does operations on these data type or
help their construction.
"""
module Acquisition

import GeoPhyInv.Models
import GeoPhyInv.Utils
using Distributions
using DataFrames
using Random
using SparseArrays
using CSV

"""
Acquisiton has supersources, sources and receivers.
Each supersource has `ns` multiple sources that are
injected (or active) simultaneously.
For each supersource,
a set of `nr` receivers are
recording waves.

# Fields

* `sx::Vector{Vector{Float64,1},1}` : ``x`` positions of sources
* `sz::Vector{Vector{Float64,1},1}` : ``z`` positions of sources
* `rx::Vector{Vector{Float64,1},1}` : ``x`` positions of receivers
* `rz::Vector{Vector{Float64,1},1}` : ``z`` positions of receivers
* `nss::Int64` : number of supersources
* `ns::Vector{Int64,1}` : number of sources for every supersource
* `nr::Vector{Int64,1}` : number of receivers for every supersource
"""
mutable struct Geom
	sx::Vector{Vector{Float64}}
	sz::Vector{Vector{Float64}}
	rx::Vector{Vector{Float64}}
	rz::Vector{Vector{Float64}}
	nss::Int64 # change to +ve only integer later
	ns::Vector{Int64} # change to +ve only integer later
	nr::Vector{Int64} # change to +ve only integer later
	"adding conditions that are to be false while construction"
	Geom(sx, sz, rx, rz, nss, ns, nr) =
		any([
                  broadcast(length, sx) != ns,
                  broadcast(length, sz) != ns,
                  broadcast(length, rx) != nr,
                  broadcast(length, rz) != nr,
       		  length(sx) != nss,
       		  length(sz) != nss,
       		  length(rz) != nss,
       		  length(rx) != nss
		  ]) ?
		error("Geom construct") : new(sx, sz, rx, rz, nss, ns, nr)
end # type

"Compare if two `TD`'s  are equal"
function Base.isequal(geom1::Geom, geom2::Geom, attrib=:all)
	srcfnames = [:sx, :sz, :nss, :ns]
	recfnames = [:rx, :rz, :nr, :nss]
	srcvec=[(isequal(getfield(geom1, name),getfield(geom2, name))) for name in srcfnames]
	recvec=[(isequal(getfield(geom1, name),getfield(geom2, name))) for name in recfnames]
	if(attrib == :sources)
		return all(srcvec)
	elseif(attrib == :receivers)
		return all(recvec)
	else
		return all(vcat(srcvec, recvec))
	end
end

"""
Print information about `Geom`
"""
function Base.print(geom::Geom, name::String="")
	println("\tAcquisition Geometry:\t",name)
	println("\t> number of supersources:\t",geom.nss)
	println("\t> sources per supersource:\t","min\t",minimum(geom.ns[:]), "\tmax\t", maximum(geom.ns[:]))
	println("\t> receivers per supersource:\t","min\t",minimum(geom.nr[:]), "\tmax\t", maximum(geom.nr[:]))
	println("\t> number of unique positions of:\t","sources\t",Geom_get([geom],:nus), "\treceivers\t", Geom_get([geom],:nur))
end

"""
Return some derived fields of `Geom`

# Arguments

* `acq::Vector{Geom}` : a vector of `Geom`
* `attrib::Symbol` : attribute to determine the return object
  * `=:nus` number of unique source positions in acquisition
  * `=:nur` number of unique receiver positions in acquisition
  * `=:uspos` a tuple of x and z positions of all the unique sources
  * `=:urpos` a tuple of x and z position of all the unique receivers
  * `=:geomurpos` a `Geom` vector as if all the unique receiver positions are used for each supersource
  * `=:geomuspos` a `Geom` vector as if all the unique source positions are used for each supersource
"""
function Geom_get(acq::Vector{Geom}, attrib::Symbol)
	ngeom = length(acq)
	sposx = Geom_getvec(acq,:sx);	sposz = Geom_getvec(acq,:sz)
	isequal(length(sposx), length(sposz)) ?
		spos = [[sposz[is],sposx[is]] for is in 1:length(sposx)] : error("input acq corrupt")

	rposx = Geom_getvec(acq,:rx);	rposz = Geom_getvec(acq,:rz)
	isequal(length(rposx), length(rposz)) ?
		rpos = [[rposz[ir],rposx[ir]] for ir in 1:length(rposx)] : error("input acq corrupt")

	uspos = unique(spos); nus=length(uspos)
	urpos = unique(rpos); nur=length(urpos)
	uspost = ([uspos[iu][1] for iu in 1:nus], [uspos[iu][2] for iu in 1:nus]);
        urpost = ([urpos[iu][1] for iu in 1:nur], [urpos[iu][2] for iu in 1:nur]);

	if(attrib == :nus)
		return nus
	elseif(attrib == :nur)
		return nur
	elseif(attrib == :uspos)
		return uspost
	elseif(attrib == :urpos)
		return urpost
	elseif(attrib == :geomurpos)
		return [Geom(acq[igeom].sx, acq[igeom].sz,
	       fill(urpost[2],acq[igeom].nss),
	       fill(urpost[1],acq[igeom].nss), acq[igeom].nss, acq[igeom].ns,
	       fill(nur,acq[igeom].nss)) for igeom=1:ngeom]
	elseif(attrib == :geomuspos)
		return [Geom(fill(uspost[2],acq[igeom].nss),
	       fill(uspost[1],acq[igeom].nss),
	       acq[igeom].rx, acq[igeom].rz,
	       acq[igeom].nss, fill(nus,acq[igeom].nss),
	       acq[igeom].nr) for igeom=1:ngeom]
	else
		error("invalid attrib")
	end
end

"""
Given receiver positions `rpos` and `rpos0`.
Returns an array Int indices of the dimension of number of supersources
with `true` at indices, if the waves due to that particular source are
recorded.
"""
function Geom_find(acq::Geom; rpos::Array{Float64,1}=nothing, rpos0::Array{Float64,1}=nothing)
	rpos==nothing ? error("need rpos") : nothing
	sson = Array{Vector{Int64}}(undef,acq.nss);
	for iss = 1:acq.nss
		rvec = [[acq.rz[iss][ir],acq.rx[iss][ir]] for ir=1:acq.nr[iss]]
		ir=findfirst(rvec, rpos)
		if(rpos0==nothing)
			ir0=1;
		else
			ir0=findfirst(rvec, rpos0);
		end
		sson[iss] = ((ir != 0) && (ir0 != 0)) ? [ir,ir0] : [0]
	end
	return sson
end

#=
"""
Modify input `Geom` such that the output `Geom` has
either sources or receivers on the boundary of
`mgrid`.

# Arguments

* `acqgeom::Geom` : input geometry
* `mgridi` : grid to determine the boundary
* `attrib::Symbol` : decide return
  * `=:srcborder` sources on boundary (useful for back propagation)
  * `=:recborder` receivers on boundary
"""
function Geom_boundary(acqgeom::Geom,
	      mgrid
	      attrib::Symbol
	     )

	if((attrib == :recborder) | (attrib == :srcborder))
		bz, bx, nb = border(mgrid, 3, :inner)
	end
	"change the position of receivers to the boundary"
	if(attrib == :recborder)
		nr = nb; 		sx = acqgeom.sx;
		sz = acqgeom.sz;		ns = acqgeom.ns;
		nss = acqgeom.nss;
		return Geom(sx, sz, fill(bx,nss), fill(bz,nss), nss, ns, fill(nr,nss))
	"change the position of source (not supersources) to the boundary for back propagation"
	elseif(attrib == :srcborder)
		ns = nb;		rx = acqgeom.rx;
		rz = acqgeom.rz;		nr = acqgeom.nr;
		nss = acqgeom.nss;
		return Geom(fill(bx,nss), fill(bz,nss), rx, rz, nss, fill(ns,nss), nr)
	else
		error("invalid attrib")
	end
end
=#

"""
Check if all the sources and receivers in `Geom` are within the model

# Return

* `true` if all the positions within the model, `false` otherwise
"""
function Geom_check(geom::Geom, mgrid)
	xmin, zmin, xmax, zmax = mgrid[2][1], mgrid[1][1], mgrid[2][end], mgrid[1][end]

	checkvec=fill(false, geom.nss)
	for iss=1:geom.nss
		checkvec[iss] = any(
		      vcat(
		     (((geom.sx[iss] .- xmin).*(xmax .- geom.sx[iss])) .< 0.0),
		     (((geom.sz[iss] .- zmin).*(zmax .- geom.sz[iss])) .< 0.0),
		     (((geom.rx[iss] .- xmin).*(xmax .- geom.rx[iss])) .< 0.0),
		     (((geom.rz[iss] .- zmin).*(zmax .- geom.rz[iss])) .< 0.0),
		     ) )
	end
	return !(any(checkvec))
end

"""
Check if the input acquisition geometry is fixed spread.
"""
function Geom_isfixed(geom::Geom)
	isfixed=false
	# can have different source positions, but also different number of sources?
	for field in [:rx, :rz, :nr, :ns]
		f=getfield(geom, field)
		if(all(f .== f[1:1]))
			isfixed=true
		else
			isfixed=false
			return isfixed
		end
	end
	return isfixed
end

"""
A fixed spread acquisition has same set of sources and
receivers for each supersource.
This method constructs a
fixed spread acquisition geometry using either a
horizontal or vertical array of supersources/ receivers.
Current implementation has only one source for every supersource.

# Arguments

* `ssmin::Float64` : minimum coordinate for sources
* `ssmax::Float64` : maximum coordinate for sources
* `ss0::Float64` : consant coordinate for sources
* `rmin::Float64` : minimum coordinate for receivers
* `rmax::Float64` : maximum coordinate for receivers
* `r0::Float64` : consant coordinate for receivers
* `nss::Int64` : number of supersources
* `nr::Int64` : number of receivers
* `ssattrib::Symbol=:horizontal` : supersource array kind
  `=:vertical` : vertical array of supersources
  `=:horizontal` horizontal array of supersources
* `rattrib::Symbol=:horizontal` : receiver array kind
  `=:vertical` : vertical array of receivers
  `=:horizontal` horizontal array of receivers
* `rand_flags::Vector{Bool}=[false, false]` : decide placement of supersources and receivers
  `=[true, false]` : randomly place supersources for regularly spaced receivers
  `=[true, true]` : randomly place supersources and receivers
  `=[false, false]` : regularly spaced supersources and receivers
  `=[false, true]` : randomly place receivers for regularly spaced supersources

# Return
* a fixed spread acquisition geometry `Geom`
"""
function Geom_fixed(
	      ssmin::Real,
	      ssmax::Real,
	      ss0::Real,
	      rmin::Real,
	      rmax::Real,
	      r0::Real,
	      nss::Int64,
	      nr::Int64,
	      ssattrib::Symbol=:horizontal,
	      rattrib::Symbol=:horizontal,
 	      rand_flags::Vector{Bool}=[false, false];
	      ssα::Real=0.0,
	      rα::Real=0.0,
	      ns::Vector{Int64}=ones(Int,nss),
	      srad::Real=0.0
	     )

	ssmin=Float64(ssmin); rmin=Float64(rmin)
	ssmax=Float64(ssmax); rmax=Float64(rmax)
	ss0=Float64(ss0); r0=Float64(r0)
	ssα=Float64(ssα)*pi/180.
	rα=Float64(rα)*pi/180.

	ssarray = isequal(nss,1) ? fill(ssmin,1) : (rand_flags[1] ?
					     Random.rand(Uniform(minimum([ssmin,ssmax]),maximum([ssmin,ssmax])),nss) : range(ssmin,stop=ssmax,length=nss))
	if(ssattrib==:horizontal)
		ssz = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssx=ssarray
	elseif(ssattrib==:vertical)
		ssx = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssz=ssarray
	else
		error("invalid ssattrib")
	end

	rarray = isequal(nr,1) ? fill(rmin,1) : (rand_flags[2] ?
					  Random.rand(Uniform(minimum([rmin,rmax]),maximum([rmin,rmax])),nr) : range(rmin,stop=rmax,length=nr))
	if(rattrib==:horizontal)
		rz = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rx = rarray
	elseif(rattrib==:vertical)
		rx = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rz = rarray
	else
		error("invalid rattrib")
	end
	rxall = [rx for iss=1:nss];
	rzall = [rz for iss=1:nss];
	ssxall=[zeros(ns[iss]) for iss=1:nss];
	sszall=[zeros(ns[iss]) for iss=1:nss];
	for iss in 1:nss
		for is in 1:ns[iss]
			θ=Random.rand(Uniform(-Float64(pi),Float64(pi)))
			r = iszero(srad) ? 0.0 : Random.rand(Uniform(0,srad))
			x=r*cos(θ)
			z=r*sin(θ)
			ssxall[iss][is]=ssx[iss]+x
			sszall[iss][is]=ssz[iss]+z
		end
	end
	return Geom(ssxall, sszall, rxall, rzall, nss, ns, fill(nr,nss))
end

"""
Randomly place source and receivers in a given model
"""
function Geom_fixed(mod::Models.Seismic,
	      nss::Int64=5,
	      nr::Int64=5,
	      ns::Vector{Int64}=ones(Int,nss))

	xmin=mod.mgrid[2][1]
	xmax=mod.mgrid[2][end]
	zmin=mod.mgrid[1][1]
	zmax=mod.mgrid[1][end]

	rxall = [Random.rand(Uniform(xmin,xmax),nr) for iss=1:nss];
	rzall = [Random.rand(Uniform(xmin,xmax),nr) for iss=1:nss];
	ssxall = [Random.rand(Uniform(xmin,xmax),ns[iss]) for iss=1:nss];
	sszall = [Random.rand(Uniform(xmin,xmax),ns[iss]) for iss=1:nss];

	return Geom(ssxall, sszall, rxall, rzall, nss, ns, fill(nr,nss))
end

"""
Circular acquisition. The sources and receivers can also be placed on a circle of radius
`rad`. The origin of the circle is at `loc`.
This geometry is unrealistic, but useful for testing.
Receivers are placed such that the limits
of the angular offset are given by `θlim`

# Arguments

* `nss::Int64=10` : number of supersources
* `nr::Int64=10` : number receivers for each super source
* `loc::Vector{Float64}=[0.,0.]` : location of origin
* `rad::Vector{Float64}=[100.,100.]` : radius for source and receiver circles, for example,
  * `=[0.,100.]` for sources at the center of circle
* `θlim::Vector{Float64}=[0.,π]` : acquisition is limited to these angular offsets between 0 and π

# Return

* a circular acquisition geometry `Geom`
"""
function Geom_circ(;
		   nss::Int64=10,
		   nr::Int64=10,
		   loc::Vector=[0.,0.],
		   rad::Vector=[100.,100.],
		   θlim::Vector=[0.,π]
	       )

	# modify nr such that approximately nr receivers are present in θlim
	nra = nr * convert(Int64, floor(π/(abs(diff(θlim)[1]))))

	(nra==0) && error("θlim should be from 0 to π")

	sxx, szz = circ_coord(loc..., nss, rad[1])
	rxxa, rzza = circ_coord(loc..., nra, rad[2])
	sx=Array{Array{Float64}}(undef,nss)
	sz=Array{Array{Float64}}(undef,nss)
	rx=Array{Array{Float64}}(undef,nss)
	rz=Array{Array{Float64}}(undef,nss)
	ns=Array{Int64}(undef,nss)
	nr=Array{Int64}(undef,nss)

	for iss=1:nss
		rangles = [angle(complex(rxxa[ira]-loc[2], rzza[ira]-loc[1])) for ira in 1:nra]
		sangles = angle(complex(sxx[iss]-loc[2], szz[iss]-loc[1]))
		diff = abs.(rangles.-sangles)
		diff = [min(diff[ira], (2π - diff[ira])) for ira in 1:nra]
		angles = findall(minimum(θlim) .<= diff .<= maximum(θlim))
		sx[iss] = [sxx[iss]]
		sz[iss] = [szz[iss]]
		ns[iss] = 1
		nr[iss] = length(angles)
		rx[iss] = rxxa[angles]
		rz[iss] = rzza[angles]
	end
	geom =  Geom(sx, sz, rx, rz, nss, ns, nr)

	print(geom, "circular acquisition geometry")
	return geom
end

function circ_coord(z, x, n, rad)
	θ = (n==1) ? [0.0] : collect(range(0, stop=2π, length=n+1))
	xx = (rad .* cos.(θ) .+ x)[1:n]
	zz = (rad .* sin.(θ) .+ z)[1:n]
	return xx, zz
end

"""
Appends the input a vector of acquisition geometries.
"Merge" `Geom` objects in `geomv` to `geom1`.
Typically you have one acquisition,
you change it and want to create the macro acquisition
where the supersources of both acquisitions are appended.
"""
function Geom_add!(geom1::Geom, geomv::Vector{Geom})
	for igeom=1:length(geomv)
		append!(geom1.sx, geomv[igeom].sx)
		append!(geom1.sz, geomv[igeom].sz)
		append!(geom1.rx, geomv[igeom].rx)
		append!(geom1.rz, geomv[igeom].rz)
		append!(geom1.ns, geomv[igeom].ns)
		append!(geom1.nr, geomv[igeom].nr)
		geom1.nss += geomv[igeom].nss
	end
	return geom1
end

"""
Adds input positions as either sources or receivers of every supershot.
"""
function Geom_add!(geom::Geom, zpos::Vector{Float64}, xpos::Vector{Float64}, attrib::Symbol)
	length(xpos) == length(zpos) ? nothing : error("same length vectors needed")
	for iss=1:geom.nss
		if(attrib == :receivers)
			geom.rx[iss] = vcat(geom.rx[iss][:], xpos)
			geom.rz[iss] = vcat(geom.rz[iss][:], zpos)
			geom.nr[iss] += length(xpos)
		elseif(attrib == :sources)
			geom.sx[iss] = vcat(geom.sx[iss][:], xpos)
			geom.sz[iss] = vcat(geom.sz[iss][:], zpos)
			geom.ns[iss] += length(xpos)
		end
	end
end

"""
Advance either source or receiver array in an acquisition geometry in horizontal or vertical directions.

# Arguments
* `geom::Geom` : acquisition geometry that is updated
* `advances::Vector{Float64}=[[0.,0.], [0.,0.,]]` : source and receiver advancements
"""
function Geom_advance(geom::Geom, advances::Vector{Vector{Float64}}=[[0.,0.], [0.,0.]])
	geom_out=deepcopy(geom)
	for iss=1:geom.nss
		geom_out.sx[iss][:] .+= advances[1][2]
		geom_out.sz[iss][:] .+= advances[1][1]
		geom_out.rx[iss][:] .+= advances[2][2]
		geom_out.rz[iss][:] .+= advances[2][1]
	end
	return geom_out
end

"""
return a vector of the order
"""
function Geom_getvec(geom::Vector{Geom}, field::Symbol)
	np = length(geom);
	vect = [getfield(geom[ip],field) for ip=1:np]
	return vec(vcat(vcat(vect...)...))
end

"""
Data type for the source related parameters during acquisiton.

# Fields
* `nss::Int64` : number of supersources
* `ns::Array{Int64}` : number of sources for each supersource
* `fields::Vector{Symbol}` : number of fields
* `wav::Array{Float64}` : wavelets in time domain
* `tgrid` : time grid
"""
mutable struct Src
	nss::Int64
	ns::Vector{Int64}
	fields::Vector{Symbol}
	wav::Array{Array{Float64,2},2}
	tgrid::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
	Src(nss, ns, fields, wav, tgrid) =
		any([
		  any([fields[ifield] ∉ [:P, :Vx, :Vz] for ifield in 1:length(fields)]),
		  length(fields) == 0,
		  broadcast(size,wav) ≠ [(length(tgrid),ns[iss]) for iss=1:nss, ifield=1:length(fields)]
		  ]) ?
		error("error in Src construction") : new(nss, ns, fields, wav, tgrid)
end

function Base.isapprox(src1::Src, src2::Src)
	vec=([(src1.nss==src2.nss), (src1.fields==src2.fields), (src1.ns==src2.ns),
       		(isequal(src1.tgrid, src2.tgrid)),
		])
	return all(vec)
end

function Base.copyto!(srco::Src, src::Src)
	if(isapprox(srco, src))
		for f in [:wav]
			srcowav=getfield(srco, f)
			srcwav=getfield(src, f)
			for iss=1:src.nss, ifield=1:length(src.fields)
				srcowavv=srcowav[iss,ifield]
				srcwavv=srcwav[iss,ifield]
				copyto!(srcowavv, srcwavv)
			end
		end
		return srco
	else
		error("attempt to copy dissimilar sources")
	end
end

"""
Allocate `Src` with zeros depending on the acquisition geometry.
"""
function Src_zeros(acqgeom::Geom,  fields::Vector{Symbol}, tgrid::StepRangeLen)
	wavsrc = [zeros(length(tgrid),acqgeom.ns[iss]) for iss=1:acqgeom.nss, ifield=1:length(fields)]
	return Src(acqgeom.nss, acqgeom.ns, fields, wavsrc, deepcopy(tgrid))
end

"""
Priiiint information about `Src`
"""
function Base.print(src::Src, name::String="")
	println("\tSource Acquisition:\t",name)
	println("\t> number of supersources:\t",src.nss)
	println("\t> sources per supersource:\t","min\t",minimum(src.ns[:]), "\tmax\t", maximum(src.ns[:]))

	freqmin, freqmax, freqpeak = freqs(src)
	println("\t> frequency:\t","min\t",freqmin, "\tmax\t",freqmax,"\tpeak\t",freqpeak)
	println("\t> time:\t","min\t",src.tgrid[1], "\tmax\t", src.tgrid[end])
	println("\t> samples:\t",length(src.tgrid))
end

"""
Return minimum, maximum and peak frequencies of `Src`
"""
function freqs(src::Src)
	freqmin = minimum([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:min) for i in 1:src.nss, j in 1:length(src.fields)])
	freqmax = maximum([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:max) for i in 1:src.nss, j in 1:length(src.fields)])
	freqpeak = mean([Utils.findfreq(src.wav[i,j][:,:],src.tgrid,attrib=:peak) for i in 1:src.nss, j in 1:length(src.fields)])
	return freqmin, freqmax, freqpeak
end

"""
Constructor for `Src` data type.
Uses same source wavelet, i.e., `wav` for all sources and supersources

# Arguments

* `nss::Int64` : number of supersources
* `ns::Int64` : number of sources
* `fields::Vector{Symbol}` : number of fields the sources are exciting
* `wav::Array{Float64}` : a source wavelet that is used for all sources and supersources
* `tgrid` : time grid for the wavelet
"""
function Src_fixed(nss::Int64,
	     ns::Int64,
	     fields::Vector{Symbol},
	     wav::Array{Float64},
	     tgrid::StepRangeLen
	     )

	wavsrc = [repeat(wav,inner=(1,ns)) for iss=1:nss, ifield=1:length(fields)]
	return Src(nss, fill(ns, nss), fields, wavsrc, deepcopy(tgrid))
end

"""
Constructor of `Src`, which is typical for a input model such that
the model has `nλ` wavelengths.

# Arguments

* `nss::Int64` : number of supersources
* `ns::Int64` : number of source per each supersource
* `fields::Vector{Symbol}` :

# Keyword Arguments

* `mod::Models.Seismic` :
* `nλ::Int64=10` : number of wavelengths in the mod
* `wav_func::Function=(fqdom, tgrid)->Utils.Wavelets.ricker(fqdom,tgrid)` : which wavelet to generate, see Utils.Wavelets.jl
* `tmaxfrac::Float64=1.0` : by default the maximum modelling time is computed using the average velocity and the diagonal distance of the model, use this fraction to increase of reduce the maximum time
"""
function Src_fixed_mod(
		nss::Int64, ns::Int64, fields::Vector{Symbol};
		mod::Models.Seismic=nothing,
		nλ::Int64=10,
		wav_func::Function=(fqdom, tgrid)->Utils.Wavelets.ricker(fqdom,tgrid),
		tmaxfrac::Float64=1.0
		)

	x=mod.mgrid[2]; z=mod.mgrid[1]
	# dominant wavelength using mod dimensions
	λdom=mean([(abs(x[end]-x[1])), (abs(z[end]-z[1]))])/real(nλ)
	# average P velocity
	vavg=Models.χ([mean(mod.χvp)], mod.ref.vp, -1)[1]

	fqdom = vavg/λdom

	# maximum distance the wave travels
	d = sqrt((x[1]-x[end]).^2+(z[1]-z[end]).^2)

	# use two-way maximum distance to get tmax
	tmax=2.0*d*inv(vavg)*tmaxfrac

	# choose sampling interval to obey max freq of source wavelet
	δmin = minimum([step(mod.mgrid[2]), step(mod.mgrid[1])])
	vpmax = mod.vp0[2]
	δt=0.5*δmin/vpmax

	# check if δt is reasonable
	#(δt > 0.1/fqdom) : error("decrease spatial sampling or nλ")

	tgrid=range(0.0, stop=tmax, step=δt)
	wav = wav_func(fqdom, tgrid);

	wavsrc = [hcat([wav_func(fqdom, tgrid) for is in 1:ns]...) for iss=1:nss, ifield=1:length(fields)]
	src=Src(nss, fill(ns, nss), fields, wavsrc, deepcopy(tgrid))
	print(src)
	# choose ricker waveletes of fdom
	return src
end

"""
Generate band-limited random source signals
"""
function Src_fixed_random(nss::Int64, ns::Int64, fields::Vector{Symbol};
			  distvec=[Normal() for iss in 1:nss],
			  sparsepvec=[1. for iss in 1:nss],
			  fmin::Float64=0.0,
			  fmax::Float64=0.0,
			  tgrid::StepRangeLen=nothing,
			  tmaxfrac::Float64=1.0)
	wavsrc = [repeat(zeros(length(tgrid)),inner=(1,ns)) for iss=1:nss, ifield=1:length(fields)]
	for ifield in 1:length(fields), iss in 1:nss, is in 1:ns
		wavsrc[iss, ifield][:,is] = Utils.get_tapered_random_tmax_signal(tgrid, fmin=fmin, fmax=fmax, tmaxfrac=tmaxfrac, dist=distvec[iss],
								 sparsep=sparsepvec[iss])
	end
	src=Src(nss, fill(ns, nss), fields, wavsrc, deepcopy(tgrid))
	print(src)
	return src
end

"""
Function that returns Src after time reversal
"""
function Src_tr(src::Src)
	return Src(src.nss,src.ns,src.fields,
	    [reverse(src.wav[i,j],dims=1) for i in 1:src.nss, j in 1:length(src.fields)],deepcopy(src.tgrid))
end

"""
Pad `Src`
tgrids should be same in all Src
"""
function Src_uspos(src::Vector{Src}, acq::Vector{Geom})
	np = length(src) == length(acq) ? length(src) : error("unequal sizez")

	# unique source positions
	nus=Geom_get(acq,:nus)
	uspos=Geom_get(acq,:uspos)
	# all zeros for all unique positions
	wavout = [[zeros(src[ip].length(tgrid),nus) for iss=1:src[ip].nss, ifield=1:length(src[ip].fields)] for ip=1:np]
	# fill source wavelets when necessary
	for ip=1:np
		for ifield=1:length(src[ip].fields), iss=1:src[ip].nss, is=1:acq[ip].ns[iss]
			is0=findall([[uspos[1][i]-acq[ip].sz[iss][is],
		      uspos[2][i]-acq[ip].sx[iss][is]] == [0., 0.,] for i in 1:nus])

			wavout[ip][iss, ifield][:,is0] = src[ip].wav[iss, ifield][:,is]
		end
	end
	# output src
	return [Src(src[ip].nss,fill(nus, src[ip].nss),
	     src[ip].fields,wavout[ip],src[ip].tgrid) for ip=1:np]
end

"""
Return a sparse ACQ matrix, for a given acqgeom
data=ACQ*snapshot
"""
function ACQmat(acqgeom::Geom,mgrid,iss=1)
	nz,nx=length.(mgrid)
	ACQ=spzeros(acqgeom.nr[iss],prod(length.(mgrid)))
	for ir = 1:acqgeom.nr[iss]
		irx=argmin(abs.(mgrid[2].-acqgeom.rx[iss][ir]))
		irz=argmin(abs.(mgrid[1].-acqgeom.rz[iss][ir]))
		ACQ[ir,irz+(irx-1)*nz]=1.0
	end
	return ACQ
end

"""
return a vector of the order

"""
function Src_getvec(src::Vector{Src}, field::Symbol)
	np = length(src);
	vect = [getfield(src[ip],field) for ip=1:np]
	return vec(hcat(hcat(vect...)...))
end

function save(geom, folder)
	!(isdir(folder)) && error("invalid directory")
	for iss in 1:geom.nss
		file=joinpath(folder, string("s", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(geom.sx[iss], geom.sz[iss])))
		file=joinpath(folder, string("r", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(geom.rx[iss], geom.rz[iss])))
	end
end

end # module


text:
## PART 5

println("\nPart 5");
println("\n\tCalculating σ̂ for each a");

# Calculate squared standard deviations
σsq = Float64[];
σ = Float64[];
for a in [1:5]
  σsqa = dot( v[:,a] .- v̄̂[a] , v[:,a] .- v̄̂[a] ) / (M-1);
  σa = sqrt(σsqa);
  println("\t\tσ̂[$a] = $σa");
  push!(σ, σa);
  push!(σsq, σsqa);
end

println("\n\tRecall, from Part 2:");
println("\n\tN1 = $N1\n\tN2 = $N2\n");
println("\n\tWe expect:");
println("\n\t\t2 * τ̂[a] * ( σ̂[a]^2 / σ̂1[a]^2 ) = N1 = $N1")
println("\n\tand")
println("\n\t\t2 * τ̂[a] * ( σ̂[a]^2 / σ̂2[a]^2 ) = N2 = $N2")

println("\n\tTest that hypothesis...");
println("\n\tfor N1:");

for a in [1:5]
  expectn = 2 * tau[a] * σsq[a] / σ̂1[a]^2;
  println("\n\t\t2 * τ̂[a] * ( σ̂[a]^2 / σ̂1[a]^2 ) = N1 = $N1");
end

println("\n\t...and for N2:");
for a in [1:5]
  expectn = 2 * tau[a] * σsq[a] / σ̂2[a]^2;
  println("\n\t\t2 * τ̂[a] * (  σ̂[a]^2 / σ̂2[a]^2 ) = N2 = $N2");
end

println("\n\tWhich is exactly what we would expect.");


text:
export ConstVariable, constvar, getconst, isconnected

import Rocket: SingleObservable, AsapScheduler
import Base: getindex, show

mutable struct ConstVariable{C, M} <: AbstractVariable
    name            :: Symbol
    collection_type :: AbstractVariableCollectionType
    constant        :: C
    messageout      :: M
    nconnected      :: Int
end

Base.show(io::IO, constvar::ConstVariable) = print(io, "ConstVariable(", indexed_name(constvar), ")")

constvar(name::Symbol, constval, collection_type::AbstractVariableCollectionType = VariableIndividual())                 = ConstVariable(name, collection_type, constval, of(Message(constval, true, false)), 0)
constvar(name::Symbol, constval::Real, collection_type::AbstractVariableCollectionType = VariableIndividual())           = constvar(name, PointMass(constval), collection_type)
constvar(name::Symbol, constval::AbstractVector, collection_type::AbstractVariableCollectionType = VariableIndividual()) = constvar(name, PointMass(constval), collection_type)
constvar(name::Symbol, constval::AbstractMatrix, collection_type::AbstractVariableCollectionType = VariableIndividual()) = constvar(name, PointMass(constval), collection_type)

function constvar(name::Symbol, fn::Function, dims::Tuple)
    return constvar(name, fn, dims...)
end

function constvar(name::Symbol, fn::Function, length::Int)
    vars = Vector{ConstVariable}(undef, length)
    @inbounds for i in 1:length
        vars[i] = constvar(name, fn(i), VariableVector(i))
    end
    return vars
end

function constvar(name::Symbol, fn::Function, dims::Vararg{Int})
    vars = Array{ConstVariable}(undef, dims)
    @inbounds for i in CartesianIndices(axes(vars))
        vars[i] = constvar(name, fn(convert(Tuple, i)), VariableArray(i))
    end
    return vars
end

degree(constvar::ConstVariable)          = nconnected(constvar)
name(constvar::ConstVariable)            = constvar.name
collection_type(constvar::ConstVariable) = constvar.collection_type

Base.getindex(constvar::ConstVariable, index) = Base.getindex(getconstant(constvar), index)

isconnected(constvar::ConstVariable) = constvar.nconnected !== 0
nconnected(constvar::ConstVariable)  = constvar.nconnected

getconst(constvar::ConstVariable{ <: PointMass }) = getpointmass(constvar.constant)
getconst(constvar::ConstVariable)                 = constvar.constant

getlastindex(::ConstVariable) = 1

messageout(constvar::ConstVariable, ::Int) = constvar.messageout
messagein(constvar::ConstVariable, ::Int)  = error("It is not possible to get a reference for inbound message for constvar")

get_pipeline_stages(::ConstVariable) = EmptyPipelineStage()

_getmarginal(constvar::ConstVariable)      = of(Marginal(constvar.constant, true, false))
_setmarginal!(::ConstVariable, observable) = error("It is not possible to set a marginal stream for `ConstVariable`")
_makemarginal(::ConstVariable)             = error("It is not possible to make marginal stream for `ConstVariable`")

# For _getmarginal
function Rocket.getrecent(observable::SingleObservable{ <: Marginal })
    return observable.value
end

function setmessagein!(constvar::ConstVariable, ::Int, messagein)
    constvar.nconnected += 1
    return nothing
end


text:
using LinearAlgebra
using FermiCG
using Printf
using Test
using Arpack

@testset "fci" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,0,1]))
    #push!(atoms,Atom(3,"H",[0,0,2]))
    #push!(atoms,Atom(4,"H",[0,0,3]))
    #basis = "6-31g"
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    n_elec_a = 1
    n_elec_b = 1

    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,1,1)
    @printf(" FCI Energy: %12.8f\n", e_fci)

    norb = size(ints.h1)[1]
    problem = FermiCG.StringCI.FCIProblem(norb, n_elec_a, n_elec_b)

    display(problem)

    @time Hmat = FermiCG.StringCI.build_H_matrix(ints, problem)
    Ssq = FermiCG.StringCI.build_S2_matrix(problem)
    println("Ssq")
    println(Ssq)
    F = eigen(Hmat+0.01*Ssq)
    e = F.values
    v = F.vectors
    e = diag(v'*Hmat*v)
    ss = diag(v'*Ssq*v)

    e = real(e)
    for (i,ei) in enumerate(e)
        @printf(" Energy: %12.8f   <S2>: %12.8f\n",ei+ints.h0,ss[i])
    end
    @test isapprox(e[1], e_fci , atol=1e-10)
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
module PolyJuMP

# MultivariatePolynomials extension

using MultivariatePolynomials
import MultivariateBases
const MB = MultivariateBases
using MultivariateMoments
using SemialgebraicSets

# MOI extension

using MathOptInterface
const MOI = MathOptInterface

include("attributes.jl")
include("zero_polynomial.jl")
include("functions.jl")

# Bridges
const MOIB = MOI.Bridges
include("zero_polynomial_bridge.jl")
include("zero_polynomial_in_algebraic_set_bridge.jl")
include("plus_minus_bridge.jl")

# JuMP extension

using JuMP

include("variable.jl")
include("constraint.jl")

include("data.jl")
include("default.jl")

end # module


text:
@polly function kernel_adi(tsteps, u, v, p, q)
    n = size(u,1)

    DX = 1.0 / n
    DY = 1.0 / n
    DT = 1.0 / tsteps
    B1 = 2.0
    B2 = 1.0
    mul1 = B1 * DT / (DX * DX)
    mul2 = B2 * DT / (DY * DY)

    a = -mul1 /  2.0;
    b = 1.0 + mul1
    c = a
    d = -mul2 / 2.0
    _e = 1.0 + mul2
    f = d

    for t = 1:tsteps
        # Column Sweep
        for i = 2:(n-1)
            v[1,i] = 1.0
            p[i,1] = 0.0
            q[i,1] = v[1,i]
            for j = 2:(n-1)
                p[i,j] = -c / (a * p[i,j-1] + b)
                q[i,j] = (-d*u[j,i-1]+(1.0+2.0*d)*u[j,i] - f*u[j,i+1]-a*q[i,j-1])/(a*p[i,j-1]+b)
            end

            v[n,i] = 1.0
            for j = (n-1):-1:2
                v[j,i] = p[i,j] * v[j+1,i] + q[i,j]
            end
        end

        # Row Sweep
        for i = 2:(n-1)
            u[i,1] = 1.0
            p[i,1] = 0.0
            q[i,1] = u[i,1]
            for j = 2:(n-1)
                p[i,j] = -f / (d*p[i,j-1] + _e)
                q[i,j] = (-a*v[i-1,j]+(1.0+2.0*a)*v[i,j] - c*v[i+1,j]-d*q[i,j-1])/(d*p[i,j-1]+_e)
            end

            u[i,n] = 1.0
            for j = (n-1):-1:2
                u[i,j] = p[i,j] * u[i,j+1] + q[i,j]
            end
        end
    end
end

let
    tsteps = 500
    n = 1000

    u = zeros(Float32, n, n)
    v = zeros(Float32, n, n)
    p = zeros(Float32, n, n)
    q = zeros(Float32, n, n)

    for i = 1:n, j = 1:n
        u[i,j] =  (i + n-j) / n
    end

    SUITE["adi"] = @benchmarkable kernel_adi($tsteps, u, v, p, q) setup = (u = copy($u); v = copy($v); p = copy($p); q = copy($q))
end


text:
module TopologicalDataAnalysis

# Write your package code here.
using LinearAlgebra
using Combinatorics

export getpowerset

"""
    getpowerset(L)

Returns the power set of the list L.

Input:
    L (list).

Output:
    L2 (list): power set of L.

Example:
    L = collect(1:3)
    getpowerset(L)
    [Int64[], [1], [2], [3], [1, 2], [1, 3], [2, 3], [1, 2, 3]]
"""
function getpowerset(L)
    return collect(powerset(L))
end

end


text:
module Simpson
    using Base.Threads

    function simpson(func, x1::BigFloat, x2::BigFloat, n)
        dh = (x2 - x1) / BigFloat(n)

        # Simpsonの公式によって数値積分する
        m = div(n, 2)
        results = Array{BigFloat}(undef, m)
        @inbounds @threads for i = 1:m
            fx1 = func(x1 + dh * BigFloat(i * 2 - 1))
            fx2 = func(x1 + dh * BigFloat(i * 2))

            results[i] = BigFloat(4) * fx1 + BigFloat(2) * fx2
        end
        sum = func(x1) + reduce(+, results) - func(x2)

        return dh / BigFloat(3) * sum
    end
end


text:
using DelimitedFiles

function support_serial_qtest(mdp::Union{MDP,POMDP}, file::AbstractString; niter::Int64=100, res::Float64=1e-3)
    qt = readdlm(file)
    solver = ValueIterationSolver(max_iterations=niter, belres=res, verbose=true)
    policy = solve(solver, mdp)
    (q, u, p, am) = locals(policy)
    npolicy = ValueIterationPolicy(mdp, deepcopy(q))
    nnpolicy = ValueIterationPolicy(mdp, deepcopy(q), deepcopy(u), deepcopy(p))
    s = GridWorldState(1,1)
    a1 = action(policy, s)
    @inferred action(policy, s)
    v1 = value(policy, s)
    a2 = action(npolicy, s)
    v2 = value(npolicy, s)
    return (isapprox(qt, q, rtol=1e-5)) && (policy.policy == nnpolicy.policy)
end

function test_complex_gridworld()
    # Load correct policy from file and verify we can reconstruct it
    rstates = [GridWorldState(4,3), GridWorldState(4,6), GridWorldState(9,3), GridWorldState(8,8)]
    rvals = [-10.0, -5.0, 10.0, 3.0]
    xs = 10
    ys = 10
    mdp = LegacyGridWorld(sx=xs, sy=ys, rs = rstates, rv = rvals)
    file = "grid-world-10x10-Q-matrix.txt"
    niter = 100
    res = 1e-3

    return support_serial_qtest(mdp, file, niter=niter, res=res)
end

function test_simple_grid()
    # Simple test....
    # GridWorld(sx=2,sy=3) w reward at (2,3):
    # Here's our grid:
    # |state (x,y)____available actions__|
    # ----------------------------------------------
    # |5 (1,3)__u,d,l,r__|6 (2,3)__u,d,l,r+REWARD__|
    # |3 (1,2)__u,d,l,r__|4 (2,2)______u,d,l,r_____|
    # |1 (1,1)__u,d,l,r__|2 (2,1)______u,d,l,r_____|
    # ----------------------------------------------
    # 7 (0,0) is absorbing state
    mdp = LegacyGridWorld(sx=2, sy=3, rs = [GridWorldState(2,3)], rv = [10.0])

    solver = ValueIterationSolver(verbose=true)
    policy = solve(solver, mdp)

    # up: 1, down: 2, left: 3, right: 4
    correct_policy = [1,1,1,1,4,1,1] # alternative policies
    # are possible, but since they are tied & the first
    # action is always 1, we will always return 1 for tied
    # actions
    return policy.policy == correct_policy
end

function test_init_solution()
    # Initialize the value to the solution
    rstates = [GridWorldState(4,3), GridWorldState(4,6), GridWorldState(9,3), GridWorldState(8,8)]
    rvals = [-10.0, -5.0, 10.0, 3.0]
    xs = 10
    ys = 10
    mdp = LegacyGridWorld(sx=xs, sy=ys, rs = rstates, rv = rvals)
    qt = readdlm("grid-world-10x10-Q-matrix.txt")
    ut = maximum(qt, dims=2)[:]
    solver = ValueIterationSolver(verbose=true, init_util=ut, belres=1e-3)
    policy = solve(solver, mdp)
    return isapprox(ut, policy.util, rtol=1e-5)
end

function test_not_include_Q()
    # Load correct policy from file and verify we can reconstruct it
    rstates = [GridWorldState(4,3), GridWorldState(4,6), GridWorldState(9,3), GridWorldState(8,8)]
    rvals = [-10.0, -5.0, 10.0, 3.0]
    xs = 10
    ys = 10
    mdp = LegacyGridWorld(sx=xs, sy=ys, rs = rstates, rv = rvals)
    qt = readdlm("grid-world-10x10-Q-matrix.txt")
    ut = maximum(qt, dims=2)[:]
    niter = 100
    res = 1e-3
    solver = ValueIterationSolver(verbose=true, init_util=ut, belres=1e-3, include_Q=false)
    policy = solve(solver, mdp)
    return isapprox(ut, policy.util, rtol=1e-3)
end

function test_warning()
    mdp = LegacyGridWorld()
    solver = ValueIterationSolver()
    println("There should be a warning bellow: ")
    solve(solver, mdp, verbose=true)
end

function test_type_stability()
    mdp = LegacyGridWorld()
    solver = ValueIterationSolver()
    @inferred solve(solver, mdp)
    return true
end

@test test_complex_gridworld() == true
@test test_simple_grid() == true
@test test_init_solution() == true
@test test_not_include_Q() == true
test_warning()
@test_throws String solve(ValueIterationSolver(), TigerPOMDP())
@test test_type_stability() == true


text:
Network = NetLib.SimpleNet

netparams = NetLib.SimpleNetHP(
  width=200,
  depth_common=6,
  use_batch_norm=true,
  batch_norm_momentum=1.)

self_play = SelfPlayParams(
  sim=SimParams(
    num_games=1000,
    num_workers=128,
    batch_size=128,
    use_gpu=false,
    reset_every=4,
    flip_probability=0.,
    alternate_colors=false),
  mcts=MctsParams(
    num_iters_per_turn=400,
    cpuct=1.0,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.2,
    dirichlet_noise_α=1.0))

arena = ArenaParams(
  sim=SimParams(
    num_games=100,
    num_workers=100,
    batch_size=100,
    use_gpu=false,
    reset_every=1,
    flip_probability=0.,
    alternate_colors=true),
  mcts = MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.3),
    dirichlet_noise_ϵ=0.1),
  update_threshold=0.00)

learning = LearningParams(
  use_gpu=false,
  samples_weighing_policy=LOG_WEIGHT,
  l2_regularization=1e-4,
  optimiser=CyclicNesterov(
    lr_base=1e-3,
    lr_high=1e-2,
    lr_low=1e-3,
    momentum_high=0.9,
    momentum_low=0.8),
  batch_size=32,
  loss_computation_batch_size=2048,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=0,
  max_batches_per_checkpoint=5_000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=2,
  memory_analysis=MemAnalysisParams(
    num_game_stages=4),
  ternary_rewards=true,
  use_symmetries=true,
  mem_buffer_size=PLSchedule(80_000))

benchmark_sim = SimParams(
  arena.sim;
  num_games=400,
  num_workers=100,
  batch_size=100)

benchmark = [
  Benchmark.Duel(
    Benchmark.Full(self_play.mcts),
    Benchmark.MctsRollouts(self_play.mcts),
    benchmark_sim)
  ]

experiment = Experiment(
  "ospiel_ttt", GameSpec(), params, Network, netparams, benchmark)


text:
abstract type Product end
abstract type OperatingProduct <: Product end
abstract type ReserveDirection end
abstract type ReserveUpEMIS <: ReserveDirection end
abstract type ReserveDownEMIS <: ReserveDirection end
abstract type ReserveEMIS{T<:ReserveDirection} <: OperatingProduct end

get_name(prod::Product) = prod.name
get_marginal_cost(prod::OperatingProduct) = prod.marginal_cost


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: codestar

using AWS.Compat
using AWS.UUIDs
"""
    AssociateTeamMember()

Adds an IAM user to the team for an AWS CodeStar project.

# Required Parameters
- `projectId`: The ID of the project to which you will add the IAM user.
- `projectRole`: The AWS CodeStar project role that will apply to this user. This role determines what actions a user can take in an AWS CodeStar project.
- `userArn`: The Amazon Resource Name (ARN) for the IAM user you want to add to the AWS CodeStar project.

# Optional Parameters
- `clientRequestToken`: A user- or system-generated token that identifies the entity that requested the team member association to the project. This token can be used to repeat the request.
- `remoteAccessAllowed`: Whether the team member is allowed to use an SSH public/private key pair to remotely access project resources, for example Amazon EC2 instances.
"""

associate_team_member(projectId, projectRole, userArn; aws_config::AWSConfig=global_aws_config()) = codestar("AssociateTeamMember", Dict{String, Any}("projectId"=>projectId, "projectRole"=>projectRole, "userArn"=>userArn); aws_config=aws_config)
associate_team_member(projectId, projectRole, userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("AssociateTeamMember", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId, "projectRole"=>projectRole, "userArn"=>userArn), args)); aws_config=aws_config)

"""
    CreateProject()

Creates a project, including project resources. This action creates a project based on a submitted project request. A set of source code files and a toolchain template file can be included with the project request. If these are not provided, an empty project is created.

# Required Parameters
- `id`: The ID of the project to be created in AWS CodeStar.
- `name`: The display name for the project to be created in AWS CodeStar.

# Optional Parameters
- `clientRequestToken`: A user- or system-generated token that identifies the entity that requested project creation. This token can be used to repeat the request.
- `description`: The description of the project, if any.
- `sourceCode`: A list of the Code objects submitted with the project request. If this parameter is specified, the request must also include the toolchain parameter.
- `tags`: The tags created for the project.
- `toolchain`: The name of the toolchain template file submitted with the project request. If this parameter is specified, the request must also include the sourceCode parameter.
"""

create_project(id, name; aws_config::AWSConfig=global_aws_config()) = codestar("CreateProject", Dict{String, Any}("id"=>id, "name"=>name); aws_config=aws_config)
create_project(id, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("CreateProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "name"=>name), args)); aws_config=aws_config)

"""
    CreateUserProfile()

Creates a profile for a user that includes user preferences, such as the display name and email address assocciated with the user, in AWS CodeStar. The user profile is not project-specific. Information in the user profile is displayed wherever the user's information appears to other users in AWS CodeStar.

# Required Parameters
- `displayName`: The name that will be displayed as the friendly name for the user in AWS CodeStar.
- `emailAddress`: The email address that will be displayed as part of the user's profile in AWS CodeStar.
- `userArn`: The Amazon Resource Name (ARN) of the user in IAM.

# Optional Parameters
- `sshPublicKey`: The SSH public key associated with the user in AWS CodeStar. If a project owner allows the user remote access to project resources, this public key will be used along with the user's private key for SSH access.
"""

create_user_profile(displayName, emailAddress, userArn; aws_config::AWSConfig=global_aws_config()) = codestar("CreateUserProfile", Dict{String, Any}("displayName"=>displayName, "emailAddress"=>emailAddress, "userArn"=>userArn); aws_config=aws_config)
create_user_profile(displayName, emailAddress, userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("CreateUserProfile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("displayName"=>displayName, "emailAddress"=>emailAddress, "userArn"=>userArn), args)); aws_config=aws_config)

"""
    DeleteProject()

Deletes a project, including project resources. Does not delete users associated with the project, but does delete the IAM roles that allowed access to the project.

# Required Parameters
- `id`: The ID of the project to be deleted in AWS CodeStar.

# Optional Parameters
- `clientRequestToken`: A user- or system-generated token that identifies the entity that requested project deletion. This token can be used to repeat the request.
- `deleteStack`: Whether to send a delete request for the primary stack in AWS CloudFormation originally used to generate the project and its resources. This option will delete all AWS resources for the project (except for any buckets in Amazon S3) as well as deleting the project itself. Recommended for most use cases.
"""

delete_project(id; aws_config::AWSConfig=global_aws_config()) = codestar("DeleteProject", Dict{String, Any}("id"=>id); aws_config=aws_config)
delete_project(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("DeleteProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    DeleteUserProfile()

Deletes a user profile in AWS CodeStar, including all personal preference data associated with that profile, such as display name and email address. It does not delete the history of that user, for example the history of commits made by that user.

# Required Parameters
- `userArn`: The Amazon Resource Name (ARN) of the user to delete from AWS CodeStar.

"""

delete_user_profile(userArn; aws_config::AWSConfig=global_aws_config()) = codestar("DeleteUserProfile", Dict{String, Any}("userArn"=>userArn); aws_config=aws_config)
delete_user_profile(userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("DeleteUserProfile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("userArn"=>userArn), args)); aws_config=aws_config)

"""
    DescribeProject()

Describes a project and its resources.

# Required Parameters
- `id`: The ID of the project.

"""

describe_project(id; aws_config::AWSConfig=global_aws_config()) = codestar("DescribeProject", Dict{String, Any}("id"=>id); aws_config=aws_config)
describe_project(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("DescribeProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    DescribeUserProfile()

Describes a user in AWS CodeStar and the user attributes across all projects.

# Required Parameters
- `userArn`: The Amazon Resource Name (ARN) of the user.

"""

describe_user_profile(userArn; aws_config::AWSConfig=global_aws_config()) = codestar("DescribeUserProfile", Dict{String, Any}("userArn"=>userArn); aws_config=aws_config)
describe_user_profile(userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("DescribeUserProfile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("userArn"=>userArn), args)); aws_config=aws_config)

"""
    DisassociateTeamMember()

Removes a user from a project. Removing a user from a project also removes the IAM policies from that user that allowed access to the project and its resources. Disassociating a team member does not remove that user's profile from AWS CodeStar. It does not remove the user from IAM.

# Required Parameters
- `projectId`: The ID of the AWS CodeStar project from which you want to remove a team member.
- `userArn`: The Amazon Resource Name (ARN) of the IAM user or group whom you want to remove from the project.

"""

disassociate_team_member(projectId, userArn; aws_config::AWSConfig=global_aws_config()) = codestar("DisassociateTeamMember", Dict{String, Any}("projectId"=>projectId, "userArn"=>userArn); aws_config=aws_config)
disassociate_team_member(projectId, userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("DisassociateTeamMember", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId, "userArn"=>userArn), args)); aws_config=aws_config)

"""
    ListProjects()

Lists all projects in AWS CodeStar associated with your AWS account.

# Optional Parameters
- `maxResults`: The maximum amount of data that can be contained in a single set of results.
- `nextToken`: The continuation token to be used to return the next set of results, if the results cannot be returned in one response.
"""

list_projects(; aws_config::AWSConfig=global_aws_config()) = codestar("ListProjects"; aws_config=aws_config)
list_projects(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("ListProjects", args; aws_config=aws_config)

"""
    ListResources()

Lists resources associated with a project in AWS CodeStar.

# Required Parameters
- `projectId`: The ID of the project.

# Optional Parameters
- `maxResults`: The maximum amount of data that can be contained in a single set of results.
- `nextToken`: The continuation token for the next set of results, if the results cannot be returned in one response.
"""

list_resources(projectId; aws_config::AWSConfig=global_aws_config()) = codestar("ListResources", Dict{String, Any}("projectId"=>projectId); aws_config=aws_config)
list_resources(projectId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("ListResources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId), args)); aws_config=aws_config)

"""
    ListTagsForProject()

Gets the tags for a project.

# Required Parameters
- `id`: The ID of the project to get tags for.

# Optional Parameters
- `maxResults`: Reserved for future use.
- `nextToken`: Reserved for future use.
"""

list_tags_for_project(id; aws_config::AWSConfig=global_aws_config()) = codestar("ListTagsForProject", Dict{String, Any}("id"=>id); aws_config=aws_config)
list_tags_for_project(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("ListTagsForProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    ListTeamMembers()

Lists all team members associated with a project.

# Required Parameters
- `projectId`: The ID of the project for which you want to list team members.

# Optional Parameters
- `maxResults`: The maximum number of team members you want returned in a response.
- `nextToken`: The continuation token for the next set of results, if the results cannot be returned in one response.
"""

list_team_members(projectId; aws_config::AWSConfig=global_aws_config()) = codestar("ListTeamMembers", Dict{String, Any}("projectId"=>projectId); aws_config=aws_config)
list_team_members(projectId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("ListTeamMembers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId), args)); aws_config=aws_config)

"""
    ListUserProfiles()

Lists all the user profiles configured for your AWS account in AWS CodeStar.

# Optional Parameters
- `maxResults`: The maximum number of results to return in a response.
- `nextToken`: The continuation token for the next set of results, if the results cannot be returned in one response.
"""

list_user_profiles(; aws_config::AWSConfig=global_aws_config()) = codestar("ListUserProfiles"; aws_config=aws_config)
list_user_profiles(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("ListUserProfiles", args; aws_config=aws_config)

"""
    TagProject()

Adds tags to a project.

# Required Parameters
- `id`: The ID of the project you want to add a tag to.
- `tags`: The tags you want to add to the project.

"""

tag_project(id, tags; aws_config::AWSConfig=global_aws_config()) = codestar("TagProject", Dict{String, Any}("id"=>id, "tags"=>tags); aws_config=aws_config)
tag_project(id, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("TagProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "tags"=>tags), args)); aws_config=aws_config)

"""
    UntagProject()

Removes tags from a project.

# Required Parameters
- `id`: The ID of the project to remove tags from.
- `tags`: The tags to remove from the project.

"""

untag_project(id, tags; aws_config::AWSConfig=global_aws_config()) = codestar("UntagProject", Dict{String, Any}("id"=>id, "tags"=>tags); aws_config=aws_config)
untag_project(id, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("UntagProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "tags"=>tags), args)); aws_config=aws_config)

"""
    UpdateProject()

Updates a project in AWS CodeStar.

# Required Parameters
- `id`: The ID of the project you want to update.

# Optional Parameters
- `description`: The description of the project, if any.
- `name`: The name of the project you want to update.
"""

update_project(id; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateProject", Dict{String, Any}("id"=>id); aws_config=aws_config)
update_project(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateProject", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    UpdateTeamMember()

Updates a team member's attributes in an AWS CodeStar project. For example, you can change a team member's role in the project, or change whether they have remote access to project resources.

# Required Parameters
- `projectId`: The ID of the project.
- `userArn`: The Amazon Resource Name (ARN) of the user for whom you want to change team membership attributes.

# Optional Parameters
- `projectRole`: The role assigned to the user in the project. Project roles have different levels of access. For more information, see Working with Teams in the AWS CodeStar User Guide.
- `remoteAccessAllowed`: Whether a team member is allowed to remotely access project resources using the SSH public key associated with the user's profile. Even if this is set to True, the user must associate a public key with their profile before the user can access resources.
"""

update_team_member(projectId, userArn; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateTeamMember", Dict{String, Any}("projectId"=>projectId, "userArn"=>userArn); aws_config=aws_config)
update_team_member(projectId, userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateTeamMember", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectId"=>projectId, "userArn"=>userArn), args)); aws_config=aws_config)

"""
    UpdateUserProfile()

Updates a user's profile in AWS CodeStar. The user profile is not project-specific. Information in the user profile is displayed wherever the user's information appears to other users in AWS CodeStar.

# Required Parameters
- `userArn`: The name that will be displayed as the friendly name for the user in AWS CodeStar.

# Optional Parameters
- `displayName`: The name that is displayed as the friendly name for the user in AWS CodeStar.
- `emailAddress`: The email address that is displayed as part of the user's profile in AWS CodeStar.
- `sshPublicKey`: The SSH public key associated with the user in AWS CodeStar. If a project owner allows the user remote access to project resources, this public key will be used along with the user's private key for SSH access.
"""

update_user_profile(userArn; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateUserProfile", Dict{String, Any}("userArn"=>userArn); aws_config=aws_config)
update_user_profile(userArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar("UpdateUserProfile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("userArn"=>userArn), args)); aws_config=aws_config)