text:
#------------------------------
#   Solve for the equilibrium interest rate in the Hugget economy with a general diffusion process.
#
#  First, solve the HJB equation to get the A ("transition") matrix, then
#   solve the discretized KF equations. Then iterate over interest rates to
#   get the asset supply function. Update interest rate until convergence.
#
#   This code uses an implicit method for the solution: use
#   an 'upwind' scheme consisting of forwards and backwards difference
#   approximation to the derivative of the value function.
#
#   Assume: CRRA utility function u(c) = c^(1-s)/(1-s)
#   Note that r < rho
#
#   Written by Ben Moll et al (2015). Optimized for speed by SeHyoun Ahn. Translation into Julia by James Graham
#------------------------------

type Huggett_model
  gamma::Float64            # risk aversion coefficient
  rho::Float64              # discount rate.
  ou_var::Float64           # Ornstein Uhlenbeck variance
  zmean::Float64            # Mean of log-normal distriution for the Ornstein Uhlenbeck process
  ou_corr::Float64          # Ornstein Uhlenbeck correlation
  theta::Float64            # Ornstein Uhlenbeck param
  sig2::Float64            # Ornstein Uhlenbeck variance
  I::Int64                  # Number of grids for assets
  J::Int64                  # Number of grids for income
  Ir::Int64                 # Number of interest rate iterations
  amin::Float64             # Min assets = borrowing constraint
  amax::Float64             # Max assets
  zmin::Float64             # Min income
  zmax::Float64             # Max income
  rmin::Float64             # Min interest rate
  rmax::Float64             # Max interest rate
  agrid::Array{Float64}    # Grid space for assets
  zgrid::Array{Float64}    # Grid space for incoomes
  aa::Array{Float64}        # Matrix for asset grid space
  zz::Array{Float64}        # Matrix for income grid space
  da::Float64               # Asset grid step size
  dz::Float64               # Income grid step size
  dz2::Float64              # Income grid step size squared
  Delta::Float64            # Time step-size
  mu::Array{Float64}       # Drift from Ito's lemma
  s2::Array{Float64}       # Std Dev from Ito's lemma
  chi::Array{Float64}      # elements in the transition matrix
  nu::Array{Float64}       # elements in the transition matrix
  zeta::Array{Float64}     # elements in the transition matrix
  Cswitch::SparseMatrixCSC{Float64,Int64}  # Transition matrix
  v0::Array{Float64}       # Initial guess for value function
  r0::Float64                 # Intial guess for interest rate
  maxit::Int64              # Max iterations for solving HJB
  crit::Float64
  crit_S::Float64           # Tolerance for asset mkt clearing
end

# Create model function
function Huggett_model(;gamma::Float64=2.0,
  rho::Float64=0.05,
  ou_var::Float64=0.1,
  ou_corr::Float64=0.5,
  I::Int64=100,
  J::Int64=100,
  Ir::Int64=40,
  amin::Float64=-0.15,
  amax::Float64=2.0,
  rmin::Float64=-0.05,
  rmax::Float64=0.045,
  maxit::Int64=100,
  crit::Float64=1e-6,
  crit_S::Float64=1e-5,
  Delta::Float64=1000.0,
  r0::Float64=0.03)

  ## Ornstein-Uhlenbeck: dlog(z) = -theta*log(z)dt + sig*dW
  zmean   = exp(ou_var)
  theta   = -log(ou_corr)
  sig2    = theta*ou_var
  ## Ornstein-Uhlenbeck: dz = theta(zmean - z)dt + sig*dW
  #theta   = -log(ou_corr)
  #ou_var  = sig2/(theta)
  #zmean   = 1

  zmin    = zmean*0.65
  zmax    = zmean*1.25
  agrid   = reshape(collect(linspace(amin,amax,I)),I,1)
  zgrid   = reshape(collect(linspace(zmin,zmax,J)),J,1)
  aa      = repmat(agrid,1,J)
  zz      = repmat(zgrid',I,1)
  da      = (amax-amin)/(I-1)
  dz      = (zmax-zmin)/(J-1)
  dz2     = dz^2
  ## Ornstein-Uhlenbeck: in logs
  mu      = (-theta*log(zgrid) + sig2/2).*zgrid
  s2      = sig2.*zgrid.^2
  ## Ornstein-Uhlenbeck: in LEVELS
  #mu      = theta*(zmean - zgrid)
  #s2      = sig2.*ones(1,J)

  chi     = s2/(2*dz2)
  nu      = - s2/dz2 - mu/dz
  zeta    = mu/dz + s2/(2*dz2)

  v0      = (zz + r0.*aa).^(1-gamma)/(1-gamma)/rho

  #This will be the upperdiagonal of the C_switch
  updiag    = zeros(I,1)
  for j in 1:J-1
      updiag=vcat(updiag, repmat([zeta[j]],I,1))
  end

  #This will be the center diagonal of the B_switch
  centdiag = repmat([chi[1]+nu[1]],I,1)
  for j in 2:J-1
      centdiag = vcat(centdiag, repmat([nu[j]],I,1))
  end
  centdiag = vcat(centdiag, repmat([nu[J]+zeta[J]],I,1))

  #This will be the lower diagonal of the B_switch
  lowdiag = repmat([chi[2]],I,1)
  for j in 3:J
      lowdiag = vcat(lowdiag, repmat([chi[j]],I,1))
  end

  # Add upper, center, and lower diagonal into a sparse matrix
  Cswitch     = spdiagm(centdiag[:,1],0,I*J,I*J) + spdiagm(lowdiag[:,1],-I,I*J,I*J) + spdiagm(updiag[I+1:end,1],I,I*J,I*J)

  Huggett_model(gamma, rho, ou_var, zmean, ou_corr, theta, sig2, I, J, Ir, amin, amax, zmin, zmax, rmin, rmax, agrid, zgrid, aa, zz, da, dz, dz2, Delta, mu, s2, chi, nu, zeta, Cswitch, v0, r0, maxit, crit, crit_S)
end

# Creates an instance of the model
mod = Huggett_model()

# Define Utility funciton
function util(x; gamma::Float64=mod.gamma)
  u = x.^(1-gamma)/(1-gamma)
  return u
end

# Define derivative of Utility funciton
function utilprime(x; gamma::Float64=mod.gamma)
  uprime = x.^(-gamma)
  return uprime
end

# Define inverse of derivative of utility function
function utilprimeinv(x; gamma::Float64=mod.gamma)
  uprimeinv = x.^(-1/gamma)
  return uprimeinv
end

# Function to compute the KFE
function KFE(A::SparseMatrixCSC{Float64,Int64}; I::Int64=mod.I, J::Int64=mod.J, da::Float64=mod.da, dz::Float64=mod.dz)
  AT          = A'
  b           = zeros(I*J,1)
  i_fix       = 1    # fix one value so that matrix is not singular
  b[i_fix]    = 0.1
  row         = hcat(zeros(1,i_fix-1), 1, zeros(1,I*J-i_fix))
  AT[i_fix,:] = row
  #Solve linear system
  gg          = AT\b          # A'*gg = b --> gg = (A')^(-1)*b
  g_sum       = gg'*ones(I*J,1)*da*dz
  gg          = gg./g_sum     # Normalize so sum(gg)*da) = 1
  g           = reshape(gg,I,J)
  return g
end

# Initialization of value functions and consumption policy
Vaf         = zeros(mod.I,mod.J)
Vab         = zeros(mod.I,mod.J)
Vzf         = zeros(mod.I,mod.J)
Vzb         = zeros(mod.I,mod.J)
Vzz         = zeros(mod.I,mod.J)
c           = zeros(mod.I,mod.J)
r_r         = zeros(mod.Ir,1)
rmin_r      = zeros(mod.Ir,1)
rmax_r      = zeros(mod.Ir,1)
g_r         = zeros(mod.I,mod.J,mod.Ir)     # Stationary distribution
adot        = zeros(mod.I,mod.J,mod.Ir)    # Policy function
V_r         = zeros(mod.I,mod.J,mod.Ir)    # Value function
S           = zeros(mod.Ir,1)      # Aggregate asset supply
V           = zeros(mod.I,2)        # Value function
A           = sparse(zeros(mod.J*mod.I,mod.J*mod.I))   # Transition matrix
r           = mod.r0
dist        = zeros(mod.maxit,1)
r_end       = 0                   # Last iteration value

for ir in 1:mod.Ir
  r_r[ir]     = r
  rmin_r[ir]  = mod.rmin
  rmax_r[ir]  = mod.rmax

  if ir>1
    mod.v0 = V_r[:,:,ir-1]
  end

  v = mod.v0

  for n in 1:mod.maxit
    V               = v
    # forward difference
    Vaf[1:mod.I-1,:]    = (V[2:mod.I,:]-V[1:mod.I-1,:])/mod.da
    Vaf[mod.I,:]        = utilprime(mod.zgrid + r.*mod.amax)
    # backward difference
    Vab[1,:]            = utilprime(mod.zgrid + r.*mod.amin)  #state constraint boundary condition
    Vab[2:mod.I,:]      = (V[2:mod.I,:]-V[1:mod.I-1,:])/mod.da

    #consumption and savings with forward difference
    cf              = utilprimeinv(Vaf)
    sf             = mod.zz + r.*mod.aa - cf
    #consumption and savings with backward difference
    cb              = utilprimeinv(Vab)
    sb             = mod.zz + r.*mod.aa - cb
    #consumption and derivative of value function at steady state
    c0              = mod.zz + r.*mod.aa
    Va0             =utilprime(c0)

    # dV_upwind makes a choice of forward or backward differences based on the sign of the drift
    If              = sf .> 0    # positive drift --> forward difference
    Ib              = sb .< 0    # negative drift --> backward difference
    I0              = (1-If-Ib)   # at steady state
    #make sure backward difference is used at amax
    #Ib(I,:) = 1; If(I,:) = 0;
    #STATE CONSTRAINT at amin: USE BOUNDARY CONDITION UNLESS muf > 0:
    #already taken care of automatically

    Va_Upwind       = Vaf.*If + Vab.*Ib + Va0.*I0    #important to include third term
    c               = utilprimeinv(Va_Upwind)
    u               = util(c)

    #CONSTRUCT MATRIX
    X               = - min(sb,0)/mod.da
    Y               = - max(sf,0)/mod.da + min(sb,0)/mod.da
    Z               = max(sf,0)/mod.da

    updiag          = 0
    for j = 1:mod.J
        updiag      = vcat(updiag, Z[1:mod.I-1,j], 0)      # Note: Z_I,j = 0 for all j
    end
    updiag = updiag[2:end-1,1]

    centdiag = reshape(Y,mod.I*mod.J,1)

    lowdiag=X[2:mod.I,1]
    for j=2:mod.J
        lowdiag=vcat(lowdiag, 0, X[2:mod.I,j])          # Note: X_1,j = 0 for all j
    end

    B              = spdiagm(centdiag[:,1],0,mod.I*mod.J,mod.I*mod.J) + spdiagm(updiag[:,1],1,mod.I*mod.J,mod.I*mod.J) + spdiagm(lowdiag[:,1],-1,mod.I*mod.J,mod.I*mod.J)

    A              = B + mod.Cswitch
    InvMat         = (1/mod.Delta + mod.rho)*speye(mod.I*mod.J) - A

    u_stacked      = reshape(u,mod.I*mod.J,1)
    V_stacked      = reshape(V,mod.I*mod.J,1)

    b              = u_stacked + V_stacked/mod.Delta
    V_stacked      = InvMat\b  #SOLVE SYSTEM OF EQUATIONS
    V              = reshape(V_stacked,mod.I,mod.J)
    Vchange        = V - v
    v              = V

    dist[n] = maximum(maximum(abs(Vchange)))
    if dist[n]<mod.crit
      @printf("-----------------------\n")
      @printf("Interest rate = %f\n", r)
      @printf("Value Function Converged, Iteration = %i\n", n)
      break
    end
  end

  # Solve Kolmogorov Forward equation (Fokker-Plankck equation)
  g = KFE(A)
  g_r[:,:,ir]     = g
  adot[:,:,ir]    = mod.zz + r.*mod.aa - c
  V_r[:,:,ir]     = V

  S[ir]           = ((g'*mod.agrid*mod.da)'*ones(mod.I,1)*mod.dz)[1,1]

  #UPDATE INTEREST RATE
  r_end = ir

  if S[ir]>mod.crit_S
    @printf("Excess Supply at r = %0.2f\n", r)
    mod.rmax = r
    r = 0.5*(r+mod.rmin)
  elseif S[ir]<-mod.crit_S
    @printf("Excess Demand at r = %0.2f\n", r)
    mod.rmin = r
    r = 0.5*(r+mod.rmax)
  elseif abs(S[ir])<mod.crit_S
    @printf("Equilibrium Found, Interest rate = %0.2f\n", r)
    break
  end

end

## EXPORT DATA FOR GRAPHING IN PYTHON
using DataFrames
df = DataFrame(adot[:,:,r_end])
writetable("adot_diff.csv", df)

df = DataFrame(c)
writetable("cons_diff.csv", df)

df = DataFrame(g_r[:,:,r_end])
writetable("dist_diff.csv", df)

## GRAPHS
using Gadfly

# Distribution using a contour plot
plot(z=adot[:,:,r_end], x = mod.agrid, y = mod.zgrid, Geom.contour, Theme(line_width=1mm),Guide.title("Stationary distribution"),
Guide.xlabel("Wealth, a"),
Guide.ylabel("Income, z"),
Coord.cartesian(xmin=mod.amin, xmax=2), Theme(major_label_font_size=6mm,
minor_label_font_size=4mm))

# Savings Policy function
plot(layer(x=mod.agrid,y=adot[:,1,r_end],Geom.line,Theme(line_width=1.5mm)), layer(x=mod.agrid,y=adot[:,2,r_end],Geom.line, Theme(line_width=1.5mm,default_color=colorant"red")),
xintercept=[mod.amin], yintercept=[0], Geom.hline(color=colorant"black", size=0.5mm), Geom.vline(color=colorant"black", size=0.5mm),
Guide.title("Savings policy function"),
Guide.xlabel("Wealth, a"),
Guide.ylabel("Savings, si(a)"),
Coord.cartesian(xmin=mod.amin, xmax=mod.amax), Theme(major_label_font_size=6mm,
minor_label_font_size=4mm))

# Distributions
plot(layer(x=mod.agrid,y=g_r[:,50,r_end],Geom.line,Theme(line_width=1.5mm)), layer(x=mod.agrid,y=g_r[:,60,r_end],Geom.line, Theme(line_width=1.5mm,default_color=colorant"red")),
xintercept=[mod.amin], Geom.vline(color=colorant"black", size=0.5mm), yintercept=[0], Geom.hline(color=colorant"black", size=0.5mm),
Guide.title("Stationary distributions"),
Guide.xlabel("Wealth, a"),
Guide.ylabel("Density, gi(a)"),
Coord.cartesian(xmin=mod.amin, xmax=1), Theme(major_label_font_size=6mm,
minor_label_font_size=4mm))


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV
using JuMP
include("SupportFunction.jl")

#Set path
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/Diabetes64/"
cd(path)

#Set amount of columns
nCols = 262

#Set dataname
fName = "Diabetes64"

totalLog = zeros(30,nCols)

#Load AALR data, calculate max coorelation and output to file
for i = 1:10
    println(i)
    mainData = Array(CSV.read(string(i)*"Xtrain.csv", delim = ',', nullable=false))
    best3Beta = CSV.read(string(i)*"_"*fName*".csv", delim = ',', nullable=false, types = fill(Float64,nCols), header=false, datarow=2)

    best3BetaArr = Array(best3Beta)

    best3Beta[1,1] = findMaxCor(mainData, best3BetaArr[1,7:end])
    best3Beta[2,1] = findMaxCor(mainData, best3BetaArr[2,7:end])
    best3Beta[3,1] = findMaxCor(mainData, best3BetaArr[3,7:end])

    CSV.write(string(i)*fName*".csv", best3Beta)

    totalLog[1+(i-1)*3:3+(i-1)*3,:]= Array(best3Beta)
end
writedlm("TotalLog.csv", totalLog, ",")

#Initilise LASSO arrays
lassoBestK = zeros(nCols,nCols)
lassoSummary = zeros(nCols,nCols)
totalLog = zeros(30,nCols-1)

dataT = [Float64]
for j = 1:nCols-2
    dataT = vcat(dataT, Float64)
end

for i = 1:10
    #Load data
    println(i)
    mainData = Array(CSV.read(""*string(i)*"Xtrain.csv", delim = ',', nullable=false))
    dataInput = CSV.read(string(i)*"_Lasso"*fName*".csv", delim = ',', nullable=false, types = dataT, datarow=2, header=false)

    dataInputArr = Array{Float64}(dataInput)
    #Check if any lambda regulations results in j active prediciton terms
    for j = 1:nCols
        lassoBestK[j,1] = j
        lassoSummary[j,1] = j
        currentK = find(k -> (k==j), dataInputArr[:,4])
        if isempty(currentK)
            continue
        end

        #Find max calue within the places we have j active prediction terms
        val, index = findmax(dataInputArr[currentK,2])

        #Summary to be average over all datasets
        lassoSummary[j,2:end] += dataInputArr[currentK[index],:]

        #summary for current dataset
        lassoBestK[j,2:end] = dataInputArr[currentK[index],:]
        lassoBestK[j,2] = findMaxCor(mainData, lassoBestK[j,7:end])
    end

    #Find 3 highest R^2 values
    for j = 1:3
        val, index = findmax(dataInputArr[:,3])
        totalLog[j+(i-1)*3:j+(i-1)*3,:]= Array(dataInputArr[index, :])
        dataInputArr[index, 3] = 0
    end

    writedlm(string(i)*"LassoSummary.csv", lassoBestK,",")
end

writedlm("TotalLasso.csv", totalLog, ",")

lassoSummary[:,2:end] = lassoSummary[:,2:end]/10
writedlm("TotalLassoSummary.csv", lassoSummary,",")

ext: jl

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/DataAnalysis_AALRStatic.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
struct SRTM <: RasterDataSource end

# SRTM Mirror with 5x5 degree tiles
const SRTM_URI = URI(scheme = "https", host = "srtm.csi.cgiar.org", path = "/wp-content/uploads/files/srtm_5x5/TIFF")

const HAS_SRTM_TILE = BitArray([
    0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  1  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0
    0  1  1  1  1  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  1  1  1  1  1
    0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0
    0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
    0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  1  1  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0
    0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  1  0  0  0  0  0  0  0  0  1  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0
    1  1  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  0  0  0  0  0  0  0
    0  0  1  1  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0
    0  0  1  0  1  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  1  0  0  0  1  0  0
    0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  1  1  1  1  1  1  1  1  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  1  1  1  0  1  1  1  1  1  1  1  1  0  1  1  1  0  0  1  1  1  0
    0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  1  1  1  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0
    1  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  1  0  1  1  1  1  1  1  1  1  1  0  0  0  0  1  0  0  0  0  1  1  1  1  1  1  1  1  0  0  0  1  0  0  1  1  0
    0  1  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  0  1  1  0  0  1  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1  1  1
    0  1  0  1  1  1  0  1  1  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  1  0  0  0  0  1  1  1  1  1  1  0  1  1  1  0  0  1  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  0  1
    1  1  1  1  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  1  0  1  0  1  1  1  1  1  1  1  1  1  1  1  1  1
    1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  0  0  0  1  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1
    1  1  0  0  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  0  0  1  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  0  1  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  0
    1  0  0  0  0  0  0  1  0  1  0  0  0  0  1  0  0  0  0  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  0  0  1  0  0
    1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  1  0
    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  1  1  1  1  0  0  0  1  1
    1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  0  0  0  1  1  1
    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  1  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1
    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  0  0  1  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  1  0  0
    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
])

function _raster_tile_stem(tile_index::CartesianIndex)
    y, x = tile_index.I
    "srtm_$(lpad(x, 2, '0'))_$(lpad(y, 2, '0'))"
end

_rastername(::Type{SRTM}, tile_index::CartesianIndex{2}) = _raster_tile_stem(tile_index) * ".tif"
_rasterpath(T::Type{SRTM}, tile_index::CartesianIndex{2}) = joinpath(rasterpath(), "SRTM", _rastername(T, tile_index))

_zipname(::Type{SRTM}, tile_index::CartesianIndex{2}) = _raster_tile_stem(tile_index) * ".zip"
_zipurl(T::Type{SRTM}, tile_index::CartesianIndex{2}) = joinpath(SRTM_URI, _zipname(T, tile_index))
_zippath(T::Type{SRTM}, tile_index::CartesianIndex{2}) = joinpath(rasterpath(), "SRTM", "zips", _zipname(T, tile_index))

function _getraster(T::Type{SRTM}, tile_index::CartesianIndex{2})
    raster_path = _rasterpath(T, tile_index)
    if !isfile(raster_path)
        @info "Note: not all oceanic tiles exist in the SRTM dataset."
        zip_path = _zippath(T, tile_index)
        _maybe_download(_zipurl(T, tile_index), zip_path)
        mkpath(dirname(raster_path))
        raster_name = _rastername(T, tile_index)
        zf = ZipFile.Reader(zip_path)
        write(raster_path, read(_zipfile_to_read(raster_name, zf)))
        close(zf)
    end
    return raster_path
end

# Adapted from https://github.com/centreborelli/srtm4/blob/master/src/srtm4.c#L87-L117
function _wgs84_to_tile_x(x)
    # tiles longitude indexes go from 1 to 72,
    # covering the range from -180 to +180
    tile_x = (1 + floor(Int, (x + 180) / 5)) % 72
    return tile_x == 0 ? 72 : tile_x
end

function _wgs84_to_tile_y(y)
    y = clamp(y, -60, 60)
    tile_y = 1 + floor(Int, (60 - y) / 5)
    return tile_y == 25 ? 24 : tile_y
end

function bounds_to_tile_indices(::Type{SRTM}, bounds::NTuple{4,Real})
    bounds_to_tile_indices(SRTM, ((bounds[1], bounds[3]), (bounds[2], bounds[4])))
end
function bounds_to_tile_indices(::Type{SRTM}, (xs, ys)::NTuple{2,NTuple{2,Real}})
    _check_order(xs)
    _check_order(ys)
    t_xs = _wgs84_to_tile_x.(xs)
    t_ys = reverse(_wgs84_to_tile_y.(ys))
    @show t_xs t_ys
    return CartesianIndices((t_ys[1]:(t_ys[2]), t_xs[1]:(t_xs[2])))
end

_check_order((a, b)) = a > b && throw(ArgumentError("Upper bound $b less than lower bound $a"))

for op in (:getraster, :rastername, :rasterpath, :zipname, :zipurl, :zippath)
    _op = Symbol('_', op) # Name of internal function
    @eval begin
        # Broadcasting function dispatch
        function $_op(T::Type{SRTM}, tile_index::CartesianIndices)
            broadcast(tile_index) do I
                HAS_SRTM_TILE[I] ? $_op(T, I) : missing
            end
        end
        # Bounds to tile indices dispatch
        $_op(T::Type{SRTM}, bounds::Tuple) = $_op(T, bounds_to_tile_indices(T, bounds))

        # Public function definition with key-word arguments
        function $op(T::Type{SRTM}; bounds=nothing, tile_index=nothing)
            if isnothing(bounds) & isnothing(tile_index)
                :op === :getraster || return joinpath(rasterpath(), "SRTM")
                throw(ArgumentError("One of `bounds` or `tile_index` kwarg must be specified"))
            elseif !isnothing(bounds) & !isnothing(tile_index)
                throw(ArgumentError("Only one of `bounds` or `tile_index` should be specified. " *
                                    "found `bounds`=$bounds and `tile_index`=$tile_index"))
            else
                # Call the internal function without key-word arguments
                return $_op(T, isnothing(tile_index) ? bounds : tile_index)
            end
        end
    end
end

max_stars_repo_name: cesaraustralia/RasterDataSources

max_issues_repo_name: cesaraustralia/RasterDataSources

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
3 #run

reload("IBFS")

Sim = IBFS.Simulation()
Sim.configure()
Sim.setup()


text:
# ----------------------------------------------------------------------------------------------------------
# DRO models
# ----------------------------------------------------------------------------------------------------------
function test_generate(d::Projections.Divergence; atol::Real = 1e-4)
    Test.@testset "$(Projections.name(d))" begin
        ϕ = Projections.generate(d)
        Test.@test isapprox(ϕ(1), 0; atol = atol)
    end
    return
end

function test_feasibility(m::Projections.DRO, p; atol::Real = 1e-6)
    Test.@testset "Feasibility" begin
        Test.@test isapprox(sum(p), 1; atol = atol)
        Test.@test minimum(p) >= - atol
        if typeof(m.d) <: Projections.Divergence
            ϕ     = Projections.generate(m.d)
            Test.@test sum(m.q .* ϕ.(p ./ m.q)) <= 1.01 * m.ε
        else
            k     = Projections.normtype(m.d)
            Test.@test LinearAlgebra.norm(p - m.q, k) <= 1.01 * m.ε
        end
    end
    return
end

function isfeasible(m::Projections.DRO, p; atol::Real = 1e-6)

    cond1 = isapprox(sum(p), 1; atol = atol)
    cond2 = minimum(p) >= - atol
    if typeof(m.d) <: Projections.Divergence
        ϕ     = Projections.generate(m.d)
        cond3 = sum(m.q.*ϕ.(p./m.q)) <= 1.01 * m.ε
    else
        k     = Projections.normtype(m.d)
        cond3 = LinearAlgebra.norm(p - m.q, k) <= 1.01 * m.ε
    end
    return all([cond1, cond2, cond3])
end

function test_optimality(m::Projections.DRO, p, psolver; atol = 1e-4)
    isfeasible(m, psolver; atol = atol) || return

    Test.@testset "Optimality" begin
        Test.@test m.c'*p >= m.c'*psolver - atol
    end
    return
end

function test_model(m::Projections.DRO)
    Test.@testset "Comparison with the general solver" begin
        p       = Projections.solve(Projections.Our(), m);
        psolver = Projections.solve(Projections.General(), m);

        test_feasibility(m, p; atol = 1e-6)
        test_optimality(m, p, psolver; atol = 1e-4)

        if typeof(m.d) <: Projections.Ltwo
            p = Projections.solve(Projections.Philpott(), m);

            test_feasibility(m, p; atol = 1e-6)
            test_optimality(m, p, psolver; atol = 1e-4)
        end
    end
    return
end

# ----------------------------------------------------------------------------------------------------------
# Simplex models
# ----------------------------------------------------------------------------------------------------------
function test_feasibility(m::Projections.Simplex, p; atol::Real = 1e-6)
    Test.@testset "Feasibility" begin
        Test.@test all(p .>= m.lb .- atol)
        Test.@test all(p .<= m.ub .+ atol)
        Test.@test isapprox(sum(p), 1; atol = atol)
    end
    return
end

function isfeasible(m::Projections.Simplex, p; atol::Real = 1e-6)

    cond1 = all(p .>= m.lb .- atol)
    cond2 = all(p .<= m.ub .+ atol)
    cond3 = isapprox(sum(p), 1; atol = atol)

    return all([cond1, cond2, cond3])
end

function test_optimality(m::Projections.Simplex, p, psolver; atol = 1e-4)
    isfeasible(m, psolver; atol = atol) || return

    Test.@testset "Optimality" begin
        Test.@test LinearAlgebra.norm(p .- m.q)^2 <= LinearAlgebra.norm(psolver .- m.q)^2 + atol
    end
    return
end

function test_model(m::Projections.Simplex)
    Test.@testset "Comparison with the general solver" begin
        p       = Projections.solve(Projections.Our(), m);
        psolver = Projections.solve(Projections.General(), m);

        test_feasibility(m, p; atol = 1e-6)
        test_optimality(m, p, psolver; atol = 1e-4)
    end
    return
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
############################################################################################
# LDOS
############################################################################################

struct Fig4_presets
    η::Float64
    angle::Float64
    iterations::Int64
    nummodes::Int64
    selfy::Bool
end

"
    ldosonlattice_averaged_sc(p, fig_p::Fig4_presets; kw...)
average ldos weights of a number: Fig4_presets.nummodes of inner eigenvalues of a system
built with parameters p and  a given rotation angle Fig4_presets.angle and disorder
model Fig4_presets.angle.  The output is averaged overa a number of disorder realizations
specified by Fig4_presets.iterations.
    see: Fig4_presets, ldosonlattice_rand_sc, rectangle_randombounds_sc "

function ldosonlattice_averaged_sc(p, fig_p::Fig4_presets; kw...)
    ldosonlattice_averaged_sc(p, fig_p.η, fig_p.angle,
        fig_p.iterations; nummodes = fig_p.nummodes, selfy = fig_p.selfy, kw...)
end

function ldosonlattice_averaged_sc(p, η, angle = 0, iterations = 5; nummodes = 4,
    selfy = false, kw...)
    println(nummodes)
    psi1, psi2 = ldosonlattice_rand_sc(p, η, angle, nummodes = nummodes, selfy = selfy; kw...)
    # for i in 1:iterations
    #     psis = ldosonlattice_rand_sc(p, η, angle, nummodes = nummodes, selfy = selfy; kw...)
    #     psi1 .+= psis[1]
    #     psi2 .+= psis[2]
    # end
    # psi ./= iterations +1
    return psi1, psi2
end

"KPM calculation of localdensity of states at a site ind for the bounded system"
function ldos(params, ind, η, angle = 0.0; order = 300000, bandrange = missing, kw...)
    h = rectangle_randombounds_sc(params, angle, η; kw...);
    si = SMatrix{4,4}(Diagonal(SA[1,1,-1,-1]))
    km = ketmodel(si, indices = ind)
    dos = dosKPM(h, ket = km, bandrange = bandrange; order = order)
    return dos
end

"""
ldosonlattice_rand_sc(p, η, angle = 0; nummodes = 4, selfy = false, kw...)
returns the wavefunctions of nev inner eigenvalues for a system, created by
       rectangle_randombounds_sc()
"""
function ldosonlattice_rand_sc(p, η, angle = 0; nummodes = 4, selfy = false, kw...)
    ho = rectangle_randombounds_sc(p, angle, 0, selfy = selfy);
    hd = rectangle_randombounds_sc(p, angle, η, selfy = selfy; kw...);
    posd = sitepositions(hd)
    positionsdisordered(r) = r in posd ? true : false
    disord_indices = siteindices(ho, region = positionsdisordered)
    s = spectrum(hd, method = ArpackPackage(sigma = 1e-7, nev = ifelse(nummodes == 4, 4, 8)))
    println(s.energies)
    if nummodes == 4
        psi = Vector{Float64}(undef, Int(size(hd,1)*4))
        psi .=  abs.(s.states[:,1])
        psi .+= abs.(s.states[:,2])
        psi .+= abs.(s.states[:,3])
        psi .+= abs.(s.states[:,4])
        psi_ord = zeros(Float64, Int(size(ho,1)*4))
        count = 1
        for i in disord_indices
            psi_ord[4*i-3:4*i] = psi[count:count+3]
            count += 4
        end
        return psi_ord, psi_ord
    else nummodes == 8
        psi1 =  Vector{Float64}(undef, Int(size(hd,1)*4))
        psi2 =  Vector{Float64}(undef, Int(size(hd,1)*4))
        psi1 .=  abs.(s.states[:,1])
        psi1 .+= abs.(s.states[:,2])
        psi1 .+= abs.(s.states[:,7])
        psi1 .+= abs.(s.states[:,8])
        psi2 .=  abs.(s.states[:,3])
        psi2 .+= abs.(s.states[:,4])
        psi2 .+= abs.(s.states[:,5])
        psi2 .+= abs.(s.states[:,6])
        psi_ord1 = zeros(Float64, Int(size(ho,1)*4))
        psi_ord2 = zeros(Float64, Int(size(ho,1)*4))
        count = 1
        for i in disord_indices
            psi_ord1[4*i-3:4*i] = psi1[count:count+3]
            psi_ord2[4*i-3:4*i] = psi2[count:count+3]
            count += 4
        end
        return psi_ord1, psi_ord2
    end
end


text:
include("mpc_helpers.jl")
T = Float64
marker = [:none, :auto, :circle, :rect, :star5,
    :diamond, :hexagon, :cross, :xcross, :utriangle,
    :dtriangle, :rtriangle, :ltriangle, :pentagon,
    :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline, :x]

# Create players
p = 2
model = DoubleIntegratorGame(p=p)
model = UnicycleGame(p=p)
n = model.n
m = model.m

# Create Roadway
roadway_opts = MergingRoadwayOptions(lane_width=0.26, merging_angle=pi/12, merging_point=1.5)
roadway = build_roadway(roadway_opts)

# Create Players
players = Vector{Player{T}}(undef, p)
players[1] = Player(model, roadway.lane[5])
players[2] = Player(model, roadway.lane[3])

# Create Scenario
sce = Scenario(model, roadway, players)

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.12 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(1*[0.0, 0.5, 1.0, 0.5]) for i=1:p] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desired state
xf = [SVector{model.ni[1],T}([1.0,-0.10,0.0,0.50]),
    SVector{model.ni[2],T}([1.0,-0.10,0.0,0.75]),
    ]

# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
radius = 0.25*ones(p)
μ = 1.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)
game_con = GameConstraintValues(N, sce)
u_max =  SVector{m,T}([1.5, 1.5, 0.75, 0.75])
u_min =  SVector{m,T}([-1.5, -1.5, -0.75, -0.75])
add_control_bound!(game_con, u_max, u_min)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.25, -0.30,
   -0.45, -0.08,
   pi/12,   0.0,
    0.50,   0.55,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-8
opts.α_dual = 0.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

# Initialize visualizers
vis = [Visualizer() for i=1:p]
open.(vis)

# Visualize trajectories
col = [:cornflowerblue, :orange]
for i = 1:p
    set_scenario!(vis[i], sce, color=col)
    set_env!(vis[i], VehicleState(1.0, 0.0, 0.0, 0.0))
    set_camera_birdseye!(vis[i], height=10.0)
end

mpc_opts = MPCOptions(M=40, h=2.5)
m_stats = MPCStatistics()
t_stats = TrajStatistics(probsize)

for i = 1:20
    xinit = random_initial_state(model, x0, amplitude=1., seed=i)
    set_state!(vis[1], model, sce, xinit)
    sleep(0.1)
end

for i = 1:1
    xinit = random_initial_state(model, x0, seed=i, amplitude=1.0)
    mpc_loop(prob, xinit, sce, vis, mpc_opts, m_stats, t_stats, seed=[3,10],
        display_traj=true, display_player=true)
end

mpc_opts = MPCOptions(M=500, h=4.0)
mpc_stats, traj_stats = run_experiment(prob, x0, sce, vis, mpc_opts,
    n_sample=100, display_traj=true, display_player=true)
x, t, e, d = process_data(model, mpc_stats, traj_stats, α=0.25, upsampling=1, atol=2.0)
get_figure(x,t,e,d)

# Plot the mismatch
mpc_opts = MPCOptions(M=400, h=2.5)
m_stats = MPCStatistics()
t_stats = TrajStatistics(probsize)

xinit = random_initial_state(model, x0, seed=1, amplitude=1.0)
mpc_loop(prob, xinit, sce, vis, mpc_opts, m_stats, t_stats, seed=[3,10],
    display_traj=true, display_player=true)

t = 5
for i = 1:2
    build_waypoint!(vis[i], players, N, color=[col[i] for j=1:p])
    set_state!(vis[i], model, sce, Algames.state(m_stats.traj[t]))
    set_line_traj!(vis[i], model, sce.player, t_stats.traj[t][i], color=[col[i] for j=1:p])
    set_waypoint_traj!(vis[i], model, sce, t_stats.traj[t][i])
    set_camera_birdseye!(vis[i], height=5.0)
end
t = 40
cumsum(t_stats.dt)[t]
for i = 1:p
    set_state!(vis[i], model, sce, Algames.state(m_stats.traj[t]))
    set_line_traj!(vis[i], model, sce.player, t_stats.traj[t][i], color=col)
    set_waypoint_traj!(vis[i], model, sce, t_stats.traj[t][i])
end

# Visualization
highway_roadway_opts = HighwayRoadwayOptions(lane_width=0.26, lane_length=1.0)
highway_vis_opts = RoadwayVisualizationOptions(bound_height=-0.001)
k = 10
for i = 1:p
    set_roadway!(vis[i], highway_roadway_opts, k=k, vis_opts=highway_vis_opts)
    settransform!(vis[i]["env/roadway$k"], Translation(-0.4, 0.0, 0.0))
end


text:
#A thin SDL wrapper based on Starlight and Gamezero.

#https://github.com/jhigginbotham64/Starlight.jl/blob/main/src/SDL.jl
#https://github.com/aviks/GameZero.jl/blob/master/src/GameZero.jl

function initSDL(;buffers=4,samples=4,flags=MIX_INIT_FLAC|MIX_INIT_MP3|MIX_INIT_OGG)
    #From gamezero, but parameterized to make it more generic.
    SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, buffers)
    SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, samples)
    r = SDL_Init(UInt32(SDL_INIT_VIDEO | SDL_INIT_AUDIO))
    if r != 0
        error("Unable to initialise SDL: $(getSDLError())")
    end
    TTF_Init()

    inited = Mix_Init(Int32(flags))
    if inited & flags != flags
        @warn "Failed to initialise audio mixer properly. All sounds may not play correctly\n$(getSDLError())"
    end

    device = Mix_OpenAudio(Int32(22050), UInt16(MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
    if device != 0
        @warn "No audio device available, sounds and music will not play.\n$(getSDLError())"
        Mix_CloseAudio()
    end
end

function getSDLError()
    x = SDL_GetError()
    return unsafe_string(x)
end

function quitSDL()
    #This part is from gamezero and starlight.
    #TODO: make the game also destroy windows/everything before actually quitting SDL.
    Mix_HaltMusic()
    Mix_HaltChannel(Int32(-1))
    Mix_CloseAudio()
    TTF_Quit()
    Mix_Quit()
    SDL_Quit()
end


text:
# Residual equation for linear systems
linear_residual!(R, A, x, b) = R .= A * x - b

## Fully-simultaneous Newton method
#==========================================================================================#

# Coupled aero-structural-load-factor residuals
function coupled_residuals!(R, all_horsies, Γs, U, Ω, speed, stiffness_matrix, δs, fem_loads, weight, load_factor, L)
    # Get residual vector views
    R_A = R.aerodynamics
    R_S = R.structures
    R_W = @view R[end]

    # Aerodynamic residuals
    @timeit "Aerodynamic Residuals" solve_nonlinear!(R_A, all_horsies, Γs / speed, U / speed, Ω / speed)

    # Structural residuals
    @timeit "Structural Residuals" linear_residual!(R_S, stiffness_matrix, δs, fem_loads)

    # Weight residual
    linear_residual!(R_W, weight, load_factor, L)

    nothing
end

# Residual setup for multiple aerostructural surfaces and multiple aerodynamic surfaces
function solve_coupled_residual!(R, x, speed, β, ρ, Ω, syms :: Vector{Symbol}, chord_meshes, camber_meshes, fem_meshes, other_horsies, stiffness_matrix, weight, load_factor)
    # Unpack aerodynamic and structural variables
    Γs = x.aerodynamics
    δs = x.structures
    α  = x.load_factor

    # Compute velocity with new angle of attack
    U = freestream_to_cartesian(-speed, α, β)

    # Compute displacements
    Δs    = map((key, fem_mesh) -> reshape(δs[key][7:end], 6, length(fem_mesh)), valkeys(δs), fem_meshes)
    # @timeit "Get Translations" dxs   =
    # @timeit "Get Rotations" Ts    =

    # New VLM variables
    @timeit "New Horseshoes" new_horsies = @. new_horseshoes(mesh_translation(Δs), mesh_rotation(Δs), chord_meshes, camber_meshes, fem_meshes)

    # Compute component forces for structural residual
    @timeit "Get Circulations" new_Γs       = getindex.(Ref(Γs), syms)
    @timeit "Get Aerodynamic Centers" new_acs      = map(horsies -> bound_leg_center.(horsies), new_horsies)
    @timeit "All Horseshoes" all_horsies  = [ mapreduce(vec, vcat, new_horsies); vec(other_horsies) ]

    @timeit "New Forces" new_forces   = map((Γ_comp, hs_comp) -> surface_forces(hs_comp, Γ_comp, all_horsies, Γs, U, Ω, ρ), new_Γs, new_horsies)

    # Compute other forces for load factor residual
    @timeit "Get Symbols" other_syms   = filter(sym -> !(sym ∈ syms), keys(Γs))
    @timeit "Get Other Circulations" other_Γs     = mapreduce(sym -> vec(getindex(Γs, sym)), vcat, other_syms)
    @timeit "Other Forces" other_forces = surface_forces(vec(other_horsies), other_Γs, all_horsies, Γs, U, Ω, ρ)

    # Compute lift
    @timeit "All Forces" vlm_forces = [ mapreduce(vec, vcat, new_forces); vec(other_forces) ]
    @timeit "Transform Summed Forces" D, Y, L = geometry_to_wind_axes(sum(vlm_forces), α, β)

    # Build force vector with constraint for structures
    @timeit "FEM Loads" fem_loads = mapreduce(vec ∘ fem_load_vector, vcat, new_acs, new_forces, fem_meshes)

    # Compute residuals
    @timeit "Compute Residuals" coupled_residuals!(R, all_horsies, Γs, U, Ω, speed, stiffness_matrix, δs, fem_loads, weight, load_factor, L * cos(α))

    return R
end

# Residual setup for single aerostructural surface
function solve_coupled_residual!(R, x, speed, β, ρ, Ω, chord_mesh, camber_mesh, fem_mesh, stiffness_matrix, weight, load_factor)
    # Unpack aerodynamic and structural variables
    Γ = x.aerodynamics
    δ = x.structures
    α = x.load_factor

    # Compute velocity with new angle of attack
    U = freestream_to_cartesian(-speed, α, β)

    # Compute displacements
    Δs  = δ.displacement
    dxs = mesh_translation(Δs)
    Ts  = mesh_rotation(Δs)

    # New VLM variables
    @timeit "New Horseshoes" new_horsies = new_horseshoes(dxs, Ts, chord_mesh, camber_mesh, fem_mesh)

    # Compute aerodynamic forces
    @timeit "Surface Forces" vlm_forces = surface_forces(new_horsies, Γ, U, Ω, ρ)

    # Compute structural residual and loads
    @timeit "FEM Loads" fem_loads = fem_load_vector(bound_leg_center.(new_horsies), vlm_forces, fem_mesh)

    # Compute lift for load factor residual
    D, Y, L = geometry_to_wind_axes(sum(vlm_forces), α, β)

    # Compute residuals
    coupled_residuals!(R, new_horsies, Γ, U, Ω, speed, stiffness_matrix, δ, fem_loads, weight, load_factor, L * cos(α))

    return R
end

# # Residual setup for single aerostructural surface and multiple aerodynamic surfaces
# function solve_coupled_residual!(R, x, speed, β, ρ, Ω, chord_mesh, camber_mesh, other_horsies, fem_mesh, stiffness_matrix, weight, load_factor)
#     # Unpack aerodynamic and structural variables
#     Γ = x.aerodynamics
#     δ = x.structures
#     α = x.load_factor

#     # Get residual vector views
#     R_A = R.aerodynamics
#     R_S = R.structures
#     R_W = @view R[end]

#     # Compute velocity with new angle of attack
#     U = freestream_to_cartesian(-speed, α, β)

#     # Compute displacements
#     δs  = @views reshape(δ[7:end], 6, length(fem_mesh))
#     dxs = mesh_translation(Δs)
#     Ts  = mesh_rotation(Δs)

#     # New VLM variables
#     new_horsies = new_horseshoes(dxs, Ts, chord_mesh, camber_mesh, fem_mesh)
#     @timeit "Combine Horseshoes" all_horsies = [ vec(new_horsies); other_horsies ]

#     # Compute aerodynamic residual and loads
#     @timeit "Surface Forces" vlm_forces = surface_forces(Γ, all_horsies, U, Ω, ρ)

#     @timeit "FEM Forces" new_forces = @views reshape(vlm_forces[1:length(new_horsies)], size(new_horsies))

#     # Compute structural residual and loads
#     @timeit "FEM Loads" fem_loads = fem_load_vector(bound_leg_center.(new_horsies), new_forces, fem_mesh)

#     # Compute lift for load factor residual
#     D, Y, L = geometry_to_wind_axes(sum(vlm_forces), α, β)

#     # Compute residuals
#     coupled_residuals!(R, all_horsies, Γ, U, Ω, speed, stiffness_matrix, δ, fem_loads, weight, load_factor, L * cos(α))

#     return R
# end

# Residual setup for multiple aerostructural surfaces (NEED TO REDUCE REDUNDANCIES)
# function solve_coupled_residual!(R, x, speed, β, ρ, Ω, syms :: Vector{Symbol}, chord_meshes, camber_meshes, fem_meshes, stiffness_matrix, weight, load_factor)
#     # Unpack aerodynamic and structural variables
#     Γs = x.aerodynamics
#     δs = x.structures
#     α  = x.load_factor

#     # Compute velocity with new angle of attack
#     U = freestream_to_cartesian(-speed, α, β)

#     # Compute displacements
#     Δs    = map((key, n) -> reshape(δs[key][7:end], 6, n), valkeys(δs), length.(fem_meshes))
#     dxs   = mesh_translation.(Δs)
#     Ts    = mesh_rotation.(Δs)

#     # New VLM variables
#     new_horsies  = new_horseshoes.(dxs, Ts, chord_meshes, camber_meshes, fem_meshes)
#     all_horsies  = mapreduce(vec, vcat, new_horsies)

#     # Compute component forces for structural residual
#     new_Γs       = getindex.(Ref(Γs), syms)
#     new_acs      = map(horsies -> bound_leg_center.(horsies), new_horsies)
#     @timeit "New Forces" new_forces = surface_forces.(new_Γs, new_horsies, Ref(Γs), Ref(all_horsies), Ref(U), Ref(Ω), Ref(ρ))

#     # Compute lift
#     D, Y, L = geometry_to_wind_axes(sum(mapreduce(vec, vcat, new_forces)), α, β)

#     # Build force vector with constraint for structures
#     @timeit "FEM Loads" fem_loads = mapreduce(fem_load_vector, vcat, new_acs, new_forces, fem_meshes)

#     # Compute residuals
#     coupled_residuals!(R, all_horsies, Γs, U, Ω, speed, stiffness_matrix, δs, fem_loads, weight, load_factor, L * cos(α))

#     return R
# end

## Nonlinear block Gauss-Seidel
#==========================================================================================#

function aerostruct_gauss_seidel(x0, speed, β, ρ, Ω, chord_mesh, camber_mesh, fem_mesh, stiffness_matrix, weight, load_factor; max_iters = 50, tol = 1e-9)
    x = deepcopy(x0)
    ε = 1e5
    for i = 1:max_iters
        xp = deepcopy(x)

        @show xp

        Γ = @views x.aerodynamics
        δ = @views x.structures
        α = @views x.load_factor

        # Compute velocity with new angle of attack
        U = freestream_to_cartesian(-speed, α, β)

        # Compute displacements
        δs  = δ.displacement
        dxs = mesh_translation(δs)
        Ts  = mesh_rotation(δs)

        # New geometric variables
        new_horsies = new_horseshoes(dxs, Ts, chord_mesh, camber_mesh, fem_mesh)

        # Solve circulations
        Γ = reshape(influence_matrix(vec(new_horsies), -U / speed) \ boundary_condition(quasi_steady_freestream(vec(new_horsies), U, Ω), horseshoe_normal.(vec(new_horsies))), size(new_horsies))

        x.aerodynamics = Γ

        @show x

        # Compute VLM forces
        vlm_forces = surface_forces(Γ, new_horsies, U, Ω, ρ)

        # Compute FEM loads
        fem_loads = fem_load_vector(bound_leg_center.(new_horsies), vlm_forces, fem_mesh)

        # Solve displacements
        δ = reshape(stiffness_matrix \ vec(fem_loads), size(fem_loads))

        x.structures = δ

        @show x

        # Compute lift
        D, Y, L = geometry_to_wind_axes(sum(vlm_forces), α, β)

        # Weight residual
        α = acos(abs(weight * load_factor - L) / (weight * load_factor))

        x.load_factor = α

        @show L

        ε    = LinearAlgebra.norm(x - xp)
        @show (i, ε)

        if ε <= tol return x end # Needs NAN checks and everything like NLsolve

    end

    return x
end


text:
module DungAnalyse

using Revise
# Dates, VideoIO, Combinatorics, UUIDs, Observables
using DungBase, Serialization, Dates, JLSO, DataDeps, ProgressMeter, StaticArrays
using DelimitedFiles
# using FixedPointNumbers, ColorTypes
# using Images, CoordinateTransformations, LinearAlgebra, StaticArrays
# import VideoIO

# export build_calibrations

register(DataDep("coffeebeetle", "the coffee beetle database", "https://s3.eu-central-1.amazonaws.com/vision-group-file-sharing/Data%20backup%20and%20storage/Yakir/coffee%20beetles/database.zip", "557ec5027d34d2641db4f2cc61c993edd2ce2bee1530254d45d26cc41c5adf7a", post_fetch_method = unpack))

run(`s3fs -o uid=1000,gid=1000,allow_other binarybeetle /home/yakir/s3-drive`)
run(`mount`)
coffeesource = "s3-drive/coffee source/"

# include("ffmpeg.jl")
include("/home/yakir/dungProject/DungAnalyse/src/ffmpeg.jl")
include("/home/yakir/dungProject/DungAnalyse/src/calibrate.jl")
include("/home/yakir/dungProject/DungAnalyse/src/smooth.jl")
include("/home/yakir/dungProject/DungAnalyse/src/common_methods.jl")
include("/home/yakir/dungProject/DungAnalyse/src/quality_report.jl")
include("/home/yakir/dungProject/DungAnalyse/src/plotresults.jl")

# _formatcalibration(x) = replace(string(filenames(x), "_", Time(0) + start(x)), r"[^\w\d\s]" => "_")

#=function build_calibrations(coffeesource)
    sourcefolder = joinpath(coffeesource, "database")
    cs = deserialize(joinpath(sourcefolder, "calibration"))
    @showprogress 1 "saving calibration images…" for c in cs
        save_calibration_images(c, coffeesource)
    end
end

function build_calibrations(coffeesource, data::Dict{String,Experiment})
    @showprogress 1 "saving calibration images…" for (k,v) in data
        for r in v.runs, (k,p) in r.pois, (k,c) in p.calib
            build_calibration(coffeesource, k, c)
        end
    end
end=#

if !isfile(joinpath(datadep"coffeebeetle", "data.jlso"))
    include("/home/yakir/dungProject/DungAnalyse/src/load_from_csv.jl")
    loaddeomcsv()
end
data = JLSO.load(joinpath(datadep"coffeebeetle", "data.jlso"))
temporal2pixel = JLSO.load(datadep"coffeebeetle/temporal2pixel.jlso")["data"]

cs = Set{Calibration}()
for (_,v) in data, r in v.runs, (_,p) in r.data
    push!(cs, p.calib)
end
if !isfile(joinpath(coffeesource, "calibration errors.jlso"))
    nameerror = Dict{Calibration, NamedTuple{(:filename, :error),Tuple{String,Float64}}}()
    @showprogress 1 "Calibrating…" for c in cs
        nameerror[c] = build_calibration(coffeesource, hash(c), c)
    end
    JLSO.save(joinpath(coffeesource, "calibration errors.jlso"), nameerror)
end
nameerror = JLSO.load(joinpath(coffeesource, "calibration errors.jlso"))["data"]

function temp2pixel(k, p::POI{C, Temporal{V, I}}) where {C <: Calibration, V <: AbstractTimeLine, I <: Instantaneous}
    xyt = readdlm(joinpath(datadep"coffeebeetle/pixel", temporal2pixel[Pair(k, p)]))
    Instantaneous(xyt)
end
function temp2pixel(k, p::POI{C, Temporal{V, P}}) where {C <: Calibration, V <: AbstractTimeLine, P <: Prolonged}
    xyt = readdlm(joinpath(datadep"coffeebeetle/pixel", temporal2pixel[Pair(k, p)]))
    Prolonged(xyt)
end

pixeldata = Dict{String, Experiment}()
@showprogress 1 "loading pixel data…" for (experimentid, v) in data
    pixeldata[experimentid] = Experiment(Run[], v.description)
    for r in v.runs
        pois = Dict{Symbol, POI}()
        for (poitype, p) in r.data
            pois[poitype] = POI(nameerror[p.calib], temp2pixel(poitype, p))
        end
        push!(pixeldata[experimentid].runs, Run(pois, r.metadata))
    end
end
cmdata = Dict{String, Experiment}()
@showprogress 1 "Calibrating…" for (experimentid, v) in pixeldata
    cmdata[experimentid] = Experiment(Run[], v.description)
    for r in v.runs
        pois = Dict{Symbol, Any}()
        for (poitype, p) in r.data
            pois[poitype] = calibrate(p.calib.filename, p.data)
        end
        push!(cmdata[experimentid].runs, Run(pois, r.metadata))
    end
end
if !isfile(joinpath(coffeesource, "trackdata.jlso"))
    trackdata = Dict{String, Experiment}()
    @showprogress 1 "Tracking…" for (experimentid, v) in cmdata
        trackdata[experimentid] = Experiment([Run(Common(r), r.metadata) for r in v.runs], v.description)
        # break
    end
    JLSO.save(joinpath(coffeesource, "trackdata.jlso"), trackdata)
end
trackdata = JLSO.load(joinpath(coffeesource, "trackdata.jlso"))

fixpath(x) = replace(x, r"[\#\:\-]" => "_")

# experimentid = "nest#closed person#therese"
# v = cmdata[experimentid]
# trackdata = Dict{String, Experiment}()
# trackdata[experimentid] = Experiment([Run(Common(r), r.metadata) for r in v.runs], v.description)

@showprogress 1 "plotting…" for (experimentid, v) in data
    targetfiles = [start(r.data[:track].data.video) for r in v.runs]
    for (i, r) in enumerate(v.runs)
        path = fixpath(joinpath(coffeesource, "results", experimentid))
        mkpath(path)
        # plotquality(r.data, pixeldata[experimentid].runs[i].data, nameerror, cmdata[experimentid].runs[i].data, trackdata[experimentid].runs[i].data.track, r.metadata.comment, joinpath(path, fixpath(string("quality_", targetfiles[i], ".pdf"))))
        plotresult(trackdata[experimentid].runs[i].data, joinpath(path, fixpath(string("result_", targetfiles[i], ".pdf"))))
        # break
    end
    # break
end

end # module

# experimentid = "transfer#far person#therese"
# v = data[experimentid]

# TODO
# change terminology to nest and imaginary nest
# make sure all the terms (like direction etc) are consistent
# define where the binary stops are for the pipeline
# consider how the code works on the experiment[run] versus just run[]
# check all runs
# filter away bad ones
# adjust turning point algorithm
# consider saving the calibration images in a tmp folder instead
# have a thorough walk through the code to fine tune stuff
# consider that calibrations might be better connected to the videos than to the temporals...
# create the regression plot for marie


text:
using CSV, DataFrames, Test, Plots, Binscatters
df = DataFrame(CSV.File(joinpath(dirname(pathof(Binscatters)), "../dataset/Cigar.csv")))
df.id1 = df.State
df.y = df.Sales
df.x1 = df.Price
df.z1 = df.Pimin
df.x2 = df.NDI
df.w = df.Pop

out = Binscatters.bin(df, @formula(y~x1))
@test minimum(out.x1) >= minimum(df.x1)
@test maximum(out.x1) <= maximum(df.x1)
#Binscatters.bin(df, @formula(y~x1), weights = :w)
#@test minimum(out.x1) >= minimum(df.x1)
#@test maximum(out.x1) <= maximum(df.x1)
Binscatters.bin(df, @formula(y~x1), 10)
Binscatters.bin(df, @formula(y~x1+x2))
Binscatters.bin(df, @formula(y ~ x1 + fe(id1)))
Binscatters.bin(df, @formula(y+x2~x1))

Binscatters.bin(groupby(df, :id1), @formula(y~x1))

binscatter(df, @formula(y ~ x1))
binscatter(df, @formula(y ~ x1), seriestype = :scatterpath)
binscatter(df, @formula(y ~ x1), seriestype = :linearfit)

binscatter(df, @formula(y + x2 ~ x1))
binscatter(df, @formula(y + x2 ~ x1), seriestype = :scatterpath)
binscatter(df, @formula(y + x2 ~ x1), seriestype = :linearfit)

binscatter(df, @formula(y ~ x1 + fe(id1)))
binscatter(df, @formula(y~x1), weights = :w)
binscatter(df, @formula(y+x2~x1))
binscatter(df, @formula(y~x1+x2))

df.dummy = df.id1 .>= 25
binscatter(groupby(df, :dummy), @formula(y~x1))
binscatter(groupby(df, :dummy), @formula(y~x1), seriestype = :scatterpath)
binscatter(groupby(df, :dummy), @formula(y~x1), seriestype = :linearfit)

binscatter(groupby(df, :dummy), @formula(y+x2~x1))
binscatter(groupby(df, :dummy), @formula(y+x2~x1), seriestype = :scatterpath)
binscatter(groupby(df, :dummy), @formula(y+x2~x1), seriestype = :linearfit)


text:
using Test, NeXLUncertainties

include("uncertainValue.jl")
include("uncertainValues.jl")
include("propagate.jl")


text:
# start in REPL or with: julia htttp_webserver_simple.jl:
using HTTP

HTTP.listen() do http                       # 1
  while !eof(http)                          # 2
      println("body data: ", String(readavailable(http)))   # 3
  end
  HTTP.setstatus(http, 200)
  HTTP.setheader(http, "Content-Type" => "text/html")  # 4
  HTTP.startwrite(http)                         # 5
  write(http, "ToDo 1: Getting groceries
")  # 6
  write(http, "ToDo 2: Visiting my therapist
")
  write(http, "ToDo 3: Getting a haircut")
end

# const host = # ip-address in string format, like "127.0.0.1"
# const port = 8081
# HTTP.listen(host, port) do http
#   	# code
# end


text:
#
# Types for orbital trajectories
#

"""
An alias for some abstract `ODESolution`
with `States` state vector and parameter
vector types.
"""
const AbstractOrbitalODESolution = SciMLBase.AbstractODESolution{T,N,<:AbstractVector{U}} where {T,N,U<:States.AbstractState}

"""
An wrapper for a `SciMLBase.ODESolution` with a `GeneralAstrodynamics.States.AbstractState`
state vector type. This represents an object's `Trajectory` in space!
"""
struct Trajectory{FR, S, P, E, O<:AbstractOrbitalODESolution}
    epoch::E
    solution::O
end

"""
Returns the start `epoch`. Typically,
This is a type defined in `AstroTime.Epochs`.
"""
initialepoch(traj::Trajectory) = traj.epoch

"""
Returns the `solution` for the `Trajectory`. Typically,
this is a `DifferentialEquations.ODESolution`.
"""
solution(traj::Trajectory) = traj.solution

"""
A wrapper for (::ODESolution-like)(args...). Returns a state
of type `initialstate)` at time `t` past `epoch`.
"""
(traj::Trajectory)(t::Number; idxs=nothing, continuity=:left) = traj.solution(t, Val{0}; idxs=idxs, continuity=continuity)
(traj::Trajectory)(t::AbstractVector{<:Number}; idxs=nothing, continuity=:left) = traj.solution(t, Val{0}; idxs=idxs, continuity=continuity)

"""
Returns the initial condition associated with the `Trajectory`.
"""
initialstate(traj::Trajectory) = solution(traj).prob.u0

"""
Returns the system associated with the `Trajectory`.
"""
States.system(traj::Trajectory) = solution(traj).prob.p

"""
Returns the `eltype` of the `Trajectory`.
"""
Base.eltype(traj::Trajectory) = eltype(solution(traj))

"""
The `length` of a `Trajectory`.
"""
Base.length(traj::Trajectory) = length(solution(traj))

"""
Returns the last index of a `Trajectory`.
"""
Base.lastindex(traj::Trajectory) = length(traj)

"""
Calls the underlying `solution`'s `getindex` function
to return the `CartesianState` of the `Trajectory`
at time `t` past the `initialepoch`.
"""
Base.getindex(traj::Trajectory, args...) = getindex(solution(traj), args...)

"""
The `size` of a `Trajectory`.
"""
Base.size(traj::Trajectory) = size(solution(traj))

"""
Returns the `OrbitalFrame` of the `Trajectory`.
"""
Base.@pure States.frame(::Trajectory{FR}) where FR = FR

"""
Returns the length unit for the `Trajectory`.
"""
States.lengthunit(traj::Trajectory) = lengthunit(initialstate(traj))

"""
Returns the time unit for the `Trajectory`.
"""
States.timeunit(traj::Trajectory) = timeunit(initialstate(traj))

"""
Returns the angular unit for the `Trajectory`.
"""
States.angularunit(traj::Trajectory) = angularunit(initialstate(traj))

"""
Returns the velocity unit for the `Trajectory`.
"""
States.velocityunit(traj::Trajectory) = velocityunit(initialstate(traj))

"""
Returns the mass unit for the `Trajectory`.
"""
States.massunit(traj::Trajectory) = massunit(system(traj))

"""
Returns the mass parameter unit for the `Trajectory`.
"""
States.massparamunit(traj::Trajectory) = massparamunit(system(traj))

"""
Returns the position of the `Trajectory` at `t` `timeunit`'s from
the `initialstate`'s `epoch`.
"""
Base.position(traj::Trajectory, t) = t isa Unitful.Time ? position(traj, t / timeunit(traj)) : traj(t, Val{0}; idxs=1:3, continuity=:left) * lengthunit(traj)

"""
Returns the scalar distance of the `Trajectory` at `t` `timeunit`'s from
the `initialstate`'s `epoch`.
"""
States.distance(traj::Trajectory, t) = t isa Unitful.Time ? distance(traj, t / timeunit(traj)) : norm(position(traj, t))

"""
Returns the position of the `Trajectory` at `t` `timeunit`'s from
the `initialstate`'s `epoch`.
"""
States.velocity(traj::Trajectory, t) = t isa Unitful.Time ? velocity(traj, t / timeunit(traj)) : traj(t, Val{0}; idxs=4:6, continuity=:left) * velocityunit(traj)

"""
Returns the scalar distance of the `Trajectory` at `t` `timeunit`'s from
the `initialstate`'s `epoch`.
"""
States.speed(traj::Trajectory, t) = t isa Unitful.Time ? speed(traj, t / timeunit(traj)) : norm(velocity(traj, t))

"""
Returns a `state` of type `typeof(initialstate)`
at time `t`.
"""
States.state(traj::Trajectory, t) = t isa Unitful.Time ? state(traj, t / timeunit(traj)) : typeof(initialstate(traj))(solution(traj)(t))

"""
Returns an `epoch` at time `t` past `initialepoch`.
"""
States.epoch(traj::Trajectory, t) = t isa Unitful.Time ? initialepoch(traj) + UnitfulToAstroTime(t) : epoch(traj, t * timeunit(traj))

"""
Converts `Unitful` types to `AstroTime` types.
Throws an `ArgumentError` if the unit of quantity
`t` is not a second, minute, hour, day, or year.
"""
function UnitfulToAstroTime(t::Unitful.Time)
    un  = Unitful.unit(t)
    val = ustrip(un, t)

    if un == u"s"
        return val * AstroTime.seconds
    elseif un == u"minute"
        return val * Astrotime.minutes
    elseif un == u"hr"
        return val * AstroTime.hours
    elseif un == u"d"
        return val * AstroTime.days
    elseif un == u"yr"
        return val * AstroTime.years
    else
        throw(ArgumentError("Invalid time unit!"))
    end
end

"""
Returns an `Orbit` at time `t`.
"""
States.Orbit(traj::Trajectory, t) = Orbit(state(traj, t), system(traj), epoch(traj, t))

"""
Show a `Trajectory`.
"""
function Base.show(io::IO, traj::Trajectory)
    println(io, "Trajectory with $(length(traj)) timesteps and eltype $(eltype(solution(traj)))")
end


text:
module WeavePlots

using Base64
import Plots
import Weave

"""Pre-execute hooks to set the plot size for the chunk """
function plots_set_size(chunk)
  w = chunk.options[:fig_width] * chunk.options[:dpi]
  h = chunk.options[:fig_height] * chunk.options[:dpi]
  Plots.default(size = (w,h))
  return chunk
end

Weave.push_preexecute_hook(plots_set_size)

#PNG or SVG is not working, output html
function Base.display(report::Weave.Report, m::MIME"image/svg+xml", data::Plots.Plot{Plots.PlotlyBackend})#
  #Remove extra spaces from start of line for pandoc
  s = repr(MIME("text/html"), data)
  splitted = split(s, "\n")
  start = split(splitted[1], r"(?=<div)")
  #script = lstrip(start[1]) #local

  div = lstrip(start[2])
  plot = join(map(lstrip, splitted[2:end]), "\n")

  if report.first_plot
    report.header_script *= "<script src=\"https://cdn.plot.ly/plotly-latest.min.js\"></script>"
    report.first_plot = false
  end

    report.rich_output *= "\n" * div * "\n" * plot
end

function Base.display(report::Weave.Report, m::MIME"image/png", data::Plots.Plot{Plots.PlotlyBackend})#
  display(report, MIME("image/svg+xml"), data)
end

#PNG or SVG is not working, output html
function Base.display(report::Weave.Report, m::MIME"image/svg+xml", plot::Plots.Plot{Plots.PlotlyJSBackend})
  body = Plots.PlotlyJS.html_body(plot.o.plot)

  if report.first_plot
    report.header_script *= "<script src=\"https://cdn.plot.ly/plotly-latest.min.js\"></script>"
    report.first_plot = false
  end

  report.rich_output *= "\n" * body
end

function Base.display(report::Weave.Report, m::MIME"image/png", plot::Plots.Plot{Plots.PlotlyJSBackend})
  display(report, MIME("image/svg+xml"), data)
end

"""Add saved figure name to results and return the name"""
function add_plots_figure(report::Weave.Report, plot::Plots.Plot, ext)
  chunk = report.cur_chunk
  full_name, rel_name = Weave.get_figname(report, chunk, ext = ext)

  Plots.savefig(plot, full_name)
  push!(report.figures, rel_name)
  report.fignum += 1
  return full_name
end

function Base.display(report::Weave.Report, m::MIME"application/pdf", plot::Plots.Plot)
    add_plots_figure(report, plot, ".pdf")
end

function Base.display(report::Weave.Report, m::MIME"image/png", plot::Plots.Plot)
    add_plots_figure(report, plot, ".png")
end

function Base.display(report::Weave.Report, m::MIME"image/svg+xml", plot::Plots.Plot)
    add_plots_figure(report, plot, ".svg")
end

# write out html to view Animated gif
function Base.display(report::Weave.Report, ::MIME"text/html", agif::Plots.AnimatedGif)
  ext = agif.filename[end-2:end]
  res = ""
  if ext == "gif"
      img = stringmime(MIME("image/gif"), read(agif.filename))
      res = "<img src=\"data:image/gif;base64,$img\" />"
  elseif ext in ("mov", "mp4")
      #Uncomment to embed mp4, make global or chunk option?
      #img = stringmime(MIME("video/$ext"), read(agif.filename))
      #res = "<video controls><source src=\"data:video/$(ext);base64,$img\" type=\"video/$ext\"></video>"
      res = "<video controls><source src=\"$(relpath(agif.filename))\" type=\"video/$ext\"></video>"
  else
      error("Cannot show animation with extension $ext: $agif")
  end

  report.rich_output *= "\n" * res * "\n"
end

end


text:
using DataFrames
using PyCall

global const pyModPandas = PyCall.PyNULL()
global const pyModxlsxwriter = PyCall.PyNULL()

copy!(pyModPandas, PyCall.pyimport_conda("pandas", "pandas"))
copy!(pyModxlsxwriter, PyCall.pyimport_conda("xlsxwriter", "xlsxwriter"))

mutable struct Chart
	sheet::String
	chartDict::Dict{AbstractString,Dict{AbstractString,Any}}
	location::String
end

mutable struct ExcelSheet
	name::String
	data::DataFrame
end

mutable struct ExcelData
	sheets::Array{ExcelSheet,1}
    charts::Array{Chart,1}
    function ExcelData()
        return new(Array{ExcelSheet}(undef,0),Array{Chart}(undef,0))
    end
    function ExcelData(a,b)
        return  new(a,b)
    end
end

function create_custom_dict(df::DataFrame)
	header = propertynames(df)
	d = Dict{AbstractString,Array{Any,1}}()
    for i = 1:length(header)
		d[string(header[i])] = df[!,i]
	end
	return d
end

function writeDFtoExcel(excelData::ExcelData, existingFile::T, row::Int, col::Int, write_header::Bool, write_index::Bool) where {T <: AbstractString}
    # http://search.cpan.org/~jmcnamara/Excel-Writer-XLSX/lib/Excel/Writer/XLSX.pm
        @assert min(row, col) >= 0
        writer = pyModPandas.ExcelWriter(existingFile, engine="xlsxwriter")

        for xlSheet in excelData.sheets
            df = xlSheet.data
            sheet = xlSheet.name
            # create python dataframe
            dataDict = create_custom_dict(df)
            pyDF = PyCall.pycall(pyModPandas.DataFrame, PyCall.PyObject, dataDict, columns=propertynames(df))
            PyCall.pycall(pyDF."to_excel", PyCall.PyAny, writer, header=write_header, index=write_index, sheet_name=sheet, startrow=row, startcol=col, encoding="utf-8")  # index=false suppress the rowcount
        end
        return writer
    end

chartDict=Dict{AbstractString,Dict{AbstractString,Any}}(
    "set_y_axis" => Dict{AbstractString,Any}("name"=>"Observed Ratio","major_gridlines"=>Dict{Any,Any}("visible"=>true)),
    "set_legend" => Dict{AbstractString,Any}("position"=>"right"),
    "set_size"   => Dict{AbstractString,Any}("y_scale"=>1.5,"x_scale"=>2.4),
    "series1"    => Dict{AbstractString,Any}("name"=>"=ModelStatistics!\$A\$1","values"=>"=ModelStatistics!\$F\$3:\$F\$6","categories"=>"=ModelStatistics!\$A\$3:\$A\$6"),
    "add_chart"  => Dict{AbstractString,Any}("type"=>"column")
  )

  aChart=Chart("ModelStatistics",chartDict,"K2")

  excelData=ExcelData()
  sheet1df=DataFrame(rand(10,10), :auto)
  xlSheet=ExcelSheet("ModelStatistics",sheet1df)

  push!(excelData.sheets,xlSheet)
  push!(excelData.charts,aChart)

  dir=mktempdir()
  statsfile=joinpath(dir,"anExcelFile.xlsx")

  writer = writeDFtoExcel(excelData, statsfile, 0, 0, true,false)
  workbook = writer.book
  chart = PyCall.pycall(workbook."add_chart", PyCall.PyAny, chartDict["add_chart"])

for x in keys(chartDict)
    fieldsWhichAreAlreadySet = [convert(String, string("series", zz)) for zz in 1:5]
    resevedKeywords = ["combChart","series_comb1","series_comb2","series_comb3", "series_comb4"] # currently limited to 1+4 series for combined charts
    append!(fieldsWhichAreAlreadySet, ["add_series","add_chart"])
    append!(fieldsWhichAreAlreadySet, resevedKeywords)
    if !in(x, fieldsWhichAreAlreadySet)
        # TBD: unclear how to fix this line... (deprecated syntax...)
        y=chartDict[x]
        @show x
        @show y
        PyCall.pycall(chart[x], PyCall.PyAny, y)  #works (unclear if this syntax is indeed deprecated or not; if so, how to fix it?)
        # TBD: unclear how to fix this line...
        #PyCall.pycall(chart.$(x), PyCall.PyAny, y) #does NOT work
        #PyCall.pycall(chart.x, PyCall.PyAny, y)  #does NOT work
        #PyCall.pycall(chart."x", PyCall.PyAny, y)  #does NOT work
        #PyCall.pycall(chart."$(x)", PyCall.PyAny, y)  #does NOT work
        #PyCall.pycall("chart.$(x)", PyCall.PyAny, y)  #does NOT work
    end
end

function addChartToWorkbook!(workbook::PyCall.PyObject, worksheet::PyCall.PyObject, chartDict::Dict{AbstractString,Dict{AbstractString,Any}}, location::AbstractString) # ;properties=["set_x_axis", "set_y_axis","set_legend"])
	chart = PyCall.pycall(workbook."add_chart", PyCall.PyAny, chartDict["add_chart"])
	stopboolean = true
	i = 1
	local thiskey
	while stopboolean
		thiskey = string("series", i)
		if haskey(chartDict, thiskey)
			PyCall.pycall(chart."add_series", PyCall.PyAny, chartDict[thiskey])
		else
			stopboolean = false
			break # if there is no series2 then we assume there is no series 3 to n either
		end
		i += 1
	end
	# check if this is a combined chart
	thiskey = "combChart"
	i2 = 1
	stopboolean = true
	if haskey(chartDict, thiskey)
		second_chart = PyCall.pycall(workbook."add_chart", PyCall.PyAny, chartDict[thiskey])
		while stopboolean
			thiskey = string("series_comb", i2)
			if haskey(chartDict, thiskey)
				PyCall.pycall(second_chart."add_series", PyCall.PyAny, chartDict[thiskey])
			else
				stopboolean = false
				break # if there is no series2 then we assume there is no series 3 to n either
			end
			i2 += 1
		end
		# combine the charts
		PyCall.pycall(chart."combine", PyCall.PyAny, second_chart)
	end
	# set other properties
	for x in keys(chartDict)
		fieldsWhichAreAlreadySet = [convert(String, string("series", zz)) for zz in 1:i]
		resevedKeywords = ["combChart","series_comb1","series_comb2","series_comb3", "series_comb4"] # currently limited to 1+4 series for combined charts
		append!(fieldsWhichAreAlreadySet, ["add_series","add_chart"])
		append!(fieldsWhichAreAlreadySet, resevedKeywords)
		if !in(x, fieldsWhichAreAlreadySet)
            # TBD: unclear how to fix this line...
            # TBD: unclear how to fix this line... (deprecated syntax...)
            @show x
            PyCall.pycall(chart[x], PyCall.PyAny, chartDict[x])	 # TBD: unclear how to fix this line...
		end
	end
	PyCall.pycall(worksheet."insert_chart", PyCall.PyAny, location, chart)
	# writer[:save]()
end

function writeStatistics(excelData::ExcelData, statsfile::T, write_header::Bool, write_index::Bool) where {T <: AbstractString}
	# writing an Excel file seems very slow if the file already exists!
	isfile(statsfile) && rm(statsfile)

	writer = writeDFtoExcel(excelData, statsfile, 0, 0, write_header, write_index)
	workbook = writer.book
	# Plot charts
	for c in excelData.charts
		sheetWhereChartIsLocated = c.sheet
		worksheet = writer.sheets[sheetWhereChartIsLocated]
		addChartToWorkbook!(workbook, worksheet, c.chartDict, c.location);
	end
	# save (=write) Excel file and close it
	writer.save()
	return nothing
end

#to write the Excel file:
  writeStatistics(excelData,statsfile,true,false)
  @show statsfile


text:
# ADMM methods

# NOTE: I was not able to run either of those methods
# one thing that wasn't done properly is to re-factorize H every-time ρ is changed
# but that makes these algorithms super inefficient.
# without modifying ρ, I was not able to get the algorithms not to explode. Interestingly
# when running a litteral translation of https://web.stanford.edu/~boyd/papers/admm/least_abs_deviations/lad.html
# it also exploded.

# function _fit(glr::GLR{L1Loss,<:L2R}, solver::ADMM, X, y)
#     n, p = size(X)
#     λ    = get_penalty_scale(glr, n)
#     φ    = 1.0 / solver.rho
#     λφ   = λ * φ
#     # pre-computations
#     H = form_XtX(X, glr.fit_intercept, λφ) # O(np²)
#     cholesky!(H) # O(p³) important assumption p << n
#     # cache
#     p_  = p + Int(glr.fit_intercept) # effective dim
#     θ   = zeros(p_)    # cache for current  θ
#     θ_  = zeros(p_)    # cache for previous θ
#     b   = zeros(p_)    # cache for right-hand-side of θ-update
#     bₐ  = view(b, 1:p) # only used in case of glr.fit_intercept
#     Xθ  = zeros(n)     # cache for over-relaxed X*θ
#     z   = zeros(n)     # cache for complimentary variable
#     u   = zeros(n)     # cache for dual variable
#     t   = zeros(n)     # cache for z .+ y .- u
#     zpy = copy(y)      # cache for z + y which is used multiple times
#     # loop-related
#     k, tol = 1, Inf
#     while k ≤ solver.max_iter && tol > solver.tol
#         # θ-update
#         # >> θ = (X'X + λ/ρ I) \ X'(y + z - u)  -- O(max{np, p²})
#         t .= zpy .- u
#         if glr.fit_intercept
#             mul!(bₐ, X', t)
#             b[end] = sum(t)
#         else
#             mul!(b, X', t)
#         end
#         ldiv!(θ, H, b)
#         # over-relaxation
#         # >> Xθ = α * (X * θ) + (1-α) * (z + y)   -- O(np)
#         apply_X!(Xθ, X, θ)
#         Xθ .= solver.alpha .* Xθ .+ (1.0 - solver.alpha) .* zpy
#         # z-update          -- O(n)
#         z  .= soft_thresh.(Xθ .- y .+ u, φ)
#         # u-update
#         zpy .= z .+ y
#         u  .+= Xθ .- zpy
#         # update tolerance
#         tol = norm(θ .- θ_) / (norm(θ) + eps())
#         copyto!(θ_, θ)
#         # update niter
#         k += 1
#     end
#     tol ≤ solver.tol || @warn "ADMM did not converge in $(solver.max_iter)."
#     return θ
# end
#
#
# # FADMM methods
#
# function _fit(glr::GLR{L1Loss,<:L2R}, solver::FADMM, X, y)
#     n, p = size(X)
#     λ    = get_penalty_scale(glr, n)
#     ρ    = solver.rho
#     η    = solver.eta   # linked to restart frequency
#     τ    = solver.tau   # linked to updating ρ
#     μ    = solver.mu    # linked to updating ρ
#     # XXX store that in solver
#     ϵ_abs = 1e-3
#     ϵ_rel = 1e-3
#     # pre-computations
#     H = form_XtX(X, glr.fit_intercept, λ/ρ) # O(np²)
#     cholesky!(H) # O(p³) important assumption p << n
#     # cache
#     p_   = p + Int(glr.fit_intercept) # effective dim
#     θ    = zeros(p_)    # cache for current  θ
#     b    = zeros(p_)    # cache for right-hand-side of θ-update
#     bₐ   = view(b, 1:p)
#     Xθ   = zeros(n)     # cache for X*θ
#     z    = zeros(n)
#     z_   = zeros(n)
#     zmz_ = zeros(n)     # cache for (z-z_) used in dual residual and accel
#     u    = zeros(n)
#     u_   = zeros(n)
#     ẑ    = zeros(n)
#     û    = zeros(n)
#     t    = zeros(n)     # cache for θ update
#     r    = zeros(n)     # cache for primal residual
#     s    = zeros(p_)    # cache for dual residual
#     sₐ   = view(s, 1:p)
#     α    = 1.0          # cache for acceleration step
#     c_   = Inf          # force accelerate on first step
#     # loop-related
#     k = 1
#     while k ≤ solver.max_iter
#         # ----------------------------------------
#         # compute ϵ-primal and ϵ-dual XXX can store precomp
#         ϵ_primal = sqrt(p_) * ϵ_abs + ϵ_rel * max(norm(Xθ), norm(z), norm(y))
#         apply_Xt!(bₐ, X, u)
#         ϵ_dual   = sqrt(n) * ϵ_abs  + ϵ_rel * norm(b) / ρ
#         # ----------------------------------------
#         # θ-update
#         # >> θ = argmin_θ L(θ) - ⟨û,Xθ⟩ + ρ/2 |y-Xθ+ẑ|₂² with L(θ) = λ|θ|₂²/2
#         # >> θ = (X'X + λ/ρ) \ X'(y + ẑ + û/ρ) = H \ X't
#         t .= y .+ ẑ .+ û ./ ρ
#         apply_Xt!(bₐ, X, t)     # X'(x + ẑ + û/ρ)
#         ldiv!(θ, H, b)          # H = (X'X+λ/ρ)
#         # ----------------------------------------
#         # z-update
#         # >> z = argmin_z P(z) + ⟨û,z⟩ + ρ/2 |y-Xθ+z|₂² with P(z) = |z|₁
#         # >> z = S_{1/ρ}(Xθ - y - û/ρ)
#         apply_X!(Xθ, X, θ)
#         z .= soft_thresh.(Xθ .- y .- û ./ ρ, 1.0/ρ)
#         # ----------------------------------------
#         # r and u-update
#         # >> u = û + ρ(y - Xθ + z)
#         r .= Xθ .- z .- y  # primal residual
#         u .= û .- ρ .* r
#         # ----------------------------------------
#         # c-update
#         # >> c = 1/ρ norm(u - û)^2 + ρ norm(z - ẑ)^2
#         # note that (u-û).^2 = ρ²r.^2
#         c = ρ * sum(abs2, r) + ρ * sum(abs2, z .- ẑ)
#         # ----------------------------------------
#         # accelerate or restart
#         zmz_ .= z .- z_
#         if c < η * c_
#             println("accel")
#             # accelerate
#             α_  = α #  store α_k
#             α   = (1.0 + sqrt(1.0 + 4.0 * α^2)) / 2.0  # α_{k+1} from α_k
#             ζ   = (α_ - 1.0) / α                       # (α_{k}-1) / α_{k+1}
#             ẑ  .= z .+ ζ .* zmz_
#             û  .= u .+ ζ .* (u .- u_)
#         else
#             println("restart")
#             # restart
#             α = 1.0
#             copyto!(ẑ, z_)
#             copyto!(û, u_)
#             c = c_ / η
#         end
#         # ----------------------------------------
#         # ρ-update depending on primal-dual gap
#         apply_Xt!(sₐ, X, zmz_) # dual residual -ρX'(z-z_) without the ρ
#         norm_r = norm(r)
#         norm_s = ρ * norm(s)
#         if norm_r / ϵ_primal > μ * norm_s / ϵ_dual
#             ρ *= τ
#         elseif norm_s / ϵ_dual > μ * norm_r / ϵ_primal
#             ρ /= τ
#         end
#         # ----------------------------------------
#         # cache update
#         c_  = c         # next c_{k-1} (used in restart step)
#         copyto!(z_, z)  # next z_{k-1} (used in acceleration step)
#         copyto!(u_, u)  # next u_{k-1} (used in acceleration step)
#         # check if convergence
#         norm_r ≤ ϵ_primal && norm_s ≤ ϵ_dual && break
#         # update niter
#         k += 1
#     end
#     k == solver.max_iter+1 && @warn "FADMM did not converge in $(solver.max_iter)."
#     return θ
# end


text:
# Visualize results that are saved in the results folder (names should match)

using PyPlot
using HDF5

# NAME AND POLICIES NAMES
fname="results/Experiment4arms"
names=["TrackAndStop","ChernoffBC","ChernoffPTS","ChernoffPTSOpt","KLLUCB","UGapE"]

# PARAMETERS
delta=0.1
N=100

clf()

# personalize bins
xdim=3
ydim=2

NBins=30
xmax=10000

xtxt=0.6*xmax
Bins=round(Int,linspace(1,xmax,NBins))

mu=h5read("$(fname)_$(names[1])_delta_$(delta)_N_$(N).h5","mu")
K=length(mu)

clf()
title("mu = $(mu)")

for j in 1:length(names)
    name="$(fname)_$(names[j])_delta_$(delta)_N_$(N).h5"
    FracNT=h5read(name,"FracNT")
    Draws=h5read(name,"Draws")
    Error=h5read(name,"Error")
    subplot(xdim,ydim,j)
    NbDraws=sum(Draws,2)'
    proportion=zeros(N,K)
    for k in 1:N
        proportion[k,:]=Draws[k,:]/sum(Draws[k,:])
    end
    prop=mean(proportion,1)
    MeanDraws=mean(NbDraws)
    StdDraws=std(NbDraws)
    histo=plt[:hist](vec(NbDraws),Bins)
    Mhisto=maximum(histo[1])
    PyPlot.axis([0,xmax,0,Mhisto])
    ytxt1=0.75*Mhisto
    ytxt2=0.6*Mhisto
    ytxt3=0.5*Mhisto
    ytxt4=0.4*Mhisto
    EmpError=round(Int,10000*mean(Error))/10000
    FracReco=round(Int,1000*prop)/1000
    axvline(MeanDraws,color="black",linewidth=2.5)
    PyPlot.text(xtxt,ytxt1,"mean = $(round(Int,MeanDraws)) (std=$(round(Int,StdDraws)))")
    PyPlot.text(xtxt,ytxt2,"delta = $(delta)")
    PyPlot.text(xtxt,ytxt3,"emp. error = $(EmpError)")
    PyPlot.text(xtxt,ytxt4,"emp. proportions = $(FracReco)")
    if (j==1)
       title("mu=$(mu), $(names[j])")
    else
       title("$(names[j])")
    end
    print("Results for $(names[j]), average on $(N) runs\n")
    print("proportion of runs that did not terminate: $(FracNT)\n")
    print("average number of draws: $(MeanDraws)\n")
    print("average proportion of draws: \n $(prop)\n")
    print("proportion of errors: $(EmpError)\n")
    print("proportion of recommendation made when termination: $(FracReco)\n\n")

end

max_stars_repo_name: jsfunc/best-arm-identification

max_issues_repo_name: ClaireVern/best-arm-identification

max_forks_repo_name: ClaireVern/best-arm-identification


text:
__precompile__()

module MultivariateTests

using Compat
using Compat.LinearAlgebra
using StatsBase
using Distributions
using HypothesisTests

using HypothesisTests: HypothesisTest

export
    genvar,
    totalvar,
    PartialCorTest,
    partialcor,
    OneSampleHotellingT2,
    EqualCovHotellingT2,
    UnequalCovHotellingT2,
    BartlettsTest

## Common utility functions

function checkdims(X::AbstractMatrix, Y::AbstractMatrix)
    nx, p = size(X)
    ny, q = size(Y)
    p == q || throw(DimensionMismatch("Inconsistent number of variables: $p, $q"))
    (nx > 0 && ny > 0) || throw(ArgumentError("Inputs must be non-empty"))
    return (p, nx, ny)
end

# Pool the given covariance matrices, trashing them both
function _poolcov!(Sx::AbstractMatrix, nxm1::Int, Sy::AbstractMatrix, nym1::Int)
    Sx .*= nxm1
    Sy .*= nym1
    Sx .+= Sy
    Sx ./= nxm1 + nym1
    return Sx
end

# Benchmarks suggest that this is faster than using pivoted Cholesky when the input
# is not positive definite
function trychol(X::AbstractMatrix)
    S = Symmetric(X)
    return isposdef(S) ? cholfact(S) : factorize(S)
end

# Helper function for computing A'B⁻¹A when B is a covariance matrix, taking advantage
# of the fact that covariance matrices are symmetric, positive (semi-)definite. We'll
# check for positive definiteness and use Cholesky, otherwise let Julia pick Bunch-
# Kaufman or LU in the case of failure.
At_Binv_A(A::AbstractArray, B::AbstractArray) = A'*(trychol(B) \ A)

include("dispersion.jl")
include("partialcor.jl")
include("hotelling.jl")
include("covariance.jl")

end # module


text:
module MimiSSPs

global g_datasets = Dict{Symbol,Any}()

include("components/SSPs.jl") # contains data for 1750 to 2500 for emissions and 2010 to 2500 for socioeconomic
# include("components/RegionAggregatorSum.jl")

export SSPs, RegionAggregatorSum

end # module


text:
using Test
using Logging
# Package Under Test
using CRCBS
using JuMP
using LightGraphs, MetaGraphs
using Parameters
using GraphUtils

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test isapprox(a,b, rtol=rtol, atol=atol)
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        @test isapprox(a, y, rtol=rtol, atol=atol)
    end
end

set_global_verbosity!(0)
# Define package tests
@time @testset "CRCBS.Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "CRCBS.ProblemDefinitionTests" begin
        include(joinpath(testdir, "unit_tests/test_problem_definitions.jl"))
    end
    @time @testset "CRCBS.InterfaceTests" begin
        include(joinpath(testdir, "unit_tests/test_interface.jl"))
    end
    @time @testset "CRCBS.CostModelTests" begin
        include(joinpath(testdir, "unit_tests/test_cost_models.jl"))
    end
    @time @testset "CRCBS.HeuristicTests" begin
        include(joinpath(testdir, "unit_tests/test_heuristics.jl"))
    end
    @time @testset "CRCBS.CommonTests" begin
        include(joinpath(testdir, "unit_tests/test_common.jl"))
    end
    @time @testset "CRCBS.UtilsTests" begin
        include(joinpath(testdir, "unit_tests/test_utils.jl"))
    end
    @time @testset "CRCBS.ImplicitGraphsTests" begin
        include(joinpath(testdir, "unit_tests/test_a_star.jl"))
    end
    @time @testset "CRCBS.SolverUtilsTests" begin
        include(joinpath(testdir, "unit_tests/test_solver_utils.jl"))
    end
    @time @testset "CRCBS.CBSTests" begin
        include(joinpath(testdir, "unit_tests/test_cbs.jl"))
    end
    @time @testset "CRCBS.CBSTests" begin
        include(joinpath(testdir, "unit_tests/test_cbs.jl"))
    end
    @time @testset "CRCBS.MultiStageCBS.ests" begin
        include(joinpath(testdir, "unit_tests/test_multi_stage_cbs.jl"))
    end
    @time @testset "CRCBS.MetaAgentCBSTests" begin
        include(joinpath(testdir, "unit_tests/test_meta_agent_cbs.jl"))
    end
    # @time @testset "CRCBS.FlowProblemsTests" begin
    #     include(joinpath(testdir, "unit_tests/test_flow_problems.jl"))
    # end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   J2 osculating orbit propagator algorithm.
#
#       This algorithm propagates the orbit considering the secular and
#       short-period perturbations introduced by the J2 gravitational term. The
#       algorithm is based on Kwok version as indicated in [1, p. 708-710].
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export j2osc_init, j2osc!

################################################################################
#                                  Functions
################################################################################

"""
    j2osc_init(epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number, ddn_o6::Number; j2_gc::J2_GravCte{T} = j2_gc_egm08) where T

Initialize the data structure of J2 osculating orbit propagator algorithm.

# Args

* `epoch`: Epoch of the orbital elements [Julian Day].
* `a_0`: Initial semi-major axis [m].
* `e_0`: Initial eccentricity.
* `i_0`: Initial inclination [rad].
* `Ω_0`: Initial right ascension of the ascending node [rad].
* `ω_0`: Initial argument of perigee [rad].
* `f_0`: Initial true anomaly [rad].
* `dn_o2`: First time derivative of the mean motion divided by two [rad/s^2].
* `ddn_o6`: Second time derivative of the mean motion divided by six [rad/s^3].

# Keywords

* `j2_gc`: J2 orbit propagator gravitational constants (see `J2_GravCte`).
           (**Default** = `j2_gc_egm08`)

# Returns

The structure `J2osc_Structure` with the initialized parameters.

# Remarks

The inputs are the mean orbital elements.

"""
function j2osc_init(epoch::Number, a_0::Number, e_0::Number, i_0::Number,
                    Ω_0::Number, ω_0::Number, f_0::Number, dn_o2::Number,
                    ddn_o6::Number; j2_gc::J2_GravCte{T} = j2_gc_egm08) where T

    j2d = j2_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0, dn_o2, ddn_o6;
                  j2_gc = j2_gc)

    # Initialize the structure.
    j2oscd = J2osc_Structure{T}(j2d, epoch, 0, 0, 0, 0, 0, 0, 0)

    # Call the propagation one first time to update the osculating elements.
    j2osc!(j2oscd, 0)

    return j2oscd
end

"""
    j2osc!(j2d::J2osc_Structure{T}, t::Number) where T

Propagate the orbit defined in `j2oscd` (see `J2osc_Structure`) until the time
`t` [s]. Notice that the values in `j2oscd` will be modified.

# Returns

* The position vector represented in the inertial frame at time `t` [m].
* The velocity vector represented in the inertial frame at time `t` [m/s]

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. Notice, however, that the
perturbation theory requires an inertial frame with true equator.

"""
function j2osc!(j2oscd::J2osc_Structure{T}, t::Number) where T
    # First, we need to propagate the mean elements since they are necessary to
    # compute the short-periodic perturbations.
    j2d = j2oscd.j2d
    j2!(j2d, t)

    @unpack_J2_GravCte j2d.j2_gc

    # Obtain the mean elements at this time instant.
    al_k = j2d.al_k
    e_k  = j2d.e_k
    i_k  = j2d.i_k
    Ω_k  = j2d.Ω_k
    ω_k  = j2d.ω_k
    M_k  = j2d.M_k
    f_k  = j2d.f_k
    p_k  = al_k*R0*(1 - e_k^2)
    u_k  = ω_k + f_k

    # Auxiliary variables to reduce the computational burden.
    KJ2 = J2*R0*R0
    sin_i_k, cos_i_k = sincos(i_k)
    sin_f_k, cos_f_k = sincos(f_k)
    sin_2u_k, cos_2u_k = sincos(2u_k)
    sin_2ω_f_k, cos_2ω_f_k = sincos(2ω_k + f_k)
    sin_2ω_3f_k, cos_2ω_3f_k = sincos(2ω_k + 3f_k)
    sin_i_k² = sin_i_k * sin_i_k
    cos_i_k² = cos_i_k * cos_i_k
    p_k² = p_k * p_k
    e_k² = e_k * e_k
    e_cos_f_k = e_k * cos_f_k
    e_sin_f_k = e_k * sin_f_k

    aux1 = 3cos(2u_k) + 3e_k*cos_2ω_f_k + e_k*cos_2ω_3f_k
    aux2 = sqrt(1 - e_k²)
    aux3 = 3cos_i_k² - 1

    # Compute the short-periodic perturbations considering only the J2
    # gravitational term.
    δisp_k = +KJ2*sin_i_k*cos_i_k/(4p_k²)*aux1

    δpsp_k = +KJ2*sin_i_k²/(2p_k)*aux1

    δΩsp_k = -KJ2*cos_i_k/(4p_k²)*(
        6*(f_k - M_k + e_sin_f_k) - 2sin_2u_k - 3e_k*sin_2ω_f_k - e_k*sin_2ω_3f_k
    )

    δrsp_k = -KJ2/(4p_k)*(
        aux3*( 2*aux2/(1 + e_cos_f_k) + e_cos_f_k/(1 + aux2) ) -
        sin_i_k²*cos_2u_k
    )

    δṙsp_k = +KJ2*sqrt(μm)/(4p_k^2.5)*(
        aux3*e_sin_f_k*( aux2 + (1 + e_cos_f_k)^2/(1 + aux2) ) -
        sin_i_k²*(1 - e_cos_f_k)^2*sin_2u_k
    )

    δusp_k = +KJ2/(8p_k²)*(
        (6 - 30cos_i_k²)*(f_k - M_k) +
        4e_k*sin_f_k*((1 - 6cos_i_k²) - aux3/(1 + aux2)) -
        aux3/(1 + aux2)*e_k²*sin(2f_k) +
        (5cos_i_k² - 2)*2e_k*sin_2ω_f_k +
        (7cos_i_k² - 1)*sin_2u_k + 2*cos_i_k²*e_k*sin_2ω_3f_k
    )

    r_k = p_k/(1 + e_k*cos(f_k))
    ṙ_k = sqrt(μm/p_k)*e_k*cos(f_k)

    r_osc_k = r_k + δrsp_k
    ṙ_osc_k = ṙ_k + δṙsp_k
    p_osc_k = p_k + δpsp_k

    A_k = p_osc_k/r_osc_k - 1
    B_k = sqrt(p_osc_k/μm)*ṙ_osc_k

    e_osc_k = sqrt(A_k^2 + B_k^2)
    a_osc_k = p_osc_k/(1-e_osc_k^2)
    i_osc_k = i_k + δisp_k
    Ω_osc_k = Ω_k + δΩsp_k
    u_osc_k = u_k + δusp_k
    f_osc_k = atan(B_k, A_k)
    ω_osc_k = u_osc_k - f_osc_k
    M_osc_k = f_to_M(e_osc_k, f_osc_k)

    # Compute the position and velocity considering the osculating elements.
    r_i_k, v_i_k = kepler_to_rv(a_osc_k, e_osc_k, i_osc_k, Ω_osc_k, ω_osc_k, f_osc_k)

    # Update the J2 orbit propagator structure.
    j2oscd.Δt  = t
    j2oscd.a_k = a_osc_k
    j2oscd.e_k = e_osc_k
    j2oscd.i_k = i_osc_k
    j2oscd.Ω_k = Ω_osc_k
    j2oscd.ω_k = ω_osc_k
    j2oscd.M_k = M_osc_k
    j2oscd.f_k = f_osc_k

    return r_i_k, v_i_k
end


text:
_big(x) = x
_big(x::Number) = big(x)
_big(x::String) = parse(BigFloat, x)

function _big(x::Expr)
    y = x
    y.args .= _big.(y.args)
    return  y
end

macro big(x)
    return esc(_big(x))
end

"Legendre polynomial P_s(x) of degree s defined on the interval [-1..+1]."
function _legendre(j::Int, x::T) where {T}
    if j <= 0
        return one(T)
    elseif j == 1
        return x
    else
        return ( (2j-1) * _legendre(j-1, x) * x - (j-1) * _legendre(j-2, x) ) / j
    end
end

"Legendre polynomial of degree s shifted to the interval [0..1], i.e., P_s(2x-1)."
function _shifted_legendre(s, T=BigFloat)
    _legendre(s, Polynomial(T[-1, 2]))
end

function istriustrict(A::AbstractMatrix)
    m, n = size(A)
    if m == n == 1
        if A[1,1] ≠ 0
            return false
        end
    else
        @inbounds for j in 1:min(n,m-1), i in j:m
            if A[i,j] ≠ 0
                return false
            end
        end
    end
    return true
end

function istrilstrict(A::AbstractMatrix)
    m, n = size(A)
    if m == n == 1
        if A[1,1] ≠ 0
            return false
        end
    else
        @inbounds for j in 2:n, i in 1:min(j,m)
            if A[i,j] ≠ 0
                return false
            end
        end
    end
    return true
end


text:
# Test framework for BondLenghtDistribution

module TestBondLengthDistribution

using Test

export test_single_distribution

include("bond_length_distribution.jl")

# using BondLengthDistributions

function triangular_distribution(min_dist::Real,
                                 dist_max_value::Real,
                                 max_dist::Real)  # ::Tuple{Vector{Float32}, Vector{Float32}}
  x_resolution = 1000

  population = zeros(Float32, x_resolution)
  x_extent = max_dist - min_dist
  peak_index = round(Int32, (dist_max_value - min_dist) / x_extent * x_resolution) + 1
  dy = 1.0 / peak_index
  for i in 1:peak_index
    population[i] = i * dy
  end

  dy = 1.0 / (x_resolution - peak_index)
  for i in peak_index:x_resolution
    population[i] = (1.0 - (i - peak_index) * dy)
  end

  dx = x_extent / x_resolution
  distances = collect(min_dist:dx:max_dist)

  return (distances, population)
end

function test_single_distribution()::Bool
  min_dist = 1.0
  mode_point = 1.4
  max_dist = 2.0
  (lengths, counts) = triangular_distribution(min_dist, mode_point, max_dist)
  distribution = BondLengthDistributions.BondLengthDistribution()

  @test BondLengthDistributions.from_arrays!(lengths, counts, distribution)
  @test BondLengthDistributions.pdf(distribution, 0.99 * min_dist) == 0.0
  @test BondLengthDistributions.pdf(distribution, 1.001 * max_dist) == 0.0

  @test BondLengthDistributions.pdf(distribution, 1.001 * min_dist) > 0.0
  @test BondLengthDistributions.pdf(distribution, 0.99  * max_dist) > 0.0

  maxval = BondLengthDistributions.pdf(distribution, mode_point)
  @test BondLengthDistributions.pdf(distribution, 0.98 * mode_point) < maxval
# @test BondLengthDistributions.pdf(distribution, 1.02 * mode_point) < maxval
  true
end

end  # module TestBondLengthDistribution

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
module Example200530d

#greet() = print("Hello world!")

using ForwardDiff

include("extra_file.jl")

export my_f, derivative_of_my_f

end


text:
__precompile__()
module julia_package_template

# Package code goes here.

end