text:
export reconstruction_simple, reconstruction_multiEcho, reconstruction_multiCoil, reconstruction_multiCoilMultiEcho, reconstruction_lowRank, RecoParameters

"""
Performs iterative image reconstruction independently for the data of all coils,
contrasts and slices

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_simple( acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = sparseTrafo

  # reconstruction
  Ireco = zeros(ComplexF64, prod(reconSize), numSl, numContr, numChan)
  @sync for k = 1:numSl
    Threads.@spawn begin
      F = encodingOps2d_simple(acqData, reconSize, slice=k; encParams...)
      for j = 1:numContr
        W = WeightingOp(weights[j])
        for i = 1:numChan
          kdata = kData(acqData,j,i,k).* weights[j]
          solver = createLinearSolver(solvername, W∘F[j]; reg=reg, params...)

          I = solve(solver, kdata, startVector=get(params,:startVector,ComplexF64[]),
                                 solverInfo=get(params,:solverInfo,nothing))

          if isCircular( trajectory(acqData, j) )
            circularShutter!(reshape(I, reconSize), 1.0)
          end
          Ireco[:,k,j,i] = I
        end
      end
    end
  end

  Ireco = reshape(Ireco, reconSize..., numSl, numContr, numChan)
  return makeAxisArray(Ireco, acqData)
end

"""
Performs a iterative image reconstruction jointly for all contrasts. Different slices and coil images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiEcho(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr)... )

  W = WeightingOp( vcat(weights...) )

  # reconstruction
  Ireco = zeros(ComplexF64, prod(reconSize)*numContr, numChan, numSl)
  @sync for i = 1:numSl
    Threads.@spawn begin
      F = encodingOp2d_multiEcho(acqData, reconSize, slice=i; encParams...)
      for j = 1:numChan
        kdata = multiEchoData(acqData, j, i) .* vcat(weights...)
        solver = createLinearSolver(solvername, W∘F; reg=reg, params...)

        Ireco[:,j,i] = solve(solver,kdata; params...)
        # TODO circular shutter
      end
    end
  end

  Ireco = reshape(Ireco, reconSize..., numContr, numChan, numSl)
  return makeAxisArray(permutedims(Ireco,[1,2,5,3,4]), acqData)
end

"""
Performs a SENSE-type iterative image reconstruction. Different slices and contrasts images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoil(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , senseMaps::Array{ComplexF64}
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = sparseTrafo

  # solve optimization problem
  Ireco = zeros(ComplexF64, prod(reconSize), numSl, numContr, 1)
  @sync for k = 1:numSl
    Threads.@spawn begin
      E = encodingOps2d_parallel(acqData, reconSize, senseMaps; slice=k, encParams...)

      for j = 1:numContr
        W = WeightingOp(weights[j],numChan)
        kdata = multiCoilData(acqData, j, k) .* repeat(weights[j], numChan)

        EFull = ∘(W, E[j], isWeighting=true)
        solver = createLinearSolver(solvername, EFull; reg=reg, params...)
        I = solve(solver, kdata; params...)

        if isCircular( trajectory(acqData, j) )
          circularShutter!(reshape(I, reconSize), 1.0)
        end
        Ireco[:,k,j] = I
      end
    end
  end

  Ireco = reshape(Ireco, reconSize..., numSl, numContr, 1)
  return makeAxisArray(Ireco, acqData)
end

"""
Performs a SENSE-type iterative image reconstruction which reconstructs all contrasts jointly.
Different slices are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Vector{Regularization}
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , senseMaps::Array{ComplexF64}
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  encParams = getEncodingOperatorParams(;params...)

  # set sparse trafo in reg
  reg[1].params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr)... )

  W = WeightingOp( vcat(weights...), numChan )

  Ireco = zeros(ComplexF64, prod(reconSize), numContr, numSl)
  @sync for i = 1:numSl
    Threads.@spawn begin
      E = encodingOp2d_multiEcho_parallel(acqData, reconSize, senseMaps; slice=i, encParams...)

      kdata = multiCoilMultiEchoData(acqData, i) .* repeat(vcat(weights...), numChan)
      solver = createLinearSolver(solvername, W∘E; reg=reg, params...)

      Ireco[:,:,i] = solve(solver, kdata; params...)
    end
  end
  return reshape( permutedims(Ireco, [1,3,2]), reconSize..., numSl, numContr )
end


text:
@testset "Born-Mayer Unit Tests" begin
    A = 1.0u"eV"
    ρ = 0.25u"bohr"
    σ = 0.25u"bohr"
    C = 1.0u"eV*Å"
    D = 1.0u"eV*Å"
    rcutoff = 2.0u"Å"
    species = [:Ar, :H]
    p = BornMayer(A, ρ, σ, C, D, rcutoff, species)

    @test p isa EmpiricalPotential{NamedTuple{(:A, :ρ, :σ, :C, :D)},NamedTuple{(:rcutoff,)}}

    @test get_rcutoff(p) == austrip(rcutoff)
    @test get_species(p) == (:Ar, :H)

    @test get_parameters(p) == (; A=austrip(A), ρ=austrip(ρ), σ = austrip(σ), C = austrip(C), D = austrip(D))
    @test set_parameters(p, (; A=2.0, ρ=3.0, σ = 4.0, C = 5.0, D = 6.0)) == BornMayer(2.0, 3.0, 4.0, 5.0, 6.0, austrip(rcutoff), species)
    @test serialize_parameters(p) == [austrip(A), austrip(ρ), austrip(σ), austrip(C), austrip(D)]
    @test deserialize_parameters(p, [2.0, 3.0, 4.0, 5.0, 6.0]) == BornMayer(2.0, 3.0, 4.0, 5.0, 6.0, austrip(rcutoff), species)

    @test get_hyperparameters(p) == (; rcutoff=austrip(rcutoff))
    @test set_hyperparameters(p, (; rcutoff=1.0)) == BornMayer(austrip(A), austrip(ρ), austrip(σ), austrip(C), austrip(D), 1.0, species)
    @test serialize_hyperparameters(p) == [austrip(rcutoff)]
    @test deserialize_hyperparameters(p, [1.0]) == BornMayer(austrip(A), austrip(ρ), austrip(σ), austrip(C), austrip(D), 1.0, species)

    r = @SVector[1.0, 1.0, 1.0]
    R = norm(r)

    @test potential_energy(R, p) isa Float64
    @test force(R, r, p) isa SVector{3,Float64}
end

ext: jl


text:
module GraphAlgorithms

using LightGraphs
using SimpleWeightedGraphs
using GLPK
using JuMP

export

# Operators
edges_of,unweighted_complement,

# Fractional numbers
fractional_covering_number, fractional_coloring_number,
fractional_matching_number,

# Sets
largest_clique_number, largest_cliques_set,
largest_independence_number, largest_independents_set,

# Diversity
diversity

GraphAlgorithms
include("diversity.jl")
include("fractional_numbers.jl")
include("operators.jl")
include("sets.jl")

end # module


text:
# Use
#
#     DOCUMENTER_DEBUG=true julia --color=yes make.jl local [fixdoctests]
#
# for local builds.

using Documenter
using LegendHDF5IO

makedocs(
    sitename = "LegendHDF5IO",
    modules = [LegendHDF5IO],
    format = :html,
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
        "LICENSE" => "LICENSE.md",
    ],
    doctest = ("fixdoctests" in ARGS) ? :fix : true,
    html_prettyurls = !("local" in ARGS),
    html_canonical = "https://legend-exp.github.io/LegendHDF5IO.jl/stable/",
)

deploydocs(
    repo = "github.com/legend-exp/LegendHDF5IO.jl.git"
)


text:
export marginalrule

@marginalrule Bernoulli(:out_p) (m_out::PointMass, m_p::Beta) = begin
    r = mean(m_out)
    p = prod(ProdAnalytical(), Beta(one(r) + r, 2one(r) - r), m_p)
    return (out = m_out, p = p)
end


text:
# diric.jl Dirichlet kernel tests

@testset "diric" begin
    @test_throws ArgumentError diric(0, -2)
    @test @inferred(diric(0, 4)) ≈ 1
    @test @inferred(diric(0 // 1, 5)) == 1
    @test @inferred(diric(4π, 4)) ≈ 1
    @test @inferred(diric(2π, 4)) ≈ -1
    @test @inferred(diric(0, 5)) ≈ 1
    @test @inferred(diric(2π, 5)) ≈ 1
    @test @inferred(diric(π, 5)) ≈ 1 // 5
    @test @inferred(diric(π/2, 5)) ≈ diric(-π/2, 5) ≈ -0.2
    @test abs(@inferred(diric(π/2, 4))) < eps()
    @test abs(@inferred(diric(3π/2, 4))) < eps()

    # check type inference
    @inferred diric(2, 4)
    @inferred diric(Float16(1), 4)
    @inferred diric(0.1f0, 4)
    @inferred diric(0., 4)
    @inferred diric(BigFloat(1), 4)

    # accuracy test near 2π for even n
    dirics(Ω,n) = real(exp(-1im * (n-1) * Ω/2)/n * sum(exp.(1im*(0:(n-1))*Ω)))
    Ω = 2π .+ 4π * (-101:2:101)/100 * 1e-9
    n = 400
    @test dirics.(Ω,n) ≈ diric.(Ω,n)
end


text:
using Robotlib
import Robotlib: Rt2T, T2R, T2t
include(joinpath(dirname(@__FILE__),"..","src","calibNAXP.jl"))

type Plane
    N
    N1
    d
    Pn
    Qn
end

function generateRandomPlane(no)
    if no == 1
        N_RB = [1, 0, 0]
    elseif no == 2
        N_RB = [0, 1, 0]
    elseif no == 3
        N_RB = [0, 0, 1]
    else
        N_RB = normalize(randn(3)) # Normal of calibration plane, given in RB
    end
    d_RB = 1 #abs(0.001*randn);
    Pn = N_RB*N_RB' # Projection matrix onto normal
    Qn = eye(3)-Pn # Projection matrix onto plane

    Plane(N_RB, [N_RB;1], d_RB, Pn, Qn)
end

function generateRandomPoses(N_poses, sigma_n,T_TF_S, plane)
    LN_S    = [0, 0, 1] # Normal of laser plane in sensor frame
    N_RB    = plane.N
    d_RB    = plane.d
    T_RB_TF = Array(Float64,4,4,N_poses)
    i       = 1
    lines_RB= zeros(3,N_poses)
    lines_S = zeros(3,N_poses)
    points_RB = zeros(3,N_poses)
    points_S = zeros(4,N_poses)
    while i <= N_poses
        T_RB_TF[:,:,i]      = Rt2T(rpy2R(1.5*randn(), 1.5*randn(), 1.5*randn()),0.6*randn(3) + 2*N_RB)
        T_RB_S              = T_RB_TF[:,:,i]*T_TF_S
        LN_RB               = T2R(T_RB_S)*LN_S
        intersection_line_RB= normalize(cross(N_RB,LN_RB))
        intersection_line_S = inv(T2R(T_RB_S))*intersection_line_RB
        lines_RB[:,i]       = intersection_line_RB
        lines_S[:,i]        = intersection_line_S
        laser_0_RB          = T_RB_S[1:3,4]
        laser_y_RB          = T_RB_S[1:3,2]
        lambda              = (d_RB-vecdot(N_RB,laser_0_RB))/(vecdot(N_RB,laser_y_RB))
        points_RB[:,i]      = laser_0_RB + lambda*laser_y_RB
        points_S[:,i]       = inv(T_RB_S)*[points_RB[:,i]; 1]
        if abs(lambda) > 2
            continue
        end
        i += 1
    end
    points_S = points_S[1:3,:]
    points_S = points_S + sigma_n*randn(size(points_S))

    # Verify poses
    C = cov(points_RB');    D,V = eig(C)
    nhat_points = V[:,1]
    C = cov(lines_RB');    D,V = eig(C)
    nhat_lines = V[:,1]
    if 1-abs(vecdot(nhat_points,nhat_lines)) > 0.001 || D[1,1] > 1e-5
        warn("Something is wrong")
        @show 1-abs(vecdot(nhat_points,nhat_lines)), D[1,1]
    end

    points = zeros(4,N_poses)
    for i = 1:N_poses
        T_RB_S = T_RB_TF[:,:,i] * T_TF_S
        #         points(:,end+1) = T_RB_S*[(points_S[:,i]);1];
        points[:,i] = T_RB_S*[(points_S[:,i]);1]
        #         points(:,end+1) = inv(T_RB_S)*[(points_S[:,i] - 0.01*lines_S[:,i]);1];
    end
    if any(abs(points[1:3,:] - points_RB) .> 1e-2)
        warn("Something is wrong 2")
    end
    points_S,lines_S, T_RB_TF
end

srand(1)
MC = 100
N_planes = 3
N_points = 2
N_calibs = 30
N_poses = 10
sigma_n = 5e-4
SSEStart = zeros(MC)
normStart = zeros(MC)
distStart = zeros(MC)
rotStart = zeros(MC)
SSEMC = zeros(MC,N_calibs)
normMC = zeros(MC,N_calibs)

for mc = 1:MC
    # Simulate measurement data
    T_TF_S = Rt2T(rpy2R(randn(3)),0.5randn(3)) #TCP to sensor
    points_S = Array(Float64,3,0);    lines_S = Array(Float64,3,0);    T_RB_TF = Array(Float64,4,4,0);    plane = Plane[]; planes = Array(Int,0)
    SSE = zeros(N_planes)
    for j = 1:N_planes
        push!(plane, generateRandomPlane(j))
        points_St,lines_St, T_RB_TFt = generateRandomPoses(N_poses, sigma_n,T_TF_S, plane[j])
        points_S = [points_S points_St]
        lines_S = [lines_S lines_St]
        T_RB_TF = cat(3,T_RB_TF, T_RB_TFt)
        planes = [planes; j*ones(Int,N_poses)]
        SSE[j] = pointDiff(T_TF_S,T_RB_TFt,points_St)[1]
    end
    @edb any(abs(SSE .> 1e-3)) && warn("Points does not seem to lie on a plane")
    dist = [norm(T_RB_TF[1:4,4,i] - T_TF_S*[points_S[:,i];1]) for i in 1:size(T_RB_TF,3)]
    #     hist(dist,15)

    ## Calibration
    # display('----------------------------------------')
    # Generate nominal transformation
    N_poses = size(points_S,2)
    T_TF_S_real = T_TF_S
    if true
        T_TF_S0 = Rt2T(rpy2R(60π/180*(rand(3)-0.5)),0.4*(rand(3)-0.5))*T_TF_S_real
        #         T_TF_S0 = rt2tr(rpy2R(2*(rand(1,3)-0.5),'deg'),0.01*(rand(3,1)-0.5))*T_TF_S_real;
    else
        T_TF_S0 = T_TF_S_real
    end

    for j = 1:N_planes
        ind = planes .== j
        ind = find(ind)
        SSE[j] = sqrt(pointDiff(T_TF_S0,T_RB_TF[:,:,ind],points_S[1:3,ind])[1])
    end

    SSEStart[mc] = mean(SSE)
    normStart[mc] = vecnorm(T_TF_S0-T_TF_S_real)
    distStart[mc] = norm(T_TF_S0[1:3,4]-T_TF_S_real[1:3,4])
    rotStart[mc] = norm(180/π*R2rpy(T_TF_S0\T_TF_S_real))

    T_TF_S, RMScalibs, norms = calibNAXP(points_S, lines_S, T_RB_TF, T_TF_S0, planes,  N_calibs)#, T_TF_S_real)
    SSEMC[mc,:] = RMScalibs
    normMC[mc,:] = norms

    if RMScalibs[end] > 1e-2
        println("Bad result")
    end

    distEnd[mc] = norm(T_TF_S[1:3,4]-T_TF_S_real[1:3,4])
    rotEnd[mc] = norm(180/π*R2rpy(T_TF_S\T_TF_S_real))
end
# normStart
# normMC

plot(0:N_calibs,[normStart', normMC]',yscale=:log10,c=:black, xlabel="Number of iterations")
plot!([0 N_calibs], [sigma_n sigma_n],l=:dash, c=:red)

plot(0:N_calibs,[SSEStart', SSEMC]',yscale=:log10,c=:black, xlabel="Number of iterations",title="RMS distance from points to plane [m]")
plot!([0 N_calibs], [sigma_n sigma_n],l=:dash,c=:red)

boxplot(["Before" "After"],([distStart', distEnd']), title="Distance error [m]")
# hold on, plot([0 3], [sigma_n sigma_n],'--r')

boxplot(["Before" "After"],([rotStart', rotEnd']),title="Rotation error [degree]")


text:
### A Pluto.jl notebook ###
# v0.19.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ a2f22b9d-13f8-4e39-b7dd-14b905d987ab
#add button to trigger presentation mode
html"<button onclick='present()'>present</button>"

# ╔═╡ 731c88b4-7daf-480d-b163-7003a5fbd41f
begin
	html"""

  FINC 462/662 -- Fixed Income Securities

	"""
end

# ╔═╡ a5de5746-3df0-45b4-a62c-3daf36f015a5
begin
	html"""

 FINC-462/662: Fixed Income Securities

  The U.S. Treasury Market

 Spring 2022

 Prof. Matt Fleckenstein

 University of Delaware,
	Lerner College of Business and Economics

	"""
end

# ╔═╡ 69dca09e-b572-436f-90c6-187af17cf027
TableOfContents(aside=true, depth=1)

# ╔═╡ 6498b10d-bece-42bf-a32b-631224857753
md"""
# Overview
"""

# ╔═╡ 95db374b-b10d-4877-a38d-1d0ac45877c4
begin
	html"""
	<fieldset>
        <legend>Our goals for today</legend>

        <input type="checkbox" value="">Know the different types of securities issued by the U.S. Department of the Treasury

        <input type="checkbox" value="">Be familiar with the Treasury auction process.

	  	</fieldset>
	"""
end

# ╔═╡ 10a73b76-b6ad-4081-b8b1-21363608618e
md"""
# The U.S. Treasury Market
- The Department of the Treasury is the largest single issuer of debt in the world.
- The large volume of total debt and the large size of any single issue have contributed to making the Treasury market the most active and hence the most liquid market in the world

[SIFMA Fixed Income Statistics](https://www.sifma.org/resources/research/fixed-income-chart)
"""

# ╔═╡ 42f7e15d-9abc-40ca-bd70-7131957907a5
md"""
- Why the Treasury Market matters ...
  - [WSJ, February 1, 2022: U.S. National Debt Exceeds $30 Trillion for First Time](https://www.wsj.com/articles/u-s-national-debt-exceeds-30-trillion-for-first-time-11643766231?mod=hp_lead_pos3)
"""

# ╔═╡ b33e1473-e374-441b-a6a4-6c2a7b079c98
md"""
## U.S. Treasury Debt Amount Outstanding
Year-end data in trillions of dollars. \
_Source: SIFMA_
"""

# ╔═╡ 7bb29c12-c0a0-4eb0-a4d6-bd9c52443751
#Treasury Market Amounts Outstanding
begin
	FI_Out = DataFrame()
	plot_Outstanding = plot()
	let
		file="./US-Fixed-Income-Securities-Statistics-SIFMA.xlsx"
		sheet="Outstanding"
		cols="A:I"
		startrow = 8
		xlsxFile = XLSX.readtable(file, sheet, cols; first_row=startrow, header=true, 			infer_eltypes=true)
		FI_Out = DataFrame(xlsxFile...)

		catSelect="Treasury"
		println(catSelect)
		plotData = select(FI_Out, "Year","$(catSelect)")
		rename!(plotData,"Year" => :x, "$(catSelect)" => :y)
		dropmissing!(plotData)
		transform!(plotData, :y => ByRow(x->x/1000), renamecols=false)
		minX = 1975
		maxX = 2025
		minY = 0.0
		maxY = 25.0
		plot!(plot_Outstanding, plotData.x,plotData.y,
			xlim=(minX,maxX), ylim=(minY, maxY),
			ylabel="Trillions of Dollars",label="$(catSelect)",
			legend = :topleft, title="Amount Outstanding")
		plot(plot_Outstanding)
	end
end

# ╔═╡ f89f6fc4-cb6f-4d7b-84c9-f8a18e764427
md"""
##
"""

# ╔═╡ 313959e0-eea6-449a-a6b8-29c9394e50b6
LocalResource("./TreasuryOutstandingSIFMA.svg",:width => 900)

# ╔═╡ e20b95fa-e038-4d3b-82ca-8a38262c64c4
md"""
[Link to SIFMA](https://www.sifma.org/resources/research/fixed-income-chart/)
"""

# ╔═╡ 2bf62523-a0f1-400e-9178-79062e0a171e
md"""
## U.S. Treasury Debt Amount Issued
Annual data at year-end in trillions of dollars.\
_Source: SIFMA_
"""

# ╔═╡ 206df1a6-bf12-4957-93f7-472ee89594c4
#US Treasury Securities Amounts Issued
begin
	FI_Iss = DataFrame()
	plot_Issuance = plot()
	let
		file="./US-Fixed-Income-Securities-Statistics-SIFMA.xlsx"
		sheet="Issuance"
		cols="A:H"
		startrow = 8
		xlsxFile = XLSX.readtable(file, sheet, cols; first_row=startrow,
			header=true, infer_eltypes=true)
		FI_Iss = DataFrame(xlsxFile...)

		catSelect="Treasury"
		plotData = select(FI_Iss, "Year","$(catSelect)")
		rename!(plotData,"Year" => :x, "$(catSelect)" => :y)
		dropmissing!(plotData)
		transform!(plotData, :y => ByRow(x->x/1000), renamecols=false)
		minX = 1994
		maxX = 2021
		minY = 0.0
		maxY = 4.0
		plot!(plot_Issuance, plotData.x,plotData.y,
			xlim=(minX,maxX), ylim=(minY, maxY),
			ylabel="Trillions of Dollars",label="$(catSelect)",
			legend = :topleft, title="Amount Issued")
		plot(plot_Issuance)
	end
end

# ╔═╡ c8f703dc-5a1c-4c39-b9a4-c3dd0a2806d8
md"""
##
"""

# ╔═╡ 13a58974-2ae8-4596-933f-5d2b5b8f5a23
LocalResource("./TreasuryIssuanceSIFMA.svg",:width => 900)

# ╔═╡ c7b5fcd2-a727-4bde-aae3-cd39dfb9a43b
md"""
##
"""

# ╔═╡ e89beeee-1ad9-4440-ba95-9bd0eb6bfe72
md"""
_Note_: Issuance is long-term instruments only.\
[Link to SIFMA](https://www.sifma.org/resources/research/fixed-income-chart/)
"""

# ╔═╡ 7f91ba74-ccca-4f08-8252-9de7f956b285
LocalResource("./TreasuryIssuanceSIFMA.png",:width => 900)

# ╔═╡ 84009ee2-cfc3-437c-8782-d370b8b8c4f6
md"""
##
"""

# ╔═╡ a2df8dc1-f6da-4816-85b9-c7caadc03584
#Treasury Market Total Amount Outstanding vs Equity
begin

	plot_MktCap = plot()
	let
		file="./MarketCaps.xlsx"
		sheet="Sheet1"
		cols="A:B"
		startrow = 8
		xlsxFile = XLSX.readtable(file, sheet, cols; first_row=startrow, header=false,
			column_labels=["Date","Equity"], infer_eltypes=true)
		EquMktCap = DataFrame(xlsxFile...)
		dropmissing!(EquMktCap)
		transform!(EquMktCap,:Equity=>(x->x./1000)=>:Equity,
			:Date => (x->year.(x)) => :Year,
			:Date => (x->month.(x)) => :Month)
		filter!(:Month => ==(12),EquMktCap)
		MktCaps = leftjoin(EquMktCap,FI_Out[:,[:Year,:Treasury]], on=:Year)
		select!(MktCaps,:Year,:Treasury,:Equity)
		transform!(MktCaps,names(MktCaps,Not(:Year)) .=> (x->x./1000.0),
			renamecols=false)
		minX = minimum(MktCaps.Year)-1
		maxX = maximum(MktCaps.Year)+1
		minY = 0.0
		maxY = 50.0
		plot!(plot_MktCap, MktCaps.Year,MktCaps.Treasury,
			  xlim=(minX,maxX), ylim=(minY, maxY),
			  ylabel="Trillions of Dollars", label="Fixed Income",
			  legend = :topleft, title="US Treasury vs. Equity Markets")
		plot!(plot_MktCap, MktCaps.Year,MktCaps."Equity",
			  label="Equity")
     	plot(plot_MktCap)
	end
end

# ╔═╡ 5db80dee-c643-4cc1-bca4-00b1ffd0a948
md"""
## U.S. Treasury Securities Trading Volume
Average Daily Trading Volume. Annual data at year-end in billions of dollars.\
_Source: SIFMA_
"""

# ╔═╡ 0357258e-8133-4ef3-946f-5134a9560b9d
#US Treasury Market Trading Volumes
begin
	FI_Vol = DataFrame()
	plot_Vol = plot()
	let
		file="./US-Fixed-Income-Securities-Statistics-SIFMA.xlsx"
		sheet="Trading Volume"
		cols="A:I"
		startrow = 8
		xlsxFile = XLSX.readtable(file, sheet, cols; first_row=startrow,
			header=true, infer_eltypes=true)
		FI_Vol = DataFrame(xlsxFile...)

		catSelect="Treasury"
		plotData = select(FI_Vol, "Year","$(catSelect)")
		rename!(plotData,"Year" => :x, "$(catSelect)" => :y)
		dropmissing!(plotData)
		minX = 1995
		maxX = 2021
		minY = 0.0
		maxY = 650.0
		plot!(plot_Vol, plotData.x,plotData.y,
			xlim=(minX,maxX), ylim=(minY, maxY),
			ylabel="Billions of Dollars",label="$(catSelect)",
			legend = :topleft, title="Trading Volume")
		plot(plot_Vol)
	end
end

# ╔═╡ ae67f236-60ec-4c2e-abd3-9ea537200ca3
md"""
# Types of Treasury Securities
"""

# ╔═╡ fa1dcebd-4a54-433b-897c-de2929fb8233
md"""
>- Where to get information about U.S. Treasury securities?
>- Go to webpage of the [U.S. Treasury](https://www.treasurydirect.gov/).
>- In the middle panel, click on `Treasury securities Overview`.
"""

# ╔═╡ 2452e56b-2a61-4608-9ca0-55c8cdca89c6
LocalResource("TreasuryDirect_01.png",:width => 900)

# ╔═╡ f688e93c-e8b1-4718-b56e-d76bf70eef57
md"""
##
"""

# ╔═╡ 4914c857-5dc8-43b8-8e70-377debbc3d0c
LocalResource("TreasuryDirect_02.png",:width => 900)

# ╔═╡ d25744ad-fbe3-4f08-b17a-d00d51bcb2b8
md"""
# Treasury Bill (T-bill)
- Short-term securities with maturities of 4, 13, 26, and 52 weeks.
- Treasury bills do not pay interest before maturity.
  - This is often referred to as a _discount security_ or _zero-coupon security_.
- Instead, Treasury bills are issued at a price less than their par value and at maturity, Treasury bills pay back their par value.
  - Intuitively, the “interest” to the investor is the difference between par value and the purchase price.
  - Example: a 52-week T-bill with par value of \$100 has a price of \$98.
"""

# ╔═╡ 382a9c58-4baa-4ad2-9f32-7189d4881c50
md"""
## Example of 52-week Treasury Bill
"""

# ╔═╡ deb54591-8769-480d-98ad-32157abac2c0
LocalResource("./BloombergTBill_01.png", :width=>900)

# ╔═╡ 3578793a-e0b6-42de-941b-2db09ef99e49
md"""
>- How to get there on the Bloomberg terminal?
>  - Open a terminal and on the keyboard type `T Bill`.
>  - In the popup window,  select `B Govt - United States Treasury Bill (Multiple Matches)`.
>  - Next, click on one of the different Treasury bills in the list.
>  - Then, click on `DES` on the top-right, or type `DES` on the keyboard and press enter.
"""

# ╔═╡ 82db2847-a1a1-45a3-8d29-c3ac4213591c
md"""
##
"""

# ╔═╡ adbc8a60-c33c-44ba-a284-fa7b676b9618
md"""
- Let's look at price quotes for this Treasury bill.
- If you were to purchase this Treasury bill, would you pay \$0.8450?
- The answer is _no_. Prices for Treasury bills are quoted in a specific way, which we will discuss in the next lecture.
"""

# ╔═╡ 1a6add34-52e1-4663-9228-f0fca63f8cc2
LocalResource("./BloombergTBill_02.png")

# ╔═╡ 8de081f6-2b7f-4d03-b565-b4ac01c2773e
md"""
>- How to get there on the Bloomberg terminal?
>  - From the `Description` page where you are currently at, type `GP` and press enter.
"""

# ╔═╡ d296b99d-588d-48ac-8e0c-a0bcbde72c56
md"""
# Treasury Note (T-Note)
- Medium-term securities with maturities of 2, 3, 5, 7, and 10 years.
- Treasury notes pay interest every six months up to and including the maturity date.
  - Example: A 2-year T-note has its last interest payment in two years, and it pays interest after 6 months, 12 months, and 18 months.
- At maturity, Treasury notes pay back their par value.
"""

# ╔═╡ 1c4db2fa-7187-4966-aa3f-9a7f56036497
md"""
## Example of 10-year Treasury Note
"""

# ╔═╡ 3a3b7641-ed63-4cce-a98c-a4cf1b6e603d
LocalResource("./BloombergTNote_01.png")

# ╔═╡ 2384b88b-e8dc-4e9c-9488-27828da0289d
md"""
>- How to get there on the Bloomberg terminal?
>  - Open a terminal and on the keyboard type `T Note`.
>  - In the popup window,  select `T Govt - United States Treasury Note/Bond (Multiple Matches)`.
>  - Next, click on one of the different Treasury notes in the list.
>  - Then, click on `DES` on the top-right, or type `DES` on the keyboard and press enter.
"""

# ╔═╡ 0918383a-9a0b-4f8c-ba63-aaad1cf6d72c
md"""
- Let's look at price quotes for this Treasury note.
- If you were to purchase this Treasury note, would you pay \$95.02?
- The answer is _no_. Prices for Treasury notes are quoted in a specific way, which we will discuss in the next lecture.
"""

# ╔═╡ c9b876e2-cb29-4a70-9c1f-abb2dedf3fe2
LocalResource("./BloombergTNote_02.png")

# ╔═╡ d912c3cb-c1a6-46fe-af9d-25a75074b31d
md"""
>- How to get there on the Bloomberg terminal?
>  - From the `Description` page where you are currently at, type `GP` and press enter.
"""

# ╔═╡ 40df1177-48cf-4cb2-9860-11d2d9f0a450
md"""
# Treasury Bond (T-Bond)
- Long-term securities with maturities of 20 and 30 years.
  - Currently, the Treasury does not issue 15-year Treasury bonds.
- Treasury bonds notes pay interest every six months up to and including the maturity date. At maturity, Treasury notes pay back their par value.
  - Similar to Treasury notes.

[Treasury Marketable Securities](https://www.treasurydirect.gov/instit/marketables/marketables.htm)
"""

# ╔═╡ 4567e3a2-24fb-494f-8005-801b03ffeabf
#https://www.treasurydirect.gov/indiv/research/history/histtime/histtime_bonds.htm

# ╔═╡ 76774fcb-f806-4db2-b653-fe69720f43a7
md"""
## Coupon bonds
- Treasury notes and bonds are referred to as **coupon** securities.
- Why?
"""

# ╔═╡ a440a28f-04af-4248-8f2b-c8e6a726c07a
LocalResource("TreasuryCouponPicture.svg",:width => 900)

# ╔═╡ fef5577e-025f-4afc-9fd9-aa025146858e
md"""
# Treasury Floating Rate Note (FRN)
- First issued in 2014 by the U.S. Treasury.
- Maturity of 2 years.
- Pay interest every three months up to and including the maturity date.
  - At maturity, FRNs pay back their par value.
- The interest on an FRN varies with interest rate on 13-week Treasury bills.
"""

# ╔═╡ 90883b07-064f-4395-a805-ddaaad7f97fd
md"""
## Example of Treasury Floating Rate Note (FRN)
"""

# ╔═╡ faecc806-5f33-49c6-a463-b349cbfb8acc
LocalResource("./BloombergFRN_01.png")

# ╔═╡ c7a254a6-d877-4de0-a3e5-a64fb93ae079
md"""
>- How to get there on the Bloomberg terminal?
>  - Open a terminal and on the keyboard type `Treasury Floating Rate Note`.
>  - In the popup window,  select `TF Govt - United States Treasury Floating Rate Note (Multiple Matches)`.
>  - Next, click on one of the different Treasury FRNs in the list.
>  - Then, click on `DES` on the top-right, or type `DES` on the keyboard and press enter.
"""

# ╔═╡ 9b53686c-e362-46c5-bdd9-0c388d9cd122
LocalResource("./BloombergFRN_02.png")

# ╔═╡ 9d0017a5-6b7b-4c1b-9323-7a8fdc54db77
md"""
- Let's look at price quotes for this Treasury FRN.
- If you were to purchase this Treasury FRN, would you pay \$100.1674?
- Generally, the answer is _no_. Prices for Treasury FRNs are quoted in a specific way, which we will discuss in the next lecture.
"""

# ╔═╡ 69d8c85f-3735-4bd2-bd37-1ddadab87140
md"""
>- How to get there on the Bloomberg terminal?
>  - From the `Description` page where you are currently at, type `GP` and press enter.
"""

# ╔═╡ cf1fa1d8-da24-4f48-8762-affa1f58279a
md"""
# Treasury Inflation Protected Securities (TIPS)
- First issued in 1997 by the U.S. Treasury.
- Maturities of 5, 10, and 30 years.
- TIPS pay interest every six months up to and including the maturity date. At maturity, Treasury notes pay back their par value.
  - Similar to Treasury notes and bonds.
- Key difference is that both par value and interest go up with the rate of inflation.
"""

# ╔═╡ 3000e0f6-381f-4a9c-85c1-c517af617f14
md"""
- Why inflation matters and why one might want to hedge against inflation.
  - [WSJ, January 12, 2022: U.S. Inflation Hit 7% in December, Fastest Pace Since 1982](https://www.wsj.com/articles/us-inflation-consumer-price-index-december-2021-11641940760?mod=article_inline)
  - [WSJ, February 10, 2022: U.S. Inflation Rate Accelerates to a 40-Year High of 7.5%](https://www.wsj.com/articles/us-inflation-consumer-price-index-january-2022-11644452274?mod=hp_lead_pos1)
  - [WSJ, February 6, 2022: What Investors Should Know About TIPS](https://www.wsj.com/articles/tips-what-investors-should-know-treasury-inflation-protected-securities-11643849892?mod=hp_lista_pos1)
"""

# ╔═╡ 5c8f4910-7b54-414a-8eb9-52104c0ff36e
md"""
## Example of a Treasury TIPS
"""

# ╔═╡ 972e44da-4614-4349-bcb7-edb0518fe24d
LocalResource("./BloombergTIPS_01.png")

# ╔═╡ a1c74635-2fb9-492c-8243-4232b771718b
md"""
>- How to get there on the Bloomberg terminal?
>  - Open a terminal and on the keyboard type `Treasury TIPS`.
>  - In the popup window,  select `TII Govt - United States Treasury Inflation Indexed Bonds (Multiple Matches)`.
>  - Next, click on one of the different Treasury bills in the list.
>  - Then, click on `DES` on the top-right, or type `DES` on the keyboard and press enter.
"""

# ╔═╡ 45ffcfef-16f3-464a-b868-0d3f64cb6019
md"""
##
"""

# ╔═╡ 874afd1d-eba3-4295-9f96-7502c187b171
md"""
- Let's look at price quotes for this Treasury TIPS.
- If you were to purchase this Treasury TIPS, would you pay \$105.28?
- The answer is _no_. Prices for Treasury TIPS are quoted in a specific way, which we will discuss in the next lecture.
"""

# ╔═╡ 95110a89-1d2b-4c83-9b1c-d7dc95ab89c8
LocalResource("./BloombergTIPS_02.png")

# ╔═╡ aa70d2f6-c726-45cc-89cb-bb6afb85f065
md"""
>- How to get there on the Bloomberg terminal?
>  - From the `Description` page where you are currently at, type `GP` and press enter.
"""

# ╔═╡ d69cce2c-dce4-4369-b3db-968c0b4d3dc6
md"""
## Inflation-Linked notes date back centuries
"""

# ╔═╡ ccba5159-b6bb-4c98-8f04-8dd91c032146
md"""
> Both Principal and Interest to be paid in the then current Money of said State, in a greater or less Sum, according as Five Bushels of Corn, Sixty-eight Pounds and four-seventh Parts of a Pound of Beef, Ten Pounds of Sheeps Wool, and Sixteen Pounds of Sole Leather shall then cost more or less than One Hundred and Thirty Pounds current Money, at the then current Prices of said Articles.
_Source: “Inflation-indexed Securities: Bonds, Swaps and Other Derivatives”, 2nd Edition, M. Deacon, A. Derry, D. Mirfendereski, Wiley._

"""

# ╔═╡ f50d144b-6d62-4b3f-b85a-d76861ea2ae0
LocalResource("./SoldiersDepreciationNote.svg",:width => 900)

# ╔═╡ c4aa9ab9-5502-4676-8b6a-5b55cc645e03
md"""
## TIPS Inflation-Adjustment
"""

# ╔═╡ f481451a-6949-4604-811b-6b4c530fd644
LocalResource("./TIPSInflationAdjustment.svg",:width => 900)

# ╔═╡ 992074c6-d6d5-4417-9edc-0229a0333070
md"""
# Treasury STRIPS

- The Treasury does not issue zero-coupon notes or bonds.
- However, because of the demand for zero-coupon instruments with no credit risk, the private sector has created such securities.
- The process of separating the interest on a bond from the underlying principal is called coupon stripping
- Zero-coupon Treasury securities were first created in August 1982 by large Wall-Street firms.
  - The problem with these securities was that they were identified with particular dealers and therefore reduced liquidity.
  - Moreover, the process involved legal and insurance costs.
  - Today, all Treasury notes and bonds (fixed-principal and inflation-indexed) are eligible for stripping.
- The zero-coupon Treasury securities created under the STRIPS program are direct obligations of the U.S. government
"""

# ╔═╡ cc282e0d-9faa-4808-b506-335e9097f2fe
md"""
## Treasury STRIPS in Bloomberg
"""

# ╔═╡ d44298e3-5763-4249-8140-5eefb80adb49
LocalResource("./TreasurySTRIPSBloomberg.png",:width => 900)

# ╔═╡ 57759a33-f026-4c15-91ba-f6af45627f1e
md"""
>- How to get there on the Bloomberg terminal?
>  - Open a terminal and on the keyboard type `Treasury STRIP`.
>  - In the popup window,  select `S Govt - United States Treasury Strip Coupon (Multiple Matches)`.
>  - Next, click on one of the different Treasury bills in the list.
>  - Then, click on `DES` on the top-right, or type `DES` on the keyboard and press enter.
"""

# ╔═╡ 9c596473-76f6-4e9a-b3dd-32add39eaf92
md"""
##
"""

# ╔═╡ 9db71771-d778-4d43-af7a-80d8cfff97ff
LocalResource("./BloombergSTRIPS_01.png")

# ╔═╡ 36678746-8926-451d-b692-b57c620170c5
md"""
##
"""

# ╔═╡ 8f65e51f-01d6-4d02-bee2-f7d16f62a6c2
md"""
- Let's look at price quotes for this Treasury STRIPS.
"""

# ╔═╡ e16d631b-befe-4657-ad03-e65b5999e49a
LocalResource("./BloombergSTRIPS_02.png")

# ╔═╡ 0f801b78-6b89-48f4-ad96-e7336094e024
md"""
# Treasury Auctions
"""

# ╔═╡ c1165da4-1149-42d2-8841-a38d7cbc131d
LocalResource("./WSJ_TreasuryAuctionHeadline.png",:width => 900)

# ╔═╡ 7759dee0-da42-4796-ae60-e8bbe507f52e
md"""
_Source_: The Wall Street Journal, March 10, 2021.\
[Link to Article](https://www.wsj.com/articles/flood-of-new-debt-tests-bond-market-11615372201?mod=Searchresults_pos15&page=2)
"""

# ╔═╡ a8973960-ab9f-435b-a089-a6da9b4e24ec
LocalResource("./WSJ_TreasuryAuctionSizes.svg",:width => 900)

# ╔═╡ db0a2d77-11b1-4f54-819c-d43b0e18002b
md"""
## Treasury Auctions Basics

- Treasury securities are sold in the primary market through sealed-bid auctions.
- Each auction is **announced** several days in advance by means of a Treasury Department press release or press conference.
- The **announcement** provides details of the offering, including the _offering amount_ and the _term_ and _type_ of security being offered, and describes some of the auction rules and procedures.
"""

# ╔═╡ bb28a8c2-aaa8-4400-8601-514e156a9fe6
md"""
> Where to get information about the Treasury auction process?
> - Go to [TreasuryDirect.gov](https://www.treasurydirect.gov/)
> - In the middle panel, you see multiple links to resources about Treasury auctions.
> - For example, to learn about the auction process, click on `How Treasury auctions work`
"""

# ╔═╡ c3a14721-d694-4302-b8b2-adf887e6fb7d
LocalResource("./TreasuryDirect_Auctions_01.png")

# ╔═╡ 03393752-d39f-45a3-8157-520f361058ab
md"""
## Treasury Auction Announcement
"""

# ╔═╡ fea57e60-5e7b-4b5d-84a3-f16ca0279799
LocalResource("./TreasuryAuctionAnnouncement.png",:width => 900)

# ╔═╡ 0608bb54-8495-4365-810c-5f9010f0eed6
md"""
[Treasury Auction Announcements](https://www.treasurydirect.gov/instit/annceresult/press/press_secannpr.htm)
"""

# ╔═╡ c85dc38a-af8c-4c79-821c-2e9c90399318
md"""
## Treasury Auction Process

- The auction for Treasury securities is conducted on a competitive bid basis.
- There are two types of bids that may be submitted by a bidder:
  - noncompetitive bids
  - competitive bids
- A noncompetitive bid is submitted by an entity that is willing to purchase the auctioned security at the yield that is determined by the auction process

"""

# ╔═╡ c1c539c2-afe1-41f4-9c26-29b5a0acfba1
md"""
## Noncompetitive Bids
- When a noncompetitive bid is submitted, the bidder only specifies the **quantity** sought.
  - The quantity in a noncompetitive bid may not exceed $5 million.
- For more information, visit  [How Treasury Auctions Work](https://www.treasurydirect.gov/instit/auctfund/work/work.htm).
"""

# ╔═╡ 35f529cd-9642-4bd9-a2f7-90018d1936be
md"""
## Treasury Auction Results - Noncompetitive Bids
"""

# ╔═╡ 094387ed-e45a-4d07-a4e6-0d3f08fda046
LocalResource("./TreasuryAuctionResultsNoncompetitive.png",:width => 900)

# ╔═╡ 43e798f2-a8de-4920-9e44-9c56ba965670
md"""
## Competitive Bids
- A competitive bid specifies both the **quantity** sought _and_ the **yield** at which the bidder is willing to purchase the auctioned security.
- The competitive bids are then arranged from the __lowest yield__ bid to the __highest yield__ bid submitted.
  - The highest yield accepted by the Treasury is referred to as the **stop-out yield** (or **high yield**).
"""

# ╔═╡ d8e7db96-9471-424d-a616-2c1327d23c77
md"""
## Treasury Auction Results - Competitive Bids
"""

# ╔═╡ 96a2e1f6-b15d-446f-8828-c4c099fca07d
LocalResource("./TreasuryAuctionResultsCompetitive.png",:width => 900)

# ╔═╡ 940baa2d-4117-4961-95a9-2043f3adc374
md"""
## [Additional Treasury Auction information](https://www.treasurydirect.gov/instit/auctfund/work/work.htm)
"""

# ╔═╡ a84b0684-0388-47ac-a338-2a52b5254a75
LocalResource("./TreasuryWhenIssuedTrading.png",:width => 1200)

# ╔═╡ 57e22195-f876-4bab-a346-c1221a7607a3
md"""
## Auction Schedule

[Upcoming Treasury Auctions](http://www.treasury.gov/resource-center/data-chart-center/quarterly-refunding/Documents/auctions.pdf)
"""

# ╔═╡ 7792e937-3487-4d29-b2ba-c1b1d335e4b0
md"""
|                       Instrument                      	|              Auction Frequency             	|
|:-----------------------------------------------------:	|:------------------------------------------:	|
|                      4-week Bills                     	|             Weekly   (Tuesdays)            	|
|               13-week Bills               	|             Weekly   (Mondays)             	|
|               26-week Bills               	|             Weekly   (Mondays)             	|
|                     52-week Bills                     	|         Every   4 weeks (Tuesdays)         	|
|                      2-year Notes                     	|          Monthly   (End of month)          	|
|                      3-year Notes                     	|         Monthly   (Middle of month)        	|
|                      5-year Notes                     	|          Monthly   (End of month)          	|
|                      7-year Notes                     	|          Monthly   (End of month)          	|
|                     10-year Notes                     	|         Monthly   (Middle of month)        	|
|                     30-year Bonds                     	|         Monthly   (Middle of month)        	|
| 5-year TIPS  	|   Three times per year (Apr, Aug, Dec)   	|
|                      10-year TIPS                     	| Bimonthly   (Jan, Mar, May, Jul, Sep, Nov) 	|
|                      30-year TIPS                     	|   Three   times per year (Feb, Jun, Oct)   	|
|            2-year FRN             	|          Monthly   (End of month)          	|
"""

# ╔═╡ f2324049-7f96-4ad6-9eb6-6f8856fdf6fc
md"""
## Treasury Auction Example
"""

# ╔═╡ 22e9ecfc-72c3-43e1-9214-b907c5cd7d07
md"""
- Treasury Department announces a \$11 billion offering.
- There are \$1 billion in noncompetitive bids received.
- Thus, there is \$11 - $1 = \$10 billion left to be allocated to competitive bidders.
- Suppose 6 competitive bidders submit the following sealed bids.
"""

# ╔═╡ 86254e63-de75-42c9-9b1e-e0a414cd0f63
md"""
| Name     	| Yield 	| Amount       	|
|----------	|-------	|--------------	|
| Bidder 1 	| 2.998% 	| \$3.5 billion 	|
| Bidder 2 	| 2.999% 	| \$2.5 billion 	|
| Bidder 3 	| 3.000% 	| \$3.0 billion   	|
| Bidder 4 	| 3.000% 	| \$3.0 billion   	|
| Bidder 5 	| 3.001% 	| \$2.0 billion   	|
| Bidder 6 	| 3.002% 	| \$1.0 billion   	|

**Question**: What do Bidder 1 to Bidder 6 get awarded in the auction?
"""

# ╔═╡ adeb0d3f-602d-4401-8a15-09b8d885ef40
md"""
##
"""

# ╔═╡ 290c5c6d-f2ab-4707-a298-99e9d678a72e
md"""
→ The Treasury works its way down the list of competitive bids and accepts the total amount submitted at the lowest possible bid yields (hightest prices) until the full offering amount has been awarded.
"""

# ╔═╡ 785bff43-7920-494d-b59c-2fa6c8cd2391
md"""
| Name     	| Yield 	| Amount Bid   	| Amount Awarded   |  Competitive Amount  |
|----------	|-------	|--------------	| ---------------- |  ------------------- |
|           |           |               |                  |  **\$10 billion**    	  |
| Bidder 1 	| 2.998% 	| \$3.5 billion |                  |  -\$3.5 billion      |
|    	 	| 		 	| 				|                  |  **\$6.5 billion**   |
| Bidder 2 	| 2.999% 	| \$2.5 billion |                  |  -\$2.5 billion      |
|    	 	| 		 	| 				|                  |  **\$4.0 billion**   |
"""

# ╔═╡ 6a64ec0a-6f67-428c-afbf-b67a4ae3568b
md"""
##
"""

# ╔═╡ 1a7a0375-507e-4a5b-befc-c88208ab8394
md"""
- At this point there are \$4 billion remaining for competitive bidding.
- However, there are a total of \$6 billion in bids at the next lowest rate (3.000%) by Bidder 3 and Bidder 4.
- The highest accepted rate (3.000%) is known as the **stop-out rate**.
- In this case, each bidder at this rate is awarded a percentage of their total bid amount.
  - The **allocation percentage** is calculated by dividing the remaining competitive offering by the total amount bid at the stop-out rate.

$$\frac{\textrm{Remaining Comptetitive Offering}}{\textrm{Total Bids at Stop-Out Rate}} = \frac{\$4 \textrm{ billion}}{\$6\textrm{ billion}} = 66.7\%$$

- Bidder 3 and Bidder 4 each get a partial allocation of \$2 billion (66.67% x \$3.0 billion)
"""

# ╔═╡ 9c1a34f3-2ad5-4f14-b77b-a941fdff79bf
md"""
##
"""

# ╔═╡ 12aea81d-c470-4fd5-b769-46040c3676e2
md"""
|        Name       	|     Yield    	|       Amount Bid      	|     Amount awarded    	|     Allocation percentage    	|     Rate awarded    	|
|:-----------------:	|:------------:	|:---------------------:	|:---------------------:	|:----------------------------:	|:-------------------:	|
|     Bidder   1    	|     3.00%    	|     \$3.5   billion    	|     \$3.5   billion    	|              100%            	|         3.00%       	|
|     Bidder   2    	|     3.00%    	|     \$2.5   billion    	|     \$2.5   billion    	|              100%            	|         3.00%       	|
|     Bidder   3    	|     3.00%    	|     \$3.0   billion    	|     \$2.0   billion    	|             66.67%           	|         3.00%       	|
|     Bidder   4    	|     3.00%    	|     \$3.0   billion    	|     \$2.0   billion    	|             66.67%           	|         3.00%       	|
|-----------------	|------------	|---------------------	|---------------------	|----------------------------	|-------------------	|
|     Bidder   5    	|     3.00%    	|     \$2.0   billion    	|           \$0          	|               0%             	|          N/A        	|
|     Bidder   6    	|     3.00%    	|     \$1.0   billion    	|           \$0          	|               0%             	|          N/A        	|
"""

# ╔═╡ 4694dcaf-e8c4-4c99-9020-25f7227c4fd3
md"""
## General Auction Timing
[General Auction Timing](https://www.treasurydirect.gov/instit/auctfund/work/auctime/auctime.htm)
"""

# ╔═╡ 53c77ef1-899d-47c8-8a30-ea38380d1614
md"""
## Wrap-Up
"""

# ╔═╡ 670e45a3-9d28-47ae-a6b6-a1b1c67a0a4c
begin
	html"""
	<fieldset>
        <legend>Our goals for today</legend>

        <input type="checkbox" value="" checked>Know the different types of securities issued by the U.S. Department of the Treasury

        <input type="checkbox" value="" checked>Be familiar with the Treasury auction process.

	  	</fieldset>
	"""
end

# ╔═╡ 2ee2c328-5ebe-488e-94a9-2fce2200484c
md"""
# Reading
Fabozzi, Fabozzi, 2021, Bond Markets, Analysis, and Strategies, 10th Edition\
Chapter 7
"""

# ╔═╡ 25448e4a-345b-4bd6-9cdb-a6a280cfd22c
#Set-up packages
begin

	using DataFrames, HTTP, CSV, Dates, Plots, PlutoUI, Printf, LaTeXStrings, HypertextLiteral, XLSX

	gr();
	Plots.GRBackend()

	#Define html elements
	nbsp = html"&nbsp" #non-breaking space
	vspace = html"""

"""
	br = html"
"

	#Sets the width of cells, caps the cell width by 90% of screen width
	#(setting overwritten by cell below)
	# @bind screenWidth @htl("""
	#

	#
	#

	# """)

	# cellWidth= min(1000, screenWidth*0.9)
	# @htl("""
	#
	# """)

	#Sets the width of the cells
	#begin
	#	html"""


text:
type Oneable <: Component
	x
end
type Named <: Component
	x
end

relations=Dict{AbstractString,Vector{Equation}}()

function validated(related::Array,pat::Array)
	tmd=Dict[]
	valid=Dict[]
	for t in 1:length(related)
		ttmd=Dict[]
		pushallunique!(ttmd,matches(related[t],pat[t]))
		if isempty(tmd)
			for ttd in ttmd
				push!(valid,ttd)
			end
		else
			for td in tmd
				for ttd in ttmd
					if !clash(td,ttd)
						push!(valid,combine(td,ttd))
					end
				end
			end
		end
		tmd=valid
		valid=Dict[]
	end
	return tmd
end
function validfilter(c1,c2,mda::Array{Dict})
	filtered=Dict[]
	for md in mda
		if randeval(replace(c2,md))==randeval(c1)
			push!(filtered,md)
		end
	end
	return filtered
end
function matches(c1::Component,c2::Component)
	args1,args2=getargs(c1),getargs(c2)
	return validfilter(c1,c2,validated(args1,args2))
end
function facalloc!(termremains::Array,patremains::Array,psremains::Array,dic::Dict,dica::Array{Dict})
	lps,lpat,lterm=length(psremains),length(patremains),length(termremains)
	ldiff=lterm-lpat
	if lps==2
		for shift in 1:ldiff+1
			tdic=deepcopy(dic)
			tdic[patremains[psremains[1]]]=extract(expression(termremains[1:end-shift]))
			tdic[patremains[psremains[2]]]=extract(expression(termremains[end-shift+1:end]))
			push!(dica,tdic)
		end
	elseif lps==1
		tdic=deepcopy(dic)
		tdic[patremains[psremains[1]]]=extract(expression(termremains[1:end]))
		push!(dica,tdic)
	else
		@assert lps>2
		for shift in 0:ldiff
			tdic=deepcopy(dic)
			tdic[patremains[psremains[end]]]=extract(expression(termremains[end-shift:end]))
			npatremains=deleteat!(deepcopy(patremains),psremains[end])
			pushallunique!(dica,facalloc!(termremains[1:end-1-shift],npatremains,psremains[1:end-1],tdic,dica))
		end

	end
	return dica
end
getcoef(term::Array)=begin;i=indsin(term,Number);isempty(i)?1:sum(term[i]);end
function whenallequal(term,pat,ps)
	mda=Dict[]
	for l in 1:length(ps)
		if isa(term[l],Component)&&isa(pat[ps[l]],Component)
			pushallunique!(mda,matches(term[l],pat[ps[l]]))
		else
			tmd=Dict()
			tmd[pat[ps[l]]]=term[l]
			push!(mda,tmd)
		end
	end
	if length(mda)>1
		md=combine(mda[1],mda[2])
		for i in 3:length(mda)
			md=combine(md,mda[i])
		end
	else
		md=mda[1]
	end
	return md
end
function matches(term::Array{Factor},pat::Array{Factor})
	md=Dict[]
	ps=indsin(pat,Ex)
	lps,lpat,lterm=length(ps),length(pat),length(term)
	if lterm==lpat==lps
		pushallunique!(md,whenallequal(term,pat,ps))
	elseif lterm>lpat==lps
		facalloc!(term,pat,ps,Dict(),md)
	elseif lterm==lpat>lps
		if getcoef(term)==getcoef(pat)
			pushallunique!(md,whenallequal(term[2:end],pat[2:end],ps-1))
		end
	elseif lterm>lpat>lps
		if getcoef(term)==getcoef(pat)
			facalloc!(term[2:end],pat[2:end],ps-1,Dict(),md)
		end
	end
	return md
end
function clash(dic1::Dict,dic2::Dict)
	for key in keys(dic1)
		if haskey(dic2,key)&&dic1[key]!=dic2[key]
			return true
		elseif !haskey(dic2,key)
			for key2 in keys(dic2)
				if dic2[key2]==dic1[key]
					return true
				end
			end
		end
	end
	return false
end
function combine(dic1::Dict,dic2::Dict)
	ndic=Dict()
	for key in keys(dic1)
		ndic[key]=dic1[key]
	end
	for key in keys(dic2)
		ndic[key]=dic2[key]
	end
	return ndic
end
function matches(ex::Expression,pattern::Expression)
	if has(pattern,Pow)
		ex=simplify(ex,Pow)
	end
	md=Dict[]
	apex=terms(sort(componify(ex)))
	if !isa(apex,Array)
		apex=Term[Factor[apex]]
	end
	apat=terms(sort(componify(pattern)))
	if length(apex)==length(apat)==1
		pushallunique!(md,matches(apex[1],apat[1]))
	else
		tmd=validated(apex,apat)
		pushallunique!(md,tmd)
	end
	return md
end
function matches(n::EX,pat::Symbol)
	md=Dict[]
	push!(md,Dict(pat=>n))
	return md
end
matches(::N, ::Expression)=[]
matches(::N, ::Component)=[]
matches(::Number, ::Equation)=[]
function stageoneables(pat)
	ois=expandindices(indsin(pat.lhs,Oneable))
	pat=deepcopy(pat)
	for toi in ois
		pat.lhs[toi]=pat.lhs[toi].x
	end
	pats=Equation[]
	push!(pats,pat)
	if !isempty(ois)
		for oip in permutations(ois)
			np=deepcopy(pat)
			for oi in oip
				np.rhs=replace(np.rhs,Dict(np.lhs[oi]=>1))
				np.lhs[oi]=1
				push!(pats,simplify(np))
			end
		end
	end
	return uniquefilter(pats)
end
function matches(ex::Ex,eq::Equation)
	m=Any[]
	mda=Array{Dict}[]
	if has(eq.lhs,Oneable)
		staged=stageoneables(eq)
		for s in staged
			sm=matches(ex,s.lhs)
			push!(mda,sm)
		end
	else
		staged=[eq]
		push!(mda,matches(ex,eq.lhs))
	end
	for i in 1:length(mda)
		for md in mda[i]
			trh=replace(staged[i].rhs,md)
			push!(m,trh)
		end
	end
	return sort!(m)
end
function quadratic(eq::Equation,xlen::Integer=0,notinx::Array=[])
	eq=simplify(eq)
	if (eq.rhs!=0&&eq.lhs!=0)||(eq.rhs==0&&eq.lhs==0)
		return false
	elseif eq.rhs==0
		termses=terms(eq.lhs)
	else
		termses=terms(eq.rhs)
	end
	connections=Any[]
	matches=Equation[]
	for ti in 1:length(termses)
		matchesfound=Any[]
		for p in permutations(termses[ti])
			for l in 1:length(termses[ti])
				foundmatch=false
				cantbematch=false
				if xlen!=0&&l!=xlen
					cantbematch=true
				end
				if cantbematch
					continue
				end
				for matchi in 1:length(termses)
					if 2l>length(termses[matchi]) || matchi==ti
						continue
					else
						xsq=componify(Expression(Term[p[1:l]])^2)
						if isa(xsq,Expression)
							@assert length(xsq)==1
							xsq=xsq[1]
						end
						for mp in permutations(termses[matchi])
							notx=false
							for notthis in notinx
								if notthis∈mp[1:2l]

									notx=true
									break
								end
							end
							if notx
								continue
							end
							@assert 2l==length(xsq)
							xsqlen=2l
							matchlen=length(termses[matchi])
							numshift=matchlen-xsqlen
							for shif in 0:numshift
								if xsq==mp[1+shif:2l+shif]
									found=Integer[]
									nomatch=false
									for termi in 1:length(termses)
										if termi∈[ti,matchi]#||termi∈found
											continue
										elseif has(termses[termi],p[1:l])
											nomatch=true
											break
										end
									end
									if !nomatch
										a=Factor[]
										for tl in 1:length(mp)
											if !in(tl,collect(1+shif:2l+shif))
												push!(a,mp[tl])
											end
										end
										b=Factor[]
										for tl in l+1:length(p)
											push!(b,p[tl])
										end
										x=p[1:l]
										c=deleteat!(deepcopy(termses),sort([ti,matchi]))
										eq1=Equation(expression(x),(-expression(b)+Sqrt(expression(b)^2-4*expression(a)*expression(c)))/2expression(a))
										if !in(eq1,matches)
											push!(matches,eq1)
										end
										eq2=Equation(expression(x),(-expression(b)-Sqrt(expression(b)^2-4*expression(a)*expression(c)))/2expression(a))
										if !in(eq2,matches)
											push!(matches,eq2)
										end
										foundmatch=true
										break
									end #if
								end #if
							end #for perm xsq
							if foundmatch;break;end
						end #for mp
					end #if
					if foundmatch;break;end
				end #matchi
			end #length(ti)
		end
	end
	return unique(simplify(matches))
end


text:
#using PyCall
using Logging
using Dierckx
using IterTools
#using .snake_forward
#using .curve_utils
using Plots
using MATLAB

include("./snake_forward.jl")
include("./curve_utils.jl")
#Reimplementation of Vedranas method with modifications.

#VEDRANA MODIFIED
function displace(centerline_points, force, radius_func, w; doplot=false)
    L = size(centerline_points,1)
    (outline_xy, normal) = get_outline(centerline_points, radius_func)

    #println(size(force), size(normal), size(w))

    forces = (force.*normal).*w

    displaced_outline = outline_xy .+ forces

    mid = Int64(size(displaced_outline,1)/2)#always even
    upper_points = displaced_outline[3:mid-1,:]
    lower_points = (displaced_outline[mid+3:end-1, :])[end:-1:1,:]

    displaced_centerline = (upper_points + lower_points)./2
    head = (displaced_outline[1,:] + displaced_outline[2,:] + displaced_outline[end,:])./3
    tail = (displaced_outline[mid,:] + displaced_outline[mid+1,:] + displaced_outline[mid+2,:])./3

    displaced_centerline = cat(head', displaced_centerline, dims = 1)
    displaced_centerline = cat(displaced_centerline, tail', dims = 1)

    #enforce π/3 constraint - way to ensure some degree of smoothness
    # cks = get_ck(displaced_centerline)
    #
    # angles_1 = angle.(cks[2:end])
    # angles_2 = angle.(cks[1:end-1])
    #
    # difference = abs.(angles_1-angles_2)
    # signs = sign.(angles_2-angles_1)
    #
    # violation_index = findall(d-> d > π/3, difference)
    # if length(violation_index) > 0
    #
    #     cks[violation_index.+1] .= abs.(cks[violation_index]).*ℯ.^(im*signs[violation_index]*π/3)
    #
    #     displaced_centerline = get_xy(cks)
    # end

    if doplot
        plot!(outline_xy[:,1], outline_xy[:,2], label="original")
        quiver!(outline_xy[:,1], outline_xy[:,2], quiver=(50*forces[:,1], 50*forces[:,2]), label = "forces")
        plot!(displaced_outline[:,1], displaced_outline[:,2], label="displaced")
        plot!(displaced_centerline[:,1], displaced_centerline[:,2], label="displaced")
    end

    return displaced_centerline
end

function move_points(residual,curves,angles,N,centerline_points,r,w; doplot=false)
    (x_length, y_length) = size(residual)
    vals = zeros(Float64, N)
    if y_length > 1
        F = Spline2D(collect(1:1.0:x_length), collect(1:1.0:y_length), residual, kx=1, ky=1);
        vals = zeros(Float64, N)
        for i = 1:length(angles)
            interp = F(curves[:,i], repeat([i], N))
            vals += interp
        end
    else
        F = Spline1D(collect(1:1.0:x_length), residual[:,1], k=1);
        interp = F(curves[:,1])
        vals = interp
    end

    force = vals*(1/length(angles))

    centerline_points = displace(centerline_points, force, r, w, doplot=doplot)
    return centerline_points
end

function to_pixel_coordinates(current, angles, bins)
    N = size(current,1);
    vertex_coordinates = zeros(Float64,N,length(angles));
    a = (length(bins)-1)/(bins[end]-bins[1]); # slope
    b = 1-a*bins[1]; # intercept
    for k = 1:length(angles)
        angle = angles[k]
        projection = [cos(angle) sin(angle)]';
        #expressing vertex coordinates as coordinates in sinogram (pixel coordinates, not spatial)
        vertex_coordinates[:,k] = (current*projection)*a.+b;
    end
    return vertex_coordinates
end

function evolve_curve(sinogram_target, centerline_points, r, angles, bins, max_iter, w, smoothness, degree::Int64; doplot=false)
    (current, normal) = get_outline(centerline_points, r)
    current_sinogram = parallel_forward(current,angles,bins)

    curves = to_pixel_coordinates(current, angles, bins);
    mu = sum(sinogram_target[:].*current_sinogram[:])/sum(current_sinogram[:].^2)
    residual = sinogram_target - mu*current_sinogram
    N = size(current,1)
    centerline_start = centerline_points[1,:]
    for iter  = 1:max_iter
        centerline_points = move_points(residual,curves,angles,N,centerline_points,r, w, doplot=doplot)

        L = size(centerline_points,1)
        #HACK
        cp = eliminate_loopy_stuff(centerline_points, 2*r(0.0))
        #HACK
        if size(cp,1) > degree
            centerline_points = cp
        else
            @warn "Too loopy"
        end
        t = curve_lengths(centerline_points)

        spl = ParametricSpline(t,centerline_points',k=degree, s=0.0)
        #HACK
        if smoothness > 0.0
            try
                spl = ParametricSpline(t,centerline_points',k=degree, s=smoothness)
            catch e
                @warn e
            end
        end

        tspl = range(0, t[end], length=L)
        centerline_points = collect(spl(tspl)')
        centerline_points[1,:] = centerline_start
        (current, normal) = get_outline(centerline_points, r)
        current_sinogram = parallel_forward(current,angles,bins);
        curves = to_pixel_coordinates(current, angles, bins);
        mu = sum(sinogram_target[:].*current_sinogram[:])/sum(current_sinogram[:].^2)
        residual = sinogram_target - mu*current_sinogram;
    end
    return centerline_points
end

function recon2d_tail(centerline_points::AbstractArray{T,2}, r, angles::Array{T},bins::Array{T},sinogram_target::Array{T,2}, max_iter::Int, smoothness::T, w::Array{T}, degree::Int64; doplot=false) where T<:AbstractFloat
    current = evolve_curve(sinogram_target, centerline_points, r, angles, bins, max_iter, w, smoothness, degree, doplot=doplot)
    return current
end

# #VEDRANA ORIGINAL
function regularization_matrix(N,alpha,beta)
    cwd = @__DIR__
    println(cwd)
    mat"[$B, $A] = regularization_matrix($N,$alpha,$beta);"
    return B
end
#
function remove_crossings(curve)
    mat"$curve = remove_crossings($curve);"
    return curve
end
#
function distribute_points(curve)
    curve = cat(dims = 1, curve, curve[1,:]'); # closing the curve
    N = size(curve,1); # number of points [+ 1, due to closing]
    dist = sqrt.(sum(diff(curve, dims=1).^2, dims=2))[:,1]; # edge segment lengths
    t = prepend!(cumsum(dist, dims=1)[:,1],0.0) # total curve length

    tq = range(0,t[end],length=N); # equidistant positions
    curve_new_1 = Spline1D(t,curve[:,1], k=1).(tq); # distributed x
    curve_new_2 = Spline1D(t,curve[:,2], k=1).(tq); # distributed y
    curve_new = hcat(curve_new_1,curve_new_2); # opening the curve again
    return curve_new[1:end-1,:]
end

function move_points(residual,curves,angles,N,centerline_points,r,w; doplot=false)
    (x_length, y_length) = size(residual)
    vals = zeros(Float64, N)
    if y_length > 1
        F = Spline2D(collect(1:1.0:x_length), collect(1:1.0:y_length), residual, kx=1, ky=1);
        vals = zeros(Float64, N)
        for i = 1:length(angles)
            interp = F(curves[:,i], repeat([i], N))
            vals += interp
        end
    else
        F = Spline1D(collect(1:1.0:x_length), residual[:,1], k=1);
        interp = F(curves[:,1])
        vals = interp
    end

    force = vals*(1/length(angles))

    centerline_points = displace(centerline_points, force, r, w, doplot=doplot)
    return centerline_points
end

function move_points_original(residual,curves,angles,N,current,B,w)
    (x_length, y_length) = size(residual)
    #F = Spline2D(collect(1:1.0:x_length), collect(1:1.0:y_length), residual, kx=1, ky=1);
    vals = zeros(Float64, N)
    if y_length > 1
        F = Spline2D(collect(1:1.0:x_length), collect(1:1.0:y_length), residual, kx=1, ky=1);
        vals = zeros(Float64, N)
        for i = 1:length(angles)
            interp = F(curves[:,i], repeat([i], N))
            vals += interp
        end
    else
        F = Spline1D(collect(1:1.0:x_length), residual[:,1], k=1);
        interp = F(curves[:,1])
        vals = interp
    end
    force = vals*(1/length(angles))
    normals = snake_normals(current)
    vectors = force.*normals
    current = current + w*vectors;
    current = distribute_points(remove_crossings(B*current))
    return current
end

function evolve_curve_original(sinogram_target, current, angles, bins, B, max_iter, w)
    current_sinogram = parallel_forward(current,angles,bins)

    curves = to_pixel_coordinates(current, angles, bins);

    mu = sum(sinogram_target[:].*current_sinogram[:])/sum(current_sinogram[:].^2)
    residual = sinogram_target - mu*current_sinogram
    N = size(current,1)

    for iter  = 1:max_iter
        current = move_points_original(residual,curves,angles,N,current,B,w);
        current_sinogram = parallel_forward(current,angles,bins);
        curves = to_pixel_coordinates(current, angles, bins);
        mu = sum(sinogram_target[:].*current_sinogram[:])/sum(current_sinogram[:].^2)
        residual = sinogram_target - mu*current_sinogram;
    end
    return current
end

function vedrana(current::Array{T,2},angles::Array{T},bins::Array{T},sinogram_target::Array{T,2}, max_iter::Int, alpha::T, beta::T, w::T) where T<:AbstractFloat
    N = size(current,1)
    B = regularization_matrix(N,alpha,beta)
    current = evolve_curve_original(sinogram_target, current, angles, bins, B, max_iter, w)
    return current
end
#

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
"""
Module for downloading a user's photos from Facebook.

Copyright Peter Brescia 2016.
"""
module FacebookPhotoFetcher

export fetch_photos

using Requests
using SQLite
using JSON

"""
Entry point for this tool.

function fetch_photos(access_token::String, save_directory::String)
"""
function fetch_photos(access_token::String, save_directory::String)
    output = set_up_output(save_directory)
    url = "https://graph.facebook.com/v2.5/me/photos?access_token=$access_token&limit=50"

    while true
        url = fetch_photos_in_batches(url, access_token, output)
        if url === false
            break
        end
    end
    println("Photo fetching complete")
end

function fetch_photos_in_batches(
    url::String,
    access_token::String,
    output::Dict,
)
    res = Requests.get(url)
    parsed_response = Requests.json(res)

    data = get(parsed_response, "data", false)
    if data === false
        println(parsed_response)
    end

    for photo in parsed_response["data"]
        fetch_image(photo["id"], access_token, output)
    end
    return get(parsed_response["paging"], "next", false) === false ? false : parsed_response["paging"]["next"]
end

function fetch_image(
    photo_id::String,
    access_token::String,
    output::Dict,
)
    url = "https://graph.facebook.com/v2.8/$photo_id?access_token=$access_token&fields=backdated_time_granularity,album,images,name,created_time,picture,link,event,from,name_tags,place,updated_time,width,page_story_id,can_delete,backdated_time,can_tag,height,icon,id,comments,likes,tags"

    res = Requests.get(url)
    parsed_response = Requests.json(res)

    SQLite.query(
        output["db"],
        "INSERT INTO images (facebook_id, metadata) VALUES (?, ?)";
        values=[
            photo_id,
            JSON.json(parsed_response),
        ]
    )

    img_source = parsed_response["images"][1]["source"]
    img = Requests.get(img_source)
    filename_fragments = split(img_source, "/")
    filename_last_fragment = filename_fragments[length(filename_fragments)]
    filename_subfragments = split(filename_last_fragment, "?")
    filename = filename_subfragments[1]

    Requests.save(img, output["output_directory"]*"/$filename")
end

function set_up_output(save_directory)
    if !isdir(save_directory)
        mkdir(save_directory)
    end
    timestamp = Dates.now()
    output_directory = "$save_directory/$timestamp"
    mkdir(output_directory)

    db = SQLite.DB("$output_directory/output.sqlite")
    SQLite.query(
        db,
        "CREATE TABLE images (
            id INTEGER not null PRIMARY KEY autoincrement,
            facebook_id INTEGER not null,
            metadata TEXT not null
        )",
    )
    SQLite.query(
        db,
        "CREATE INDEX images_fid_idx ON images (facebook_id)",
    )
    Dict(
        "db"                => db,
        "output_directory"  => output_directory,
    )
end

end

max_stars_repo_name: peterbrescia/FacebookPhotoFetcher

max_issues_repo_name: peterbrescia/FacebookPhotoFetcher

max_forks_repo_name: peterbrescia/FacebookPhotoFetcher


text:
function prompt()
    println("Please enter length of short side A if possible, if not leave blank")
    a = readline()
    println("Please enter length of short side B if possible, if not leave blank")
    b = readline()
    println("Please enter length of Hypotinous C if possible, if not leave blank")
    c = readline()

    all_val = a,b,c

    all_val = [all_val[i] == "" ?  missing : all_val[i] for i in 1:length(all_val)]

    if sum(x -> ismissing(x), all_val) == 1
        nothing
    else
        println("Error: too many or too few missing values.")
    end

    all_v = []
    for i in 1:size(all_val,1)
        if  ismissing(all_val[i])
            push!(all_v, missing)
        else
            push!(all_v, parse(Float64,all_val[i]) )
        end
    end

    return all_v

end

all_v = prompt()

function Pythag(inp::Array{Any,1})
    if ismissing(inp[3])
        result = sqrt(inp[1]^2 + inp[2]^2)
    else
        vals = collect(skipmissing(inp))
        result = sqrt( vals[2]^2 - vals[1]^2 )
    end
    return result
end

println("The length of the missing side is: $(Pythag(all_v))")


text:
function Lagrange(x,xP,iP)
f=1;
for i=1:size(xP,1)
   if i≠iP
      f=f*(x-xP[i])/(xP[iP]-xP[i]);
   end
end
return f
end


text:
## Instrumental variables problems

dt = readcsv("iv.csv");

y = dt[:,1];
x = dt[:,2:4];
z = dt[:,5:end];

f(theta) = z.*(y-x*theta);

gmm2s = GMMEstimator(f, [.0, .0, .0])
estimate!(gmm2s)

Dsn(θ)    = -z'x;
Dws(θ, p) = -z'*Diagonal(p)*x;
Dsl(θ, λ) = -x.*(z*λ);
Hwsl(θ, p, λ) = zeros(3,3);

mds = MDEstimator(f, [.0, .0, .0], grad = (Dsn, Dws, Dsl, Hwsl));

## GMM using linear algebra

## Stata Coef
## -.0050209   .0708816   .0593661

## Stata vcov
## symmetric e(V)[3,3]
##             v2          v3          v4
##    .00478911
##   -.00272551   .01026473
##    .00163441  -.00201291   .00983372

##   x1 |  -.0050209   .0692034    -0.07   0.942     -.140657    .1306153
##   x2 |   .0708816    .101315     0.70   0.484    -.1276922    .2694554
##   x3 |   .0593661   .0991651     0.60   0.549     -.134994    .2537261

## Hansen J statistic (overidentification test of all instruments):         1.191
##                                                     Chi-sq(4) P-val =    0.8795

## srand(1)
## using MomentBasedEstimators
## using ModelsGenerators
## y,x,z = ModelsGenerators.sim_iv_d01(n = 10000, CP = 2000)
## xb    = randn(10000,2)
## x     = [x xb]
## z     = [z xb]

## f(theta) = z.*(y-x*theta);

## gmm = MomentBasedEstimators.GMMEstimator(f, [.1, 0, 0], initialW = z'z);
## md = MomentBasedEstimators.MDEstimator(f, [.1, 0, 0]);

## MomentBasedEstimators.initialize!(gmm);
## MomentBasedEstimators.initialize!(md);

## estimate!(gmm);
## estimate!(md);

## MomentBasedEstimators.setparLB!(md, coef(md))
## MomentBasedEstimators.setparUB!(md, coef(md))
## estimate!(md);

## MomentBasedEstimators.setparLB!(md, [0.,0,0])
## MomentBasedEstimators.setparUB!(md, [0.,0,0])

## ## Methods

## # 1. Coefficient
## coef(gmm)
## coef(md)

## ## 2. Variances
## vcov(gmm)
## vcov(md)

## ## 3. Objectiv value

## srand(1)
## y,x,z=randiv(100)
## xb = randn(100,2)
## x = [x xb]
## z = [z xb]

## z=sparse(z)

## f(theta) = broadcast(*, z, (y-x*theta))

## g = MomentBasedEstimators.GMMEstimator(f, [.1, 0, 0])
## m = MomentBasedEstimators.MDEstimator(f, [.1, 0, 0])

## MomentBasedEstimators.initialize!(g);
## MomentBasedEstimators.initialize!(m);

## @time estimate!(g);
## @time estimate!(m);

## srand(1)
## y,x,z=randiv(100)
## xb = randn(100,2)
## x = [x xb]
## z = [z xb]

## f(theta) = z.*(y-x*theta);

## g = MomentBasedEstimators.GMMEstimator(f, [.1, 0, 0])
## m = MomentBasedEstimators.MDEstimator(f, [.1, 0, 0])

## MomentBasedEstimators.initialize!(g)
## MomentBasedEstimators.initialize!(m)

## @time estimate!(g);
## @time estimate!(m);

## Before
## elapsed time: 0.129503039 seconds (69486712 bytes allocated, 40.36% gc time)


text:
### Instruction:
# Run after mainMulti.jl is finished.
# sim824.jld needs to be in the same output folder.
# For definition of inputs, see postMonteCarlo function in Post.jl
push!(LOAD_PATH,pwd())

using Revise
using TypeDef, Helper, Tools, Compute, Solve, FOCs, UserFileIO, MonteCarlo
using Post
using JLD, Printf

### Post-simulation in wholesale-retail markets
dim,eqb,tab = postMonteCarlo(15,10,false,false,0,true)

### Post-simulation under vertical integration
dim,eqb,tab = postMonteCarlo(15,10,true,false,0)

### Demand elasticity for Model 0 (Table 4 Column 3)
dim,eqb,tab = postMonteCarlo( 0,20,false,false,0,false)

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
include("dependencies_for_runtests.jl")
using Oceananigans.Models.HydrostaticFreeSurfaceModels
import Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitFreeSurface
import Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, SplitExplicitAuxiliary, SplitExplicitSettings
import Oceananigans.Models.HydrostaticFreeSurfaceModels: barotropic_mode!, barotropic_split_explicit_corrector!, set_average_to_zero!

@testset "Barotropic Kernels" begin

    for arch in archs
        FT = Float64
        topology = (Periodic, Periodic, Bounded)
        Nx = Ny = 16 * 8
        Nz = 32
        Nx = 128
        Ny = 64
        Lx = Ly = Lz = 2π
        grid = RectilinearGrid(arch, topology = topology, size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0))

        tmp = SplitExplicitFreeSurface()
        sefs = SplitExplicitState(grid)
        sefs = SplitExplicitAuxiliary(grid)
        sefs = SplitExplicitFreeSurface(grid)

        state = sefs.state
        auxiliary = sefs.auxiliary
        U, V, η̅, U̅, V̅ = state.U, state.V, state.η̅, state.U̅, state.V̅
        Gᵁ, Gⱽ = auxiliary.Gᵁ, auxiliary.Gⱽ
        Hᶠᶜ, Hᶜᶠ = sefs.auxiliary.Hᶠᶜ, sefs.auxiliary.Hᶜᶠ

        u = Field{Face,Center,Center}(grid)
        v = Field{Center,Face,Center}(grid)

        @testset "Average to zero" begin
            # set equal to something else
            η̅ .= U̅ .= V̅ .= 1.0
            # now set equal to zero
            set_average_to_zero!(sefs.state)
            # don't forget the ghost points
            fill_halo_regions!(η̅, arch)
            fill_halo_regions!(U̅, arch)
            fill_halo_regions!(V̅, arch)
            # check
            @test all(Array(η̅.data.parent) .== 0.0)
            @test all(Array(U̅.data.parent .== 0.0))
            @test all(Array(V̅.data.parent .== 0.0))
        end

        @testset "Inexact integration" begin
            # Test 2: Check that vertical integrals work on the CPU(). The following should be "inexact"
            Δz = zeros(Nz)
            Δz .= grid.Δzᵃᵃᶠ

            set_u_check(x, y, z) = cos((π / 2) * z / Lz)
            set_U_check(x, y) = (sin(0) - (-2 * Lz / (π)))
            set!(u, set_u_check)
            exact_U = similar(U)
            set!(exact_U, set_U_check)
            barotropic_mode!(U, V, grid, u, v)
            tolerance = 1e-3
            @test all((Array(interior(U) .- interior(exact_U))) .< tolerance)

            set_v_check(x, y, z) = sin(x * y) * cos((π / 2) * z / Lz)
            set_V_check(x, y) = sin(x * y) * (sin(0) - (-2 * Lz / (π)))
            set!(v, set_v_check)
            exact_V = similar(V)
            set!(exact_V, set_V_check)
            barotropic_mode!(U, V, grid, u, v)
            @test all((Array(interior(V) .- interior(exact_V))) .< tolerance)
        end

        @testset "Vertical Integral " begin
            Δz = zeros(Nz)
            Δz .= grid.Δzᵃᵃᶜ

            u .= 0.0
            U .= 1.0
            barotropic_mode!(U, V, grid, u, v)
            @test all(Array(U.data.parent) .== 0.0)

            u .= 1.0
            U .= 1.0
            barotropic_mode!(U, V, grid, u, v)
            @test all(Array(interior(U)) .≈ Lz)

            set_u_check(x, y, z) = sin(x)
            set_U_check(x, y) = sin(x) * Lz
            set!(u, set_u_check)
            exact_U = similar(U)
            set!(exact_U, set_U_check)
            barotropic_mode!(U, V, grid, u, v)
            @test all(Array(interior(U)) .≈ Array(interior(exact_U)))

            set_v_check(x, y, z) = sin(x) * z * cos(y)
            set_V_check(x, y) = -sin(x) * Lz^2 / 2.0 * cos(y)
            set!(v, set_v_check)
            exact_V = similar(V)
            set!(exact_V, set_V_check)
            barotropic_mode!(U, V, grid, u, v)
            @test all(Array(interior(V)) .≈ Array(interior(exact_V)))
        end

        @testset "Barotropic Correction" begin
            # Test 4: Test Barotropic Correction
            FT = Float64
            topology = (Periodic, Periodic, Bounded)
            Nx = Ny = 16 * 8
            Nz = 32
            Nx = 128
            Ny = 64
            Lx = Ly = Lz = 2π
            grid = RectilinearGrid(arch, topology = topology, size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0))

            sefs = SplitExplicitFreeSurface(grid)

            state = sefs.state
            auxiliary = sefs.auxiliary
            U, V, η̅, U̅, V̅ = state.U, state.V, state.η̅, state.U̅, state.V̅
            Gᵁ, Gⱽ = auxiliary.Gᵁ, auxiliary.Gⱽ
            Hᶠᶜ, Hᶜᶠ = sefs.auxiliary.Hᶠᶜ, sefs.auxiliary.Hᶜᶠ

            u = Field{Face,Center,Center}(grid)
            v = Field{Center,Face,Center}(grid)
            u_corrected = similar(u)
            v_corrected = similar(v)

            set_u(x, y, z) = z + Lz / 2 + sin(x)
            set_U̅(x, y) = cos(x) * Lz
            set_u_corrected(x, y, z) = z + Lz / 2 + cos(x)
            set!(u, set_u)
            set!(U̅, set_U̅)
            set!(u_corrected, set_u_corrected)

            set_v(x, y, z) = (z + Lz / 2) * sin(y) + sin(x)
            set_V̅(x, y) = (cos(x) + x) * Lz
            set_v_corrected(x, y, z) = (z + Lz / 2) * sin(y) + cos(x) + x
            set!(v, set_v)
            set!(V̅, set_V̅)
            set!(v_corrected, set_v_corrected)

            sefs.auxiliary.Hᶠᶜ .= Lz
            sefs.auxiliary.Hᶜᶠ .= Lz

            Δz = zeros(Nz)
            Δz .= grid.Δzᵃᵃᶜ

            barotropic_split_explicit_corrector!(u, v, sefs, grid)
            @test all(Array((interior(u) .- interior(u_corrected))) .< 1e-14)
            @test all(Array((interior(v) .- interior(v_corrected))) .< 1e-14)
        end
    end # end of architecture loop
end # end of testset


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using Documenter, FEMTruss

# copy README.md to index.md to avoid DRY
cp("../README.md", "src/index.md"; remove_destination=true)

makedocs(modules=[FEMTruss],
         format = :html,
         checkdocs = :all,
         sitename = "FEMTruss.jl",
         pages = [
                  "Introduction" => "index.md",
                  "API documentation" => "api.md"
                 ])


text:
using KMarkovGenSeq

function main()
    cnt, pred_cls, des = query_freq("data/reads.fa", "data/genomes", 9, nothing, false)
    len = length(des)
    println("Total query sequences: $cnt")
    counter = fill(0, len)
    for cls in pred_cls
        counter[cls] += 1
    end
    println()
    for threshold in [1, 5, 10, 50]
        println("Number of groups containing at least $threshold query sequences: ", count(>=(threshold), counter))
    end
    for i in eachindex(counter)
        if counter[i] >= 10
            println(des[i], "contains\t", counter[i], "\tquery sequences")
        end
    end
    println()
    @show(counter)
end

main()

"""
Total query sequences: 1876

Number of groups containing at least 1 query sequences: 10
Number of groups containing at least 5 query sequences: 10
Number of groups containing at least 10 query sequences: 10
Number of groups containing at least 50 query sequences: 8

Baumannia cicadellinicola str. Hc (Homalodisca coagulata), complete genomecontains      298     query sequences
Psychromonas ingrahamii 37 chromosome, complete genomecontains  22      query sequences
Sphingomonas wittichii RW1 chromosome, complete genomecontains  81      query sequences
Roseiflexus castenholzii DSM 13941 chromosome, complete genomecontains  375     query sequences
Hydrogenobaculum sp. Y04AAS1 chromosome, complete genomecontains        104     query sequences
Alteromonas macleodii str. 'Deep ecotype' chromosome, complete genomecontains   279     query sequences
Denitrovibrio acetiphilus DSM 12809 chromosome, complete genomecontains 111     query sequences
Frankia symbiont of Datisca glomerata chromosome, complete genomecontains       23      query sequences
Candidatus Midichloria mitochondrii IricVA chromosome, complete genomecontains  348     query sequences
Corynebacterium variabile DSM 44702 chromosome, complete genomecontains 235     query sequences

counter = [298, 22, 81, 375, 104, 279, 111, 23, 348, 235]
"""


text:
# ------------------------------------------------------------------------------
# finite element operator gallery
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# user utility constructors
# ------------------------------------------------------------------------------

"""
```julia
GalleryOperator(name, numbernodes1d, numberquadraturepoints1d, mesh)
```

Finite element operator from a gallery of options

# Arguments:
- `name`:                      string containing name of operator
- `numbernodes1d`:             polynomial order of TensorH1LagrangeBasis
- `numberquadraturepoints1d`:  number of quadrature points in one dimension for basis
- `mesh`:                      mesh for operator

# Returns:
- Finite element operator object

# Mass matrix example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# verify
println(mass)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
```

# Diffusion operator example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryOperator("diffusion", 4, 4, mesh);

# verify
println(diffusion)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
```
"""
function GalleryOperator(
    name::String,
    numbernodes1d::Int,
    numberquadraturepoints1d::Int,
    mesh::Mesh,
)
    if haskey(operatorgallery, name)
        basis = TensorH1LagrangeBasis(
            numbernodes1d,
            numberquadraturepoints1d,
            1,
            mesh.dimension,
        )
        return operatorgallery[name](basis, mesh)
    else
        throw(ArgumentError("operator name not found")) # COV_EXCL_LINE
    end
end

"""
```julia
GalleryVectorOperator(name, numbernodes1d, numberquadraturepoints1d, numberelements1d, mesh)
```

Finite element operator from a gallery of options

# Arguments:
- `name`:                      string containing name of operator
- `numbernodes1d`:             polynomial order of TensorH1LagrangeBasis
- `numberquadraturepoints1d`:  number of quadrature points in one dimension for basis
- `numbercomponents`:          number of components
- `mesh`:                      mesh for operator

# Returns:
- Finite element operator object

# Mass matrix example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryVectorOperator("mass", 4, 4, 3, mesh);

# verify
println(mass)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    interpolation
```

# Diffusion operator example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryVectorOperator("diffusion", 4, 4, 3, mesh);

# verify
println(diffusion)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    gradient
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 3
    dimension: 2
  evaluation mode:
    gradient
```
"""
function GalleryVectorOperator(
    name::String,
    numbernodes1d::Int,
    numberquadraturepoints1d::Int,
    numbercomponents::Int,
    mesh::Mesh,
)
    if haskey(operatorgallery, name)
        basis = TensorH1LagrangeBasis(
            numbernodes1d,
            numberquadraturepoints1d,
            numbercomponents,
            mesh.dimension,
        )
        return operatorgallery[name](basis, mesh)
    else
        throw(ArgumentError("operator name not found")) # COV_EXCL_LINE
    end
end

"""
```julia
GalleryMacroElementOperator(name, numbernodes1d, numberquadraturepoints1d, numberelements1d, mesh)
```

Finite element operator from a gallery of options

# Arguments:
- `name`:                      string containing name of operator
- `numbernodes1d`:             polynomial order of TensorH1LagrangeBasis
- `numberquadraturepoints1d`:  number of quadrature points in one dimension for basis
- `numberelements1d`:          number of elements in macro-element
- `mesh`:                      mesh for operator

# Returns:
- Finite element operator object

# Mass matrix example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryMacroElementOperator("mass", 4, 4, 2, mesh);

# verify
println(mass)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    interpolation
```

# Diffusion operator example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
diffusion = GalleryMacroElementOperator("diffusion", 4, 4, 2, mesh);

# verify
println(diffusion)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    gradient
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  macro-element tensor product basis:
    numbernodes1d: 7
    numberquadraturepoints1d: 8
    numbercomponents: 1
    numberelements1d: 2
    dimension: 2
  evaluation mode:
    gradient
```
"""
function GalleryMacroElementOperator(
    name::String,
    numbernodes1d::Int,
    numberquadraturepoints1d::Int,
    numberelements1d::Int,
    mesh::Mesh,
)
    if haskey(operatorgallery, name)
        basis = TensorH1LagrangeMacroBasis(
            numbernodes1d,
            numberquadraturepoints1d,
            1,
            mesh.dimension,
            numberelements1d,
        )
        return operatorgallery[name](basis, mesh)
    else
        throw(ArgumentError("operator name not found")) # COV_EXCL_LINE
    end
end

# ------------------------------------------------------------------------------
# Operator Gallery
# ------------------------------------------------------------------------------

"""
```julia
massoperator(basis, mesh)
```

Convenience constructor for mass operator

# Weak form:
- ``\\int v u``

# Arguments:
- `basis`: basis to use of all operator fields
- `mesh`:  mesh for operator

# Returns:
- Mass matrix operator with basis on mesh

# Example:
```jldoctest
# mass operator
mesh = Mesh2D(1.0, 1.0);
basis = TensorH1LagrangeBasis(3, 4, 1, mesh.dimension)
mass = LFAToolkit.massoperator(basis, mesh);

# verify
println(mass)

# output
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
```
"""
function massoperator(basis::AbstractBasis, mesh::Mesh)
    # setup
    function massweakform(u::Array{Float64}, w::Array{Float64})
        v = u * w[1]
        return [v]
    end

    # fields
    inputs = [
        OperatorField(basis, [EvaluationMode.interpolation]),
        OperatorField(basis, [EvaluationMode.quadratureweights]),
    ]
    outputs = [OperatorField(basis, [EvaluationMode.interpolation])]

    # operator
    mass = Operator(massweakform, mesh, inputs, outputs)
    return mass
end

"""
```julia
diffusionoperator(basis, mesh)
```

Convenience constructor for diffusion operator

# Weak form:
- ``\\int \\nabla v \\nabla u``

# Arguments:
- `basis`:  basis for all operator fields to use
- `mesh`:   mesh for operator

# Returns:
- Diffusion operator with basis on mesh

# Example:
```jldoctest
# diffusion operator
mesh = Mesh2D(1.0, 1.0);
basis = TensorH1LagrangeBasis(3, 4, 1, mesh.dimension)
diffusion = LFAToolkit.diffusionoperator(basis, mesh);

# verify
println(diffusion)

# output

finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 3
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
```
"""
function diffusionoperator(basis::AbstractBasis, mesh::Mesh)
    # setup
    function diffusionweakform(du::Array{Float64}, w::Array{Float64})
        dv = du * w[1]
        return [dv]
    end

    # fields
    inputs = [
        OperatorField(basis, [EvaluationMode.gradient]),
        OperatorField(basis, [EvaluationMode.quadratureweights]),
    ]
    outputs = [OperatorField(basis, [EvaluationMode.gradient])]

    # operator
    diffusion = Operator(diffusionweakform, mesh, inputs, outputs)
    return diffusion
end

# ------------------------------------------------------------------------------
# operator gallery dictionary
# ------------------------------------------------------------------------------

operatorgallery = Dict("mass" => massoperator, "diffusion" => diffusionoperator)

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


text:
#=
Ray, T. (2003). Golinski’s speed reducer problem revisited. AIAA Journal, 41(3), 556–558. https://doi.org/10.2514/2.1984
=#

function speed_reducer(solver = OCTHaGOn.SOLVER_SILENT)
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, x[1:8])
    @objective(m, Min, x[8])
    gm = GlobalModel(model = m, name = "speed_reducer")
    lbs = Dict(x .=> [2.6, 0.7, 17, 7.3, 7.3, 2.9, 5, 0])
    ubs = Dict(x .=> [3.6, 0.8, 28, 8.3, 8.3, 3.9, 5.5, 5000])
    JuMP.set_integer(x[3])
    bound!(gm, Dict(var => [lbs[var], ubs[var]] for var in gm.vars))

    # Objective
    # add_nonlinear_constraint(gm, :(x -> x[8] - (0.7854*x[1]*x[2]^2*(10. / 3. *x[3]^2 + 14.9334*x[3] - 43.0934) -
    #          1.5079*x[1]*(x[6]^2 + x[7]^2) + 7.477*(x[6]^3 + x[7]^3) +
    #          0.7854*(x[4]*x[6]^2 + x[5]*x[7]^2))), vars = x[1:8])
    add_nonlinear_constraint(gm, :(x -> 0.7854*x[1]*x[2]^2*(3.3333 *x[3]^2 + 14.9334*x[3] - 43.0934) -
    1.5079*x[1]*(x[6]^2 + x[7]^2) + 7.477*(x[6]^3 + x[7]^3) +
    0.7854*(x[4]*x[6]^2 + x[5]*x[7]^2)), vars = x[1:7], dependent_var = x[8])

    # Constraints
    add_nonlinear_constraint(gm, :(x -> -27 + x[1] * x[2]^2 * x[3]), vars = x[1:3])

    add_nonlinear_constraint(gm, :(x -> -397.5 + x[1] * x[2]^2 * x[3]^2), vars = x[1:3])

    add_nonlinear_constraint(gm, :(x -> -1.93 + x[2] * x[6]^4 * x[3] / x[4]^3),  vars = [x[2], x[3], x[4], x[6]])

    add_nonlinear_constraint(gm, :(x -> -1.93 + x[2] * x[7]^4 * x[3] / x[5]^3), vars = [x[2], x[3], x[5], x[7]])

    add_nonlinear_constraint(gm, :(x -> 110.0*x[6]^3 - (((745*x[4]/(x[2]*x[3]))^2) +
    16.9*10^6)^0.5), vars = [x[2], x[3], x[4], x[6]])

    add_nonlinear_constraint(gm, :(x -> 85.0*x[7]^3 - (((745*x[5]/(x[2]*x[3]))^2) +
                                            157.5*10^6)^0.5), vars = [x[2], x[3], x[5], x[7]])

    add_nonlinear_constraint(gm, :(x -> 40 - x[2]*x[3]), vars = [x[2], x[3]])

    add_nonlinear_or_compatible(gm, :(x -> x[1] - 5*x[2]), vars = [x[1], x[2]])

    add_nonlinear_or_compatible(gm, :(x -> 12*x[2] - x[1]), vars = [x[1], x[2]])

    add_nonlinear_or_compatible(gm, :(x -> x[4] - 1.5*x[6] - 1.9), vars = [x[4], x[6]])

    add_nonlinear_or_compatible(gm, :(x -> x[5] - 1.1*x[7] - 1.9), vars = [x[5], x[7]])

    set_param(gm, :ignore_accuracy, true)

    return gm
end


text:
export CSSWidgetWrapper

"""
    CSSWidgetWrapper(widget, style::Dict{String,String})

Modify the style properties of the first HTML node for the given widget.

# Examples
Generate a text field with a red background and white text.
```
@bind text CSSWidgetWrapper(TextField(),Dict("backgroundColor"=>"red","color"=>"white"))
```
Generate a text field with 1px wide solid red border.
```
@bind text CSSWidgetWrapper(TextField(),Dict("border"=>"1px solid red"))
```
"""
struct CSSWidgetWrapper
    widget
    style::Dict{String,String}
end

function show(io::IO, m::MIME"text/html", w::CSSWidgetWrapper)
    style = join(["widget.style.$k = '$v';" for (k,v) in w.style],"\n")
    script = """
        <script>
        var widget = currentScript.parentElement.firstElementChild
        $style
        </script>
    """
    show(io,m,w.widget)
    print(io,script)
end

Base.get(w::CSSWidgetWrapper) = get(w.widget)


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, alignmode = Outside(30))

    la = maingl[1, 1] = LAxis(scene)
    la.attributes.yautolimitmargin = (0f0, 0.05f0)

    poly!(la, BBox(0, 1, 4, 0), color=:blue)
    poly!(la, BBox(1, 2, 7, 0), color=:red)
    poly!(la, BBox(2, 3, 1, 0), color=:green)

    la.attributes.xticks = ManualTicks([0.5, 1.5, 2.5], ["blue", "red", "green"])
    la.attributes.xlabel = "Color"
    la.attributes.ylabel = "Value"
end


text:
using XLSX, LinearAlgebra

data = XLSX.readxlsx("transfer_m.xlsx")

data = data["Sheet6"]
num = length(data["F"])
R_matrix = zeros(2, num)
I_matrix = zeros(2, num)

R_matrix[1, :] = data["F"]
R_matrix[2, :] = data["H"]

I_matrix[1, :] = data["B"]
I_matrix[2, :] = data["D"]

M = I_matrix * pinv(R_matrix)

@show(M)