text:
#     These functions deal with unspecified websocket protocols.
#      Included in server_functions.jl

function ws_general(ws::WebSockets.WebSocket)
    id = cid = ""
    while isopen(ws)
        id, msg = wsmsg_listen(ws)
        cid = "ws_general     $(id)\t"
        # RECEIVED_WS_MSGS[id] = String[]
        # RECEIVED_WS_MSGS_LENGTH[id] = Vector{Int64}()
        # RECEIVED_WS_MSGS_TIME[id] = Vector{Float64}()
        # RECEIVED_WS_MSGS_ALLOCATED[id] = Vector{Int64}()
        sendmsg = "Welcome, general websocket. Our ref.: $(id)"
        clog(cid,  " Send: $sendmsg\n")
        clog(cid, "Received: \n", :yellow, :bold, msg, "\n")
        if startswith(msg, "Ping me!")
            send_ping(ws)
            sendmsg = "Your browser just got pinged (and presumably ponged back). Our ref.: $(id)"
            clog(cid,  " Sending: $sendmsg\n")
            writeto(id, sendmsg)
        elseif startswith(msg, "ws1 echo: Your browser just got pinged ")
                sendmsg = "Close, wait, and navigate to: browsertest2.html"
            clog(cid,  " Sending: $sendmsg\n")
            writeto(id, sendmsg)
        end
    end
    clog(cid, "Websocket $(id) was closed; exits handler.\n")
    nothing
end

function writeto(wsid::String, message)
    if haskey(WEBSOCKETS, wsid)
        write(WEBSOCKETS[wsid], message)
    end
    return nothing
end

"""
Listens for the next websocket message
"""
function wsmsg_listen(ws::WebSockets.WebSocket)
    id = msg = cid = ""
    try
        # Holding here.
        data , t, allocbytes = @timed JSON.parse(String(read(ws)))
        id = data["id"]
        msg = data["msg"]
        cid = "wsmsg_listen     $(id)\t"
        clog(cid,  " Listening...\n")

        if !haskey(WEBSOCKETS,id)
            WEBSOCKETS[id] = ws
            WEBSOCKETS_SUBPROTOCOL[id] = ws
            RECEIVED_WS_MSGS[id] = Vector{String}()
            RECEIVED_WS_MSGS_LENGTH[id] = Vector{Int64}()
            RECEIVED_WS_MSGS_TIME[id] = Vector{Float64}()
            RECEIVED_WS_MSGS_ALLOCATED[id] = Vector{Int64}()
        end
        push!(RECEIVED_WS_MSGS[id],msg)
        push!(RECEIVED_WS_MSGS_LENGTH[id], length(msg))
        push!(RECEIVED_WS_MSGS_TIME[id], t)
        push!(RECEIVED_WS_MSGS_ALLOCATED[id], allocbytes)
    catch e
        if typeof(e) == WebSockets.WebSocketClosedError
            clog(cid, :green, " Websocket was or is being closed.\n")
        else
            clog(cid, "Caught exception..\n", "\t\t", :red, e, "\n")
            if isopen(ws)
                clog(cid, "Closing\n")
                close(ws)
            end
        end
    end
    return id, msg
end
return nothing


text:
# Autogenerated wrapper script for CoinUtils_jll for i686-w64-mingw32-cxx03
export libCoinUtils

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CoinUtils")
JLLWrappers.@declare_library_product(libCoinUtils, "libCoinUtils-3.dll")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libCoinUtils,
        "bin\\libCoinUtils-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
export ssor, sorPrec!, ssorPrecTrans!

"""
x,flag,err,iter,resvec = ssor(A::SparseMatrixCSC, b::Vector; omega::Real=2/3, tol::Real=1e-2, maxIter::Int=1,out::Int=0)

Symmetric successive over-relaxation applied to the linear system A*x = b.

!! Note that upper/lower triangular matrix are never build to save memory !!

Input:

	A       - sparse matrix CSC
	b       - right hand side vector
	omega   - relaxation parameter (omega=1.0: Gauss Seidel)
	tol     - error tolerance (default=1e-2)
	maxIter - maximum number of iterations (default=1)
	out     - flag for output (-1 : no output, 0 : only errors, 1 : final status, 2: error at each iteration)

Output:

	x       - solution
	flag    - exit flag (  0 : desired tolerance achieved,
	                      -1 : maxIter reached without converging)
	err     - error, i.e., norm(A*x-b)/norm(b)
	iter    - number of iterations
	resvec  - error at each iteration
"""
function ssor(A::SparseMatrixCSC,b::Vector;x::Vector=[],tol::Real=1e-2,maxIter::Int=1,omega::Real=1.0,out::Int=0,storeInterm::Bool=false)
n      = length(b)
OmInvD = Array(omega./diag(A))

if storeInterm
    X = zeros(n,maxIter)	# allocate space for intermediates
end

if isempty(x)
	x = zeros(eltype(A),n)
	r = copy(b)
	rnorm0 = norm(b)
else
	r = b - A*x
	rnorm0 = norm(r)
end
if out==2
	println(@sprintf("=== sor ===\n%4s\t%7s\n","iter","relres"))
end

resvec = zeros(maxIter)
flag   = -1
iter   = 1
while iter <= maxIter
	ssorPrec!(A,x,r,OmInvD)
	if storeInterm; X[:,iter] = x; end
	resvec[iter] = norm(r)/rnorm0
	if out==2; println(@sprintf("%3d\t%1.2e",iter,resvec[iter]));end
	if resvec[iter] < tol
		flag=0; break
	end
	iter += 1
end
iter = min(iter,maxIter)

if (flag==-1) && (out>=0)
	println(@sprintf("ssor iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",maxIter,resvec[iter],tol))
elseif  (out>=1)
	println(@sprintf("ssor achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,resvec[iter]))
end
if storeInterm
	return X[:,1:iter],flag,resvec[iter],iter,resvec[1:iter]
else
	return x,flag,resvec[iter],iter,resvec[1:iter]
end
end

"""
 ssorPrec!(A::SparseMatrixCSC,x::Vector,r::Vector,OmInvD::Vector=1./diag(A))

 Applies one SSOR step to A*x = r and updates the residual.

 Input:

	A       - sparse matrix CSC
	x       - current iterate
	r       - residual
	OmInvD  - weighted diagional, that is, omega./diag(A)

 Output: This method changes x and r on the fly.

 see also sor, which uses ssorPrec! in each iteration.
 WARNING: DO NOT USE AS PRECONDITIONER, SINCE THIS METHOD OVEWRITES THE VECTOR r
"""
function ssorPrec!(A::SparseMatrixCSC,x::Vector,r::Vector,OmInvD::Vector=Array(1 ./diag(A)))

dx = 0.0
n  = length(x)

for j=n:-1:1
	dx = OmInvD[j]
	dx *= r[j]
	for gidx = A.colptr[j] : A.colptr[j+1]-1
		r[A.rowval[gidx]] -= A.nzval[gidx]*dx
	end
	x[j] += dx
end
for j=1:n
	dx  = OmInvD[j]
	dx *= r[j]
	for gidx = A.colptr[j] : A.colptr[j+1]-1
		r[A.rowval[gidx]] -= A.nzval[gidx]*dx
	end
	x[j] += dx
end
return x
end

"""
x = ssorPrecTrans!(A::SparseMatrixCSC,X::Array,B::Array,OmInvD::Vector=1 ./diag(A))

Applies SSOR steps to the linear systems A*X = B, for symmetric A. For efficiency, it works on the
transpose of the matrix A.

Input:

  A       - sparse matrix CSC, assumed to be symmetric here
  X       - current iterates
  R       - residuals
  OmInvD  - weighted diagional, that is, omega./diag(A)

Output: This method changes X and on the fly. The Array B is unchanged.

This version is well-suited as a preconditioner, for example, for CG.

Example for solving A*X=B, when A is Sparse Symmtric Positive definite:
    omega = 1.2;
    d = omega./diag(A);
    x = zeros(length(d)) # pre allocation for the preconditioner result.
    PC(r) = (x[:]=0.0; return ssorPrecTrans!(A,x,r,d));
    y = KrylovMethods.cg(A,b,tol=1e-12,maxIter=200,M=PC,out=1)[1]
"""
function ssorPrecTrans!(A::SparseMatrixCSC,X::Array,B::Array,OmInvD::Vector=Array(1 ./diag(A)))

r    = 0.0
n    = size(X,1)
nrhs = size(X,2)

for k=1:nrhs
	for j=n:-1:1
		r = B[j,k]
		for gidx = A.colptr[j] : A.colptr[j+1]-1
			r -= A.nzval[gidx]*X[A.rowval[gidx],k]
		end
		r *= OmInvD[j]
		X[j,k] += r
	end
	for j=1:n
		r = B[j,k]
		for gidx = A.colptr[j] : A.colptr[j+1]-1
			r -= A.nzval[gidx]*X[A.rowval[gidx],k]
		end
		r *= OmInvD[j]
		X[j,k] += r
	end
end
return X
end


text:
module ICGPracticalSixSDW

#prints "Hello World" when called
function printHW()
	println("Hello World")
end

#prints "My name is Simon"
function printSW()
	println("My name is Simon")
end

export printHW
export printSW
end


text:
import PredictMDFull

PredictMDFull.import_all()
PredictMDFull.import_all(Main)

PredictMDFull._import_all()
PredictMDFull._import_all(Main)


text:
#=
Associated YT tutorial: https://youtu.be/TEWy9vZcxW4
=#

using Random

Random.seed!(0)

X = [ 1    2    3    2.5;
      2.0  5.0 -1.0  2.0;
     -1.5  2.7  3.3 -0.8]

# Julia doesn't have classes, so instead we make a struct.
mutable struct Layer_Dense
    weights::Array{Float64, 2}
    biases::Array{Float64, 2}
    output::Array{Float64, 2}
end

# Make a constructor for this struct.
function Layer_Dense(n_inputs, n_neurons)
    weights = 0.10 * randn(n_inputs, n_neurons)
    biases = zeros(1, n_neurons)

    # Initialize with an undefined output.
    output = Array{Float64, 2}(undef, n_inputs, n_neurons)

    Layer_Dense(weights, biases, output)
end

# And define the forward function for this struct.
function forward(layer::Layer_Dense, inputs)
    layer.output = inputs * layer.weights .+ layer.biases
end

# Now do the example.
layer1 = Layer_Dense(4,5)
layer2 = Layer_Dense(5,2)

forward(layer1, X)
forward(layer2, layer1.output)
println(layer2.output)

max_stars_repo_name: ThendBend/NNfSiX

max_issues_repo_name: ThendBend/NNfSiX

max_forks_repo_name: ThendBend/NNfSiX


text:
module MaxHelpingHandCustomStarfield

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/CustomStarfield" CustomStarfield(
    only::String="*", exclude::String="",
    paths::String="starfield", colors::String="ffffff", alphas::String="1",
    shuffle::Bool=true,
    speed::Number=1.0)

placements = CustomStarfield

function Ahorn.canFgBg(effect::CustomStarfield)
    return true, true
end

end


text:
struct BB{T1} <: QuasiNewton{T1}
    approx::T1
end
BB(; inverse = false) = BB(inverse ? Inverse() : Direct())
summary(::BB) = "Barzilai-Borwein"
hasprecon(::BB) = NoPrecon()
update!(scheme::BB{<:Direct}, B, s, y) = _bb(scheme, B, s, y)
update(scheme::BB{<:Direct}, B, s, y) = _bb(scheme, B, s, y)
_bb(scheme::BB{<:Direct}, M, s, y) = dot(s, y) / dot(s, s)

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
module TestModel

using ..NESSie
using ..NESSie: ε0, ec, cos, potprefactor
using Distances: euclidean
using LinearAlgebra: norm, rmul!
using SpecialFunctions: besseli, besselk

#=
    Born models
=#
include("testmodel/born/model.jl")
export BornIon, bornion

include("testmodel/born/potentials.jl")
export φΩ, φΣ

#=
    Multi-charge spheres
=#
include("testmodel/xie/model.jl")
export XieModel

include("testmodel/xie/common.jl")

include("testmodel/xie/nonlocal1.jl")
export NonlocalXieModel1

include("testmodel/xie/potentials.jl")

end # module


text:
sv_constants = GPSL1Constants()

sv1_data = GPSData(
    integrity_status_flag= false,
    TOW = 34810,
    alert_flag = false,
    anti_spoof_flag = true,
    trans_week = 67,
    codeonl2 = 1,
    ura = 2.0,
    svhealth = "000000",
    IODC = "0000011011",
    l2pcode = false,
    T_GD = -1.1175870895385742e-8,
    t_oc = 215984,
    a_f2 = 0.0,
    a_f1 = -8.640199666842818e-12,
    a_f0 = -0.0003322265110909939,
    IODE_Sub_2 = "00011011",
    C_rs = -15.6875,
    Δn = 5.142714214689423e-9,
    M_0 = -1.2713497555403992,
    C_uc = -6.761401891708374e-7,
    e = 0.014137965277768672,
    C_us= 4.753470420837402e-6,
    sqrt_A= 5153.64093208313,
    t_oe = 215984,
    fit_interval = false,
    AODO =13,
    C_ic =  3.129243850708008e-7,
    Ω_0 = -1.0654053987541545,
    C_is = 4.0978193283081055e-8,
    i_0 = 0.9525395470841136,
    C_rc = 276.71875,
    ω = -2.378063276589369,
    Ω_dot = -8.233200088688144e-9,
    IODE_Sub_3 = "00011011",
    IDOT = -1.2286226056345314e-10
)
sv1_decoder_struct = GNSSDecoderState(
    PRN = 7,
    buffer = BitArray(undef, 1502),
    data = sv1_data,
    constants = sv_constants,
    preamble_found =false,
    subframe_count =0,
    prev_30=false,
    prev_29=false,
    data_integrity=true,
    new_data_needed= false,
    subframes_decoded = [1,1,1,1,1],
    nb_prev = 0,
    num_bits_buffered = 187
)
sv1_code_phase = 18126.41231251371
sv1_carrier_phase = -0.019570794887840748
sv1_struct = SatelliteState(decoder_state = sv1_decoder_struct, code_phase = sv1_code_phase, carrier_phase = sv1_carrier_phase)

sv2_data = GPSData(
    integrity_status_flag = false,
    TOW =  34810,
    alert_flag =  false,
    anti_spoof_flag =  true,
    trans_week =  67,
    codeonl2 =  1,
    ura =  2.0,
    svhealth =  "000000",
    IODC =  "0001011001",
    l2pcode =  false,
    T_GD =  5.122274160385132e-9,
    t_oc =  216000,
    a_f2 =  0.0,
    a_f1 =  -1.3642420526593924e-12,
    a_f0 =  -4.171021282672882e-5,
    IODE_Sub_2 =  "01011001",
    C_rs =  148.90625,
    Δn =  4.43518474324698e-9,
    M_0 =  1.74993761786873,
    C_uc =  7.92182981967926e-6,
    e =  0.0053556435741484165,
    C_us =  5.660578608512878e-6,
    sqrt_A =  5153.701900482178,
    t_oe =  216000,
    fit_interval =  false,
    AODO =  31,
    C_ic =  1.0617077350616455e-7,
    Ω_0 =  1.0019742934168678,
    C_is =  8.940696716308594e-8,
    i_0 =  0.9693796725429337,
    C_rc =  272.0,
    ω =  -0.060505953824413546,
    Ω_dot =  -8.216770832915124e-9,
    IODE_Sub_3 =  "01011001",
    IDOT =  2.1858053332800384e-10
)
sv2_decoder_struct = GNSSDecoderState(
    PRN= 8,
    buffer = BitArray(undef, 1502),
    data = sv2_data,
    constants = sv_constants,
    preamble_found = false,
    subframe_count =  0,
    prev_30 =  false,
    prev_29 =  false,
    data_integrity =  true,
    new_data_needed =  false,
    subframes_decoded = [1,1,1,1,1],
    nb_prev =  1,
    num_bits_buffered =  188
)
sv2_code_phase = 1731.7849623335098
sv2_carrier_phase = 0.3726482419297099
sv2_struct = SatelliteState(decoder_state = sv2_decoder_struct, code_phase = sv2_code_phase, carrier_phase = sv2_carrier_phase)

sv3_data = GPSData(
    integrity_status_flag =  false,
    TOW =  34810,
    alert_flag =  false,
    anti_spoof_flag =  true,
    trans_week =  67,
    codeonl2 =  1,
    ura =  2.0,
    svhealth =  "000000",
    IODC =  "0000001001",
    l2pcode =  false,
    T_GD =  2.3283064365386963e-9,
    t_oc =  215984,
    a_f2 =  0.0,
    a_f1 =  -1.057287590811029e-11,
    a_f0 =  -0.00040573813021183014,
    IODE_Sub_2 =  "00001001",
    C_rs =  -128.0,
    Δn =  3.82944522605042e-9,
    M_0 =  -2.5953206787807925,
    C_uc =  -6.614252924919128e-6,
    e =  0.005705653806217015,
    C_us =  1.1416152119636536e-5,
    sqrt_A =  5153.6776695251465,
    t_oe =  215984,
    fit_interval =  false,
    AODO =  11,
    C_ic =  -4.842877388000488e-8,
    Ω_0 =  3.1100423470271776,
    C_is =  -3.166496753692627e-8,
    i_0 =  0.9666236109533345,
    C_rc =  160.65625,
    ω =  -2.6449072365113864,
    Ω_dot =  -7.597102164084918e-9,
    IODE_Sub_3 =  "00001001",
    IDOT =  4.071598169835366e-11
)
sv3_decoder_struct = GNSSDecoderState(
    PRN =  10,
    buffer = BitArray(undef, 1502),
    data = sv3_data,
    constants = sv_constants,
    preamble_found =  false,
    subframe_count =  0,
    prev_30 =  true,
    prev_29 =  true,
    data_integrity =  true,
    new_data_needed =  false,
    subframes_decoded = [1,1,1,1,1],
    nb_prev =  1,
    num_bits_buffered =  188
)
sv3_code_phase = 2975.5753833700155
sv3_carrier_phase = 0.007562681101262569
sv3_struct = SatelliteState(decoder_state = sv3_decoder_struct, code_phase = sv3_code_phase, carrier_phase = sv3_carrier_phase)

sv4_data = GPSData(
    integrity_status_flag =  false,
  TOW =  34810,
  alert_flag =  false,
  anti_spoof_flag =  true,
  trans_week =  67,
  codeonl2 =  1,
  ura =  2.0,
  svhealth =  "000000",
  IODC =  "0000010100",
  l2pcode =  false,
  T_GD =  -1.0710209608078003e-8,
  t_oc =  216000,
  a_f2 =  0.0,
  a_f1 =  2.6147972675971687e-12,
  a_f0 =  -0.00021605053916573524,
  IODE_Sub_2 =  "00010100",
  C_rs =  26.09375,
  Δn =  5.252004481353425e-9,
  M_0 =  0.768362779336754,
  C_uc =  1.4789402484893799e-6,
  e =  0.012468118453398347,
  C_us =  9.216368198394775e-6,
  sqrt_A =  5153.7172775268555,
  t_oe =  216000,
  fit_interval =  false,
  AODO =  28,
  C_ic =  -2.3469328880310059e-7,
  Ω_0 =  -2.2306191090736784,
  C_is =  6.891787052154541e-8,
  i_0 =  0.9282541736265442,
  C_rc =  186.90625,
  ω =  0.8934476147650858,
  Ω_dot =  -8.434994208508933e-9,
  IODE_Sub_3 =  "00010100",
  IDOT =  5.610948004220491e-10
)
sv4_decoder_struct = GNSSDecoderState(
    PRN =  15,
    buffer = BitArray(undef, 1502),
    data = sv4_data,
    constants = sv_constants,
    preamble_found =  false,
    subframe_count =  0,
    prev_30 =  false,
    prev_29 =  false,
    data_integrity =  true,
    new_data_needed =  false,
    subframes_decoded = [1,1,1,1,1],
    nb_prev =  0,
    num_bits_buffered =  187
)
sv4_code_phase = 17790.853982632394
sv4_carrier_phase = 0.2946765711531043
sv4_struct = SatelliteState(decoder_state = sv4_decoder_struct, code_phase = sv4_code_phase, carrier_phase = sv4_carrier_phase)

sv5_data = GPSData(
    integrity_status_flag = false,
    TOW = 34810,
    alert_flag = false,
    anti_spoof_flag = true,
    trans_week = 67,
    codeonl2 = 1,
    ura = 2.0,
    svhealth = "000000",
    IODC = "1011101101",
    l2pcode = false,
    T_GD = -7.916241884231567e-9,
    t_oc = 216000,
    a_f2 = 0.0,
    a_f1 = 9.208633855450898e-12,
    a_f0 = 0.0002515711821615696,
    IODE_Sub_2 = "11101101",
    C_rs = -81.9375,
    Δn = 4.641621913612317e-9,
    M_0 = 0.01827858501071997,
    C_uc = -4.407018423080444e-6,
    e = 0.0007737367413938046,
    C_us = 1.044943928718567e-6,
    sqrt_A = 5153.718004226685,
    t_oe = 216000,
    fit_interval = false,
    AODO = 31,
    C_ic = -2.60770320892334e-8,
    Ω_0 = 2.089468318582403,
    C_is = -2.9802322387695312e-8,
    i_0 = 0.9643063355495126,
    C_rc = 363.78125,
    ω = 2.6017948651652136,
    Ω_dot = -8.454637883889717e-9,
    IODE_Sub_3 = "11101101",
    IDOT = -3.8001582918463414e-10
)
sv5_decoder_struct = GNSSDecoderState(
    PRN =  18,
    buffer = BitArray(undef, 1502),
    data = sv5_data,
    constants = sv_constants,
    preamble_found= false,
    subframe_count= 0,
    prev_30 = false,
    prev_29 = true,
    data_integrity = true,
    new_data_needed = false,
    subframes_decoded = [1,1,1,1,1],
    nb_prev = 0,
    num_bits_buffered = 188
)

sv5_code_phase =  10736.222919350244
sv5_carrier_phase = 0.23463375493884087

sv5_struct = SatelliteState(decoder_state = sv5_decoder_struct, code_phase = sv5_code_phase, carrier_phase = sv5_carrier_phase)

satellite_states = [sv1_struct, sv2_struct, sv3_struct, sv4_struct, sv5_struct]


text:
"""
    Singleton{T} <: Construct{T}

Singleton type empty construct.
"""
struct Singleton{T} <: Construct{T}
    function Singleton{T}() where {T}
        if !Base.issingletontype(T)
            throw(ArgumentError("$T is not a singleton type!"))
        end
        new()
    end
end

Singleton(::Type{T}) where {T} = Singleton{T}()
Singleton(::T) where {T} = Singleton{T}()

deserialize(::Singleton{T}, ::IO; contextkw...) where {T} = T.instance
serialize(::Singleton{T}, ::IO, ::T; contextkw...) where {T} = 0
estimatesize(::Singleton; contextkw...) = 0

Construct(::Type{Nothing}) = Singleton(nothing)
Construct(::Type{Missing}) = Singleton(missing)


text:
using PyPlot, PrintFig

x = linspace(0,2pi,100);
fig = plt.figure();
plot(x,sin(x),color="red");
title("Test plot");
xlabel(L"$x$");
ylabel(L"$\sin(x)$");
printfig(fig)	# Saves to file "FIG1.tex"
printfig(fig,filename="test.tex") # Saves to file "test.tex"


text:
using LibPQ, DataFrames, CSV

host     = "remoteHost"
dbname   = "db1"
user     = "username"
password = "userPwd"
port     = "1111"

conStr= "host=" *host *
        " port=" *port *
        " dbname=" *dbname *
        " user=" *user *
        " password=" *password
conn = LibPQ.Connection(conStr)

df = DataFrame(execute(conn, "SELECT * FROM S1.T1"))
close(conn)

CSV.write("example.csv", df);

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using CxxWrap
using libint_jll

pref = CxxWrap.prefix_path()
path = libint_jll.Libint2_path

run(`rm -rf CMakeCache.txt CMakeFiles lib/libLints.*`)
run(pipeline(`sed "s|JLLPATH|$path|g" _CMakeLists.txt`,`tee CMakeLists.txt`))
run(`cmake -DCMAKE_PREFIX_PATH=$pref .`)
run(`cmake --build .`)


text:
# ----------------------------------------------------
# --- Accessor.jl
# ----------------------------------------------------
# Function to access to radio parameters

getError(obj::UHDBinding) = UHDBindings.getError(obj);
getError(obj::RadioSim) = RadioSims.getError(obj);
getError(obj::SDROverNetwork) = SDROverNetworks.getMD(obj)[3];
getError(obj::RTLSDRBinding) = RTLSDRBindings.getError(obj);

getTimestamp(obj::UHDBinding) = UHDBindings.getTimestamp(obj);
getTimestamp(obj::RadioSim) = RadioSims.getTimestamp(obj);
getTimestamp(obj::SDROverNetwork) = SDROverNetworks.getMD(obj)[1:2];
getTimestamp(obj::RTLSDRBinding) = RTLSDRBindings.getTimestamp(obj);

"""
Get the current sampling rate of the radio device
The second parameter (optionnal) speicfies the Rx or Tx board (default : Rx)
"""
getSamplingRate(obj::AbstractSDR;mode=:rx) = ((mode == :rx) ? obj.rx.samplingRate : obj.tx.samplingRate)
getSamplingRate(obj::PlutoSDR;mode=:rx) = ((mode == :rx) ? obj.rx.effectiveSamplingRate : obj.tx.effectiveSamplingRate)

"""
Get the current carrier frequency   of the radio device
The second parameter (optionnal) speicfies the Rx or Tx board (default : Rx)
"""
getCarrierFreq(obj::AbstractSDR;mode=:rx) = (mode == :rx) ? obj.rx.carrierFreq : obj.tx.carrierFreq
getCarrierFreq(obj::PlutoSDR;mode=:rx) = (mode == :rx) ? obj.rx.effectiveCarrierFreq : obj.tx.effectiveCarrierFreq

"""
Get the current radio gain
The second parameter (optionnal) specifies the Rx or Tx board (default : Rx)
"""
getGain(obj::AbstractSDR;mode=:rx) = (mode == :rx) ? obj.rx.gain : obj.tx.gain
getGain(obj::PlutoSDR;mode=:rx) = AdalmPluto.getGain(obj)

"""
Check if a SDR has already been closed. The falg is true is the SDR ressources have been released and false otherwise.
# --- Syntax
flag = isClosed(radio)
# --- Input parameters
- radio	  : SDR device
# --- Output parameters
- flag : True is SDR is already closed, false otherwise
"""
isClosed(obj::AbstractSDR)  = Bool(obj.tx.released) || Bool(obj.rx.released)
isClosed(obj::PlutoSDR)     = Bool(obj.released)

"""
Returns the radio packet size. Each radio backend encapsulates the IQ samples into chunks of data. The `recv` command can be used with any size but it can be more efficient to match the desired size with the one provided by the radio
# --- Syntax
bufferSize = getBufferSize(radio)
# --- Input parameters
- radio	  : SDR device
# --- Output parameters
bufferSize : Size of radio internal buffer
"""
getBufferSize(obj::AbstractSDR) = obj.rx.packetSize          # We get the fields
getBufferSize(obj::PlutoSDR)    = obj.rx.buf.C_sample_size   # For Pluto this is hidden in the buffer config


text:
# completions
# -----------

using REPL.REPLCompletions
using REPL.REPLCompletions: Completion, KeywordCompletion, PathCompletion, ModuleCompletion,
    PackageCompletion, PropertyCompletion, FieldCompletion, MethodCompletion, BslashCompletion,
    ShellCompletion, DictCompletion

function repl_getcompletions_request(_, params::GetCompletionsRequestParams)
    mod, line = params.mod, params.line
    mod = module_from_string(mod)

    cs = try
        first(Base.invokelatest(completions, line, lastindex(line), mod))
    catch err
        @debug "completion error" exception=(err, catch_backtrace())
        # might error when e.g. type inference fails
        Completion[]
    end
    filter!(is_target_completion, cs)

    return completion.(cs)
end

repl_resolvecompletion_request(conn, completion_item) = completion_item # not used currently, return as is

function is_target_completion(c)
    return c isa PropertyCompletion ||
        c isa FieldCompletion ||
        c isa DictCompletion
end

completion(c) = (
    label = completion_label(c),
    detail = string("REPL completion. ", completion_detail(c)),
    kind = completion_kind(c)
)

completion_label(c) = completion_text(c)
completion_label(c::DictCompletion) = rstrip(completion_text(c), (']', '"'))

completion_detail(c::PropertyCompletion) = begin
    hasproperty(c.value, c.property) || return ""
    t = typeof(getproperty(c.value, c.property))
    string(t)
end
completion_detail(c::FieldCompletion) = string(fieldtype(c.typ, c.field))
completion_detail(c::DictCompletion) = string(valtype(c.dict))
completion_detail(c) = ""

# always mark these completions as Events to distinguish them from LS completions
completion_kind(c) = 22

max_stars_repo_name: JaapHaitsma/julia-vscode

max_issues_repo_name: JaapHaitsma/julia-vscode

max_forks_repo_name: JaapHaitsma/julia-vscode


text:
using robustml
using Documenter

DocMeta.setdocmeta!(robustml, :DocTestSetup, :(using robustml); recursive=true)

makedocs(;
    modules=[robustml],
    authors="Jonas Striaukas, Seyed Abbas Jazaeri",
    repo="https://github.com/sajazaerica/robustml.jl/blob/{commit}{path}#{line}",
    sitename="robustml.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sajazaerica.github.io/robustml.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sajazaerica/robustml.jl",
)

max_stars_repo_name: jstriaukas/robustml

max_issues_repo_name: jstriaukas/robustml

max_forks_repo_name: jstriaukas/robustml


text:
function core_ge(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, M_0, theta, mdoti_c_star, Tti_c_star, Ttj_c_star, Pti_c_star, DTt_des_c_star)

    # Initialize solution
    T = eltype(mdoti_c_star)
    msap_core = zeros(T, (n_t, settings.N_f))
    r_s_star = settings.r_0/sqrt(settings.A_e)
    A_c_star = 1.

    # Turbine transmission loss function
    # Source: Zorumski report 1982 part 2. Chapter 8.2 Equation 3
    g_TT = DTt_des_c_star.^(-4)

    # Calculate acoustic power (Pi_star)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 3
    Pi_star = reshape(8.85e-7 * (mdoti_c_star / A_c_star) .* ((Ttj_c_star .- Tti_c_star) ./ Tti_c_star).^2 .* Pti_c_star.^2 .* g_TT, (n_t,1))

    # Calculate directivity function (D)
    # Take the D function as SAE ARP876E Table 18 and all other values which are not in the table from Zorumski report 1982 part 2. Chapter 8.2 Table II
    array_1 = range(0,180,step=10)
    array_2 = [-0.85, -0.8, -0.75, -0.7, -0.65, -0.6, -0.53, -0.46, -0.39, -0.16, 0.08, 0.31, 0.5, 0.35, 0.12,-0.19,-0.51, -0.8, -0.9]
    #f_D = PCHIPInterpolation.Interpolator(array_1, array_2)
    f_D = LinearInterpolation(array_1, array_2)
    D_funct = reshape(10 .^f_D.(theta), (n_t,1))

    # Calculate the spectral function (S)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 4
    f_p = 400 ./ (1 .- M_0 .* cos.(theta * π / 180.))
    log10ffp = log10.(reshape(data.f, (1, settings.N_f)) ./ reshape(f_p, (n_t,1)))

    # Take the S function as SAE ARP876E Table 17 and all other values which are not in the table from Zorumski report 1982 part 2. Chapter 8.2 Table III
    array_1 = range(-1.1, 1.6, step=0.1)
    array_2 = [-3.87, -3.47, -3.12, -2.72, -2.32, -1.99, -1.7, -1.41, -1.17, -0.97, -0.82, -0.72, -0.82, -0.97, -1.17, -1.41, -1.7, -1.99, -2.32, -2.72, -3.12, -3.47, -3.87, -4.32, -4.72, -5.22, -5.7, -6.2]
    #f_S = PCHIPInterpolation.Interpolator(array_1, array_2)
    f_S = LinearInterpolation(array_1, array_2)
    S_funct = 10. .^f_S.(log10ffp)

    # Calculate mean-square acoustic pressure (msap)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 1
    msap_j = Pi_star * A_c_star / (4 * π * r_s_star^2) .* D_funct .* S_funct ./ reshape((1. .- M_0 .* cos.(π / 180. * theta)).^4, (n_t,1))

    # Multiply with number of engines
    msap_j = msap_j * ac.n_eng

    # Normalize msap by reference pressure
    msap_core = msap_j/settings.p_ref^2

end

function core_pw(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, M_0, theta, mdoti_c_star, Tti_c_star, Ttj_c_star, Pti_c_star, rho_te_c_star, c_te_c_star, rho_ti_c_star, c_ti_c_star)

    # Initialize solution
    T = eltype(mdoti_c_star)
    msap_core = zeros(T, (n_t, settings.N_f))
    r_s_star = settings.r_0/sqrt(settings.A_e)
    A_c_star = 1.

    # Turbine transmission loss function
    # Source: Hultgren, 2012: A comparison of combustor models Equation 6
    zeta = (rho_te_c_star .* c_te_c_star) ./ (rho_ti_c_star .* c_ti_c_star)
    g_TT = 0.8 * zeta ./ (1 .+ zeta).^2

    # Calculate acoustic power (Pi_star)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 3
    Pi_star = reshape(8.85e-7 * (mdoti_c_star / A_c_star) .* ((Ttj_c_star .- Tti_c_star) ./ Tti_c_star).^2 .* Pti_c_star.^2 .* g_TT, (n_t,1))

    # Calculate directivity function (D)
    # Take the D function as SAE ARP876E Table 18 and all other values which are not in the table from Zorumski report 1982 part 2. Chapter 8.2 Table II
    array_1 = range(0,180,step=10)
    array_2 = [-0.85, -0.8, -0.75, -0.7, -0.65, -0.6, -0.53, -0.46, -0.39, -0.16, 0.08, 0.31, 0.5, 0.35, 0.12,-0.19,-0.51, -0.8, -0.9]
    #f_D = PCHIPInterpolation.Interpolator(array_1, array_2)
    f_D = LinearInterpolation(array_1, array_2)
    D_funct = reshape(10 .^f_D.(theta), (n_t,1))

    # Calculate the spectral function (S)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 4
    f_p = 400 ./ (1 .- M_0 .* cos.(theta * π / 180.))
    log10ffp = log10.(reshape(data.f, (1, settings.N_f)) ./ reshape(f_p, (n_t,1)))

    # Take the S function as SAE ARP876E Table 17 and all other values which are not in the table from Zorumski report 1982 part 2. Chapter 8.2 Table III
    array_1 = range(-1.1, 1.6, step=0.1)
    array_2 = [-3.87, -3.47, -3.12, -2.72, -2.32, -1.99, -1.7, -1.41, -1.17, -0.97, -0.82, -0.72, -0.82, -0.97, -1.17, -1.41, -1.7, -1.99, -2.32, -2.72, -3.12, -3.47, -3.87, -4.32, -4.72, -5.22, -5.7, -6.2]
    #f_S = PCHIPInterpolation.Interpolator(array_1, array_2)
    f_S = LinearInterpolation(array_1, array_2)
    S_funct = 10. .^f_S.(log10ffp)

    # Calculate mean-square acoustic pressure (msap)
    # Source Zorumski report 1982 part 2. Chapter 8.2 Equation 1
    msap_j = Pi_star * A_c_star / (4 * π * r_s_star^2) .* D_funct .* S_funct ./ reshape((1. .- M_0 .* cos.(π / 180. * theta)).^4, (n_t,1))

    # Multiply with number of engines
    msap_j = msap_j * ac.n_eng

    # Normalize msap by reference pressure
    msap_core = msap_j/settings.p_ref^2

end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
# Automatically generated test file
function super_exponential2()

g = [
1.000000000000000,
1.064494458917860,
1.284025416687741,
1.755054656960298,
2.718281828459045,
4.770733181967603,
9.487735836358526,
21.380942759123343,
54.598150033144236,
54.598150033144236,
21.380942759123343,
9.487735836358526,
4.770733181967603,
2.718281828459045,
1.755054656960298,
1.284025416687741]

ζ_vec = [1780.25, 409.434, 94.1645, 21.6566, 4.98073, 1.1455, 0.26345, 0.06059, 0.0139349, 0.00320484]

I_sols = [
[1, 16],
[1, 10, 16],
[1, 7, 9, 16],
[1, 7, 10, 11, 16],
[1, 7, 9, 10, 11, 16],
[1, 6, 8, 9, 10, 11, 16],
[1, 6, 8, 9, 10, 11, 12, 16],
[1, 5, 7, 8, 9, 10, 11, 12, 16],
[1, 5, 7, 8, 9, 10, 11, 12, 13, 16],
[1, 4, 6, 7, 8, 9, 10, 11, 12, 13, 16],
[1, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
[1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]]

f_sols = [
4681.821409565800423,
1780.252181436468163,
852.902446264641185,
274.135935240700746,
149.682164816625118,
79.779903677689617,
16.245588913201573,
8.328684761287150,
1.929086519767225,
0.854273036165978,
0.251603511152098,
0.081666565227351,
0.041290097977253,
0.004006052672594,
0.000000000008185]

return g, ζ_vec, I_sols, f_sols

end


text:
abstract type ODEjlAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
# Making the ODE-solver functions into types lets us dispatch on them.
#  Used in the DiffEqBase interface.
struct ode23 <: ODEjlAlgorithm end
struct ode45 <: ODEjlAlgorithm end
struct ode23s <: ODEjlAlgorithm end
struct ode78 <: ODEjlAlgorithm end
struct ode4 <: ODEjlAlgorithm end
struct ode4ms <: ODEjlAlgorithm end
struct ode4s <: ODEjlAlgorithm end


text:
import Distributions
using MatrixCompletion

@testset "$(format("Estimator: MLE[construction]"))"  begin
    @test typeof(MLE{Gaussian}()) == MLE{Gaussian}
    @test typeof(MLE(Gaussian())) == MLE{Gaussian}
    @test typeof(MLE(:Gaussian))  == MLE{Gaussian}
end

@testset "$(format("Estimator: MLE[Gaussian]"))" begin
    for i = 1:10
        input_σ = rand() * 10
        input_μ = rand() * 10
        tc = rand(Distributions.Gaussian(input_μ,input_σ),1000)
        out_1 = Distributions.fit_mle(Distributions.Gaussian,tc)
        out_2 = estimator(MLE{Gaussian}(),tc)
        @test out_1.μ == out_2[:μ] && out_1.σ == out_2[:σ]
    end
end

@testset "$(format("Estimator: MLE[Gamma]"))" begin
    for i = 1:10
        input_α = rand() * 10
        input_θ = rand() * 10
        tc = rand(Distributions.Gamma(input_α,input_θ),10000)
        out_1 = Distributions.fit_mle(Distributions.Gamma,tc)
        out_2 = estimator(MLE{Gamma}(),tc)
        @test out_1.α == out_2[:α] && out_1.θ == out_2[:θ]
    end
end

@testset "$(format("Estimator: MLE[Poisson]"))" begin
    for i = 1:10
        input_λ = rand() * 10
        tc = rand(Distributions.Poisson(input_λ),10000)
        out_1 = sum(tc) / length(tc)
        out_2 = estimator(MLE{Poisson}(),tc)
        @test check(:l2diff,out_1,out_2[:λ]) < 0.5
    end
end

@testset "$(format("Estimator: MLE[Bernoulli]"))" begin
    for i = 1:10
        input_p = rand()
        tc = rand(Distributions.Bernoulli(input_p),10000)
        out_1 = sum(tc) / length(tc)
        out_2 = estimator(MLE{Bernoulli}(),tc)
        @test check(:l2diff,out_1,out_2[:p]) < 0.5
    end
end


text:
"""
```
gensys2_uncertain_altpol(prob_vec::AbstractVector{S}, T_alt::AbstractMatrix{S}, C_alt::AbstractVector{S},
                         T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

gensys2_uncertain_altpol(prob_vec::AbstractVector{S}, T_alt::Vector{Matrix{S}}, C_alt::Vector{Vector{S}},
                         T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

gensys2_uncertain_altpol(prob_vec::Vector{S}, T_alt::Vector{Matrix{S}}, C_alt::Vector{Vector{S}},
                         T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

gensys2_uncertain_altpol(prob_vec::Vector{Vector{S}}, T_alt::Vector{Matrix{S}}, C_alt::Vector{Vector{S}},
                         T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

gensys2_uncertain_altpol(prob_vec::Vector{S}, gensys2_regimes::Vector{Int64}, inds::UnitRange{Int64},
                         systems::Vector{Union{System, RegimeSwitchingSystem}}, is_altpol::Vector{Bool},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

gensys2_uncertain_altpol(prob_vec::Vector{Vector{S}}, gensys2_regimes::Vector{Int64}, inds::UnitRange{Int64},
                         systems::Vector{Union{System, RegimeSwitchingSystem}}, is_altpol::Vector{Bool},
                         Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                         C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}
```

calculates the transition matrices when there is a temporary alternative policy with imperfect awareness.
In particular, for the first four methods
- the `T_alt` and `C_alt` are the matrices associated with the other alternative policies
    agents believe may occur. If T_alt is a vector of matrices and C_alt a vector of vectors,
    then there are multiple alternative policies agents believe may occur.
- the `til` matrices are the predictable form of the gensys matrices under the
    policy the central bank actually implements during the horizon of the temporary policy.
- `T_impl`, `R_impl`, and `C_impl` are the gensys2 matrices of the temporary policy
    that is actually implemented.

As an example, the Federal Reserve in 2020-Q3 switched to flexible AIT. The "alternative" policy is
a Taylor-style rule, and the `til` matrices correspond to the predictable form of flexible AIT.
The `T_impl`, `R_impl`, and `C_impl` include temporary switches to a ZLB specified by `regime_eqcond_info`
followed by a switch to a permanent flexible AIT policy.

Note that it is still assumed that only one policy type occurs during the temporary alternative policy horizon.
For example, we cannot have a temporary ZLB, followed by temporary AIT, before switching back to a Taylor Rule.
To generalize further, we need to make `Γ0_til`, etc., a vector of matrices/vectors and figure out a way
to easily pass these conditions in. However, we can have regime dependence in the temporary rule,
e.g. a ZLB rule that has different ZLB values depending on the regime.

The last two methods allow the user to specify temporary policies as alternative policies
which agents believe may occur but are not implemented.
- `gensys2_regimes`: specifies the correct regimes (and take the form [first_gensys2_regime - 1, first_gensys2_regime, ..., last_gensys2_regime],
    where first_gensys2_regime is the first regime in which a temporary policy applies and
    last_gensys2_regime is the lift-off regime, i.e. the regime in which a permanent policy now applies
    instead of a temporary one).
- `inds`: specifies the required endogenous states for solving the model (i.e. states which are not added by augment_states)
- `systems`: specifies the alternative policies (which may or may not include state space systems
    with temporary policies) agents believe may occur but are not implemented.
- `is_altpol`: a vector specifying which of the alternative policiesare AltPolicy types
    rather than MultiPeriodAltPolicy. Note this vector excludes the implemented policy, so it has the same length of `systems`.
"""
function gensys2_uncertain_altpol(prob_vec::AbstractVector{S}, T_alt::Matrix{S}, C_alt::Vector{S},
                                  T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Rout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Cout  = Vector{Vector{S}}(undef, length(T_impl) + 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    for i in 1:length(T_impl)
        Tbars = prob_vec[1] * T_impl[i] + prob_vec[2] * T_alt # it is assumed T_impl is a vector of the
        Cbars = prob_vec[1] * C_impl[i] + prob_vec[2] * C_alt # T_{t + 1}^{(temporary altpolicy)} matrices

        # hence Tbars will be a vector of T_t matrices
        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    Tout[end] = T_impl[end]
    Rout[end] = R_impl[end]
    Cout[end] = C_impl[end]

    return Tout, Rout, Cout
end

# With time-varying credibility
function gensys2_uncertain_altpol(prob_vec::Vector{Vector{S}}, T_alt::AbstractMatrix{S}, C_alt::AbstractVector{S},
                                  T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Rout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Cout  = Vector{Vector{S}}(undef, length(T_impl) + 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    for i in 1:length(T_impl)
        Tbars = prob_vec[i][1] * T_impl[i] + (1.0-prob_vec[i][1]) * T_alt # it is assumed T_impl is a vector of the
        Cbars = prob_vec[i][1] * C_impl[i] + (1.0-prob_vec[i][1]) * C_alt # T_{t + 1}^{(temporary altpolicy)} matrices

        # hence Tbars will be a vector of T_t matrices
        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    Tout[end] = T_impl[end]
    Rout[end] = R_impl[end]
    Cout[end] = C_impl[end]

    return Tout, Rout, Cout
end

# Same 2 functions as above but where T_alt and C_alt are vectors of matrices/vectors.
function gensys2_uncertain_altpol(prob_vec::AbstractVector{S},
                                  T_alt::Vector{Matrix{S}}, C_alt::Vector{Vector{S}},
                                  T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Rout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Cout  = Vector{Vector{S}}(undef, length(T_impl) + 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    for i in 1:length(T_impl)
        Tbars = prob_vec[1] * T_impl[i] .+ sum([prob_vec[j+1] .* T_alt[j] for j in 1:length(T_alt)]) # it is assumed T_impl is a vector of the
        Cbars = prob_vec[1] * C_impl[i] .+ sum([prob_vec[j+1] .* C_alt[j] for j in 1:length(C_alt)]) # T_{t + 1}^{(temporary altpolicy)} matrices

        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    Tout[end] = T_impl[end]
    Rout[end] = R_impl[end]
    Cout[end] = C_impl[end]

    return Tout, Rout, Cout
end

# With time-varying credibility
function gensys2_uncertain_altpol(prob_vec::Vector{Vector{S}}, T_alt::Vector{Matrix{S}}, C_alt::Vector{Vector{S}},
                                  T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Rout  = Vector{Matrix{S}}(undef, length(T_impl) + 1)
    Cout  = Vector{Vector{S}}(undef, length(T_impl) + 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    for i in 1:length(T_impl)
        Tbars = prob_vec[i][1] .* T_impl[i] .+ sum([prob_vec[i][j+1] .* T_alt[j] for j in 1:length(T_alt)]) # it is assumed T_impl is a vector of the
        Cbars = prob_vec[i][1] * C_impl[i] .+ sum([prob_vec[i][j+1] .* C_alt[j] for j in 1:length(C_alt)])  # T_{t + 1}^{(temporary altpolicy)} matr

        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    Tout[end] = T_impl[end]
    Rout[end] = R_impl[end]
    Cout[end] = C_impl[end]

    return Tout, Rout, Cout
end

# Same 2 functions as above but where T_alt and C_alt replaced by
# systems, a vector of perfectly credible alternative policies, including the implemented one
function gensys2_uncertain_altpol(prob_vec::AbstractVector{S}, gensys2_regimes::UnitRange{Int64}, inds::UnitRange{Int64},
                                  systems::Vector{Union{System, RegimeSwitchingSystem}}, is_altpol::Vector{Bool},
                                  T_impl::Vector{Matrix{S}}, R_impl::Vector{Matrix{S}}, C_impl::Vector{Vector{S}},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(gensys2_regimes) - 1)
    Rout  = Vector{Matrix{S}}(undef, length(gensys2_regimes) - 1)
    Cout  = Vector{Vector{S}}(undef, length(gensys2_regimes) - 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    n_alt = length(prob_vec) - 1
    for (i, reg) in enumerate(gensys2_regimes[2:end - 1])
        # it is assumed T_impl is a vector of the T_{t + 1}^{(temporary altpolicy)} matrices,
        # hence Tbars will be a vector of T_t matrices
        Tbars = prob_vec[1] * (@view systems[1][reg + 1, :TTT][inds, inds]) .+
        sum([prob_vec[j] .* (is_altpol[j - 1] ? (@view systems[j][:TTT][inds, inds]) : (@view systems[j][min(length(systems[j].transitions), reg + 1), :TTT][inds, inds]))
                 for j in 2:n_alt])
        Cbars = prob_vec[1] * (@view systems[1][reg + 1, :CCC][inds]) .+
        sum([prob_vec[j] .* (is_altpol[j - 1] ? (@view systems[j][:CCC][inds]) : (@view systems[j][min(length(systems[j].transitions), reg + 1), :CCC][inds]))
                 for j in 2:n_alt])

        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    l_g2      = last(gensys2_regimes)
    Tout[end] = systems[1][l_g2, :TTT]
    Rout[end] = systems[1][l_g2, :RRR]
    Cout[end] = systems[1][l_g2, :CCC]

    return Tout, Rout, Cout
end

# With time-varying credibility
function gensys2_uncertain_altpol(prob_vec::Vector{Vector{S}}, gensys2_regimes::Vector{Int64}, inds::UnitRange{Int64},
                                  systems::Vector{Union{System, RegimeSwitchingSystem}}, is_altpol::Vector{Bool},
                                  Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                                  C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tout  = Vector{Matrix{S}}(undef, length(gensys2_regimes) - 1)
    Rout  = Vector{Matrix{S}}(undef, length(gensys2_regimes) - 1)
    Cout  = Vector{Vector{S}}(undef, length(gensys2_regimes) - 1)

    # Calculate "uncertain" temporary altpolicy matrices and back out the transition equation
    n_alt = length(prob_vec[1])
    for (i, reg) in enumerate(gensys2_regimes[2:end - 1]) # gensys2_regimes includes 1 extra regime at beginning for boundary condition
        # it is assumed T_impl is a vector of the T_{t + 1}^{(temporary altpolicy)} matrices,
        # hence Tbars will be a vector of T_t matrices
        Tbars = prob_vec[i][1] * (@view systems[1][reg + 1, :TTT][inds, inds]) .+
        sum([prob_vec[i][j] .* (is_altpol[j - 1] ? (@view systems[j][:TTT][inds, inds]) : (@view systems[j][min(length(systems[j].transitions), reg + 1), :TTT][inds, inds]))
                 for j in 2:n_alt])
        Cbars = prob_vec[i][1] * (@view systems[1][reg + 1, :CCC][inds]) .+
        sum([prob_vec[i][j] .* (is_altpol[j - 1] ? (@view systems[j][:CCC][inds]) : (@view systems[j][min(length(systems[j].transitions), reg + 1), :CCC][inds]))
                 for j in 2:n_alt])

        tmp     = Γ2_til * Tbars + Γ0_til
        Tout[i] = tmp \ Γ1_til
        Rout[i] = tmp \ Ψ_til
        Cout[i] = tmp \ (C_til  - Γ2_til * Cbars)
    end

    # Add terminal condition
    l_g2      = last(gensys2_regimes)
    Tout[end] = systems[1][l_g2, :TTT]
    Rout[end] = systems[1][l_g2, :RRR]
    Cout[end] = systems[1][l_g2, :CCC]

    return Tout, Rout, Cout
end


text:
using DashBootstrapComponents

list_group = dbc_listgroup([
    dbc_listgroupitem("Active item", active=true),
    dbc_listgroupitem("Item 2"),
    dbc_listgroupitem("Item 3"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text: include("NonSimpleNumField/Field.jl")


text:
import Base: get, length, show

export CosDict, CosString, CosXString, CosLiteralString, CosNumeric,
    CosBoolean, CosTrue, CosFalse, CosObject, CosNull, CosNullType,
    CosFloat, CosInt, CosArray, CosName, CosDict, CosIndirectObjectRef,
    CosStream, set!, @cn_str, createTreeNode, CosTreeNode, CosIndirectObject,
    CosDictType, IDD, IDDRef, IDDN, IDDNRef, ID, IDRef, IDN, IDNRef, showref

"""
```
    CosObject
```

PDF is a structured document format with lots of internal data
structures like dictionaries, arrays, trees. `CosObject` is the
interface to access these objects and get detailed access to the
objects and gather additional information. Although, defined in the
COS layer, objects of these type are returned from almost all the
APIs. Hence, the objects have a separate significance whether you need
to use the `Cos` layer or not. Below is the object hierarchy.

```
CosObject                           Abstract
    CosNull                         Value (CosNullType)
CosString                           Abstract
CosName                             Concrete
CosNumeric                          Abstract
    CosInt                          Concrete
    CosFloat                        Concrete
CosBoolean                          Concrete
    CosTrue                         Value (CosBoolean)
    CosFalse                        Value (CosBoolean)
CosDict                             Concrete
CosArray                            Concrete
CosStream                           Concrete (always wrapped as an indirect object)
CosIndirectObjectRef                Concrete (only useful when CosDoc is available)
```

*Note*: As a reader API you may not need to instantiate any of `CosObject` types.
They are normally populated as a result of parsing a PDF file.
"""
abstract type CosObject end

"""
```
    get(o::CosObject) -> val
    get(o::CosIndirectObjectRef) -> (objnum, gennum)
```
"""
get(o::CosObject) = o.val

"""
```
    CosString
```
Abstract type that represents a PDF string. In PDF objects are mere byte representations.
They translate to actual text strings by application of fonts and associated encodings.
"""
abstract type CosString <: CosObject end

get(o::T) where {T <: CosString} = copy(o.val)
"""
```
    CosNumeric
```
Abstract type for numeric objects. The objects can be an integer [`CosInt`](@ref) or float
[`CosFloat`](@ref).
"""
abstract type CosNumeric <: CosObject end

"""
```
    CosBoolean
```
A boolean object in PDF which is either a `CosTrue` or `CosFalse`
"""
struct CosBoolean <: CosObject
    val::Bool
end

const CosTrue=CosBoolean(true)
const CosFalse=CosBoolean(false)

struct CosNullType <: CosObject end

"""
```
    CosNull
```
PDF representation of a `null` object. Can be applied to [`CosObject`](@ref)
of any type.
"""
const CosNull=CosNullType()

"""
```
    CosFloat
```
A numeric float data type.
"""
struct CosFloat <: CosNumeric
    val::Float32
end

"""
```
    CosInt
```
An integer in PDF document.
"""
struct CosInt <: CosNumeric
    val::Int
end

"""
```
    CosIndirectObjectRef
```
A parsed data structure to ensure the object information is stored as an object.
This has no meaning without a associated CosDoc. When a reference object is hit
the object should be searched from the CosDoc and returned.
"""
struct CosIndirectObjectRef <: CosObject
    val::Tuple{Int, Int}
    CosIndirectObjectRef(num::Int, gen::Int)=new((num, gen))
end

mutable struct CosIndirectObject{T <: CosObject} <: CosObject
    num::Int
    gen::Int
    obj::T
end

CosIndirectObjectRef(obj::CosIndirectObject) =
    CosIndirectObjectRef(obj.num, obj.gen)

get(o::CosIndirectObject) = get(o.obj)

# Aliases for certain commonly used Union types
# IDD     = Indirect and Direct
# IDDRef  = Indirect, Direct and Reference
# IDDN    = Indirect, Direct, Null
# IDDNRef = Indirect, Direct, Null, Refence

const IDD{X}     = Union{X, CosIndirectObject{X}}
const IDDRef{X}  = Union{X, CosIndirectObject{X}, CosIndirectObjectRef}
const IDDN{X}    = Union{X, CosIndirectObject{X}, CosNullType}
const IDDNRef{X} = Union{X, CosIndirectObject{X}, CosNullType,
                         CosIndirectObjectRef}
const ID{X}      = CosIndirectObject{X}
const IDRef{X}   = Union{CosIndirectObject{X}, CosIndirectObjectRef}
const IDN{X}     = Union{CosIndirectObject{X}, CosNullType}
const IDNRef{X}  = Union{CosIndirectObject{X}, CosNullType, CosIndirectObjectRef}

"""
```
    CosName
```
Name objects are symbols used in PDF documents.
"""
struct CosName <: CosObject
    val::Symbol
    CosName(str::AbstractString) = new(Symbol("CosName_",str))
end

"""
```
    @cn_str(str) -> CosName
```
A string decorator for easier instantiation of a [`CosName`](@ref)

*Example*:
```
julia> cn"Name"
/Name
```
"""
macro cn_str(str)
    return CosName(str)
end

#=
```
    CosXString
```
Concrete representation of a [`CosString`](@ref) object. The underlying data is
represented as hexadecimal characters in ASCII.
=#
struct CosXString <: CosString
    val::Vector{UInt8}
    CosXString(arr::Vector{UInt8})=new(arr)
end

#=
```
    CosLiteralString
```
Concrete representation of a [`CosString`](@ref) object. The underlying data is
represented by byte representations without any encoding.
=#
struct CosLiteralString <: CosString
    val::Vector{UInt8}
    CosLiteralString(arr::Vector{UInt8}) = new(arr)
end

function CosLiteralString(str::AbstractString)
    buf = IOBuffer()
    for c in str
        print(buf, Char(c))
    end
    CosLiteralString(take!(buf))
end

"""
```
    CosArray
```
An array in a PDF file. The objects can be any combination of [`CosObject`](@ref).
"""
mutable struct CosArray <: CosObject
    val::Vector{CosObject}
    CosArray(arr::Vector{CosObject}) = new(arr)
    CosArray() = new(Vector{CosObject}())
end

"""
```
    get(o::CosArray, isNative=false) -> Vector{CosObject}
```
An array in a PDF file. The objects can be any combination of
[`CosObject`](@ref).

`isNative = true` will return the underlying native object inside the `CosArray`
by invoking get method on it.

# Example
```
julia> a = CosArray(CosObject[CosInt(1), CosFloat(2f0), CosInt(3), CosFloat(4f0)])
[1 2.0 3 4.0 ]

julia> get(a)
4-element Array{CosObject,1}:
 1
 2.0
 3
 4.0

julia> get(a, true)
4-element Array{Real,1}:
 1
 2.0f0
 3
 4.0f0
```
"""
get(o::CosArray, isNative=false) = isNative ? map(get, o.val) : o.val

get(o::CosIndirectObject{CosArray}, isNative=false) = get(o.obj, isNative)
"""
```
    length(o::CosArray) -> Int
```
Number of elements in `CosArray`

# Example
```
julia> a = CosArray(CosObject[CosInt(1), CosFloat(2f0),
                              CosInt(3), CosFloat(4f0)])
[1 2.0 3 4.0 ]

julia> length(a)
4
```
"""
length(o::CosArray) = length(o.val)
length(o::ID{CosArray}) = length(o.obj)

Base.getindex(o::CosArray, i::Int) = o.val[i]
Base.getindex(o::ID{CosArray}, i::Int) = getindex(o.obj, i)

"""
```
    CosDict
```
Name value pair of a PDF objects. The object is very similar to the `Dict`
object. The `key` has to be of a [`CosName`](@ref) type.
"""
mutable struct CosDict <: CosObject
    val::Dict{CosName, CosObject}
    CosDict()=new(Dict{CosName, CosObject}())
end

const CosDictType = IDD{CosDict}

"""
```
    get(dict::CosDict, name::CosName, defval::T = CosNull) where T ->
                                                            Union{CosObject, T}
    get(stm::CosStream, name::CosName, defval::T = CosNull) where T ->
                                                            Union{CosObject, T}
```
Returns the value as a [`CosObject`](@ref) for the key `name` or the `defval`
provided.

In case of `CosStream` objects the data is collected from the extent dictionary.

# Example
```
julia> get(catalog, cn"Version")
null

julia> get(catalog, cn"Version", cn"1.4")
/1.4

julia> get(catalog, cn"Version", "1.4")
"1.4"
```
"""
get(dict::CosDict, name::CosName, defval::T = CosNull) where T =
    get(dict.val, name, defval)

"""
```
    set!(dict::CosDict, name::CosName, obj::CosObject) -> CosDict

    set!(stm::CosStream, name::CosName, obj::CosObject) -> CosStream

```
Sets the value on a dictionary object. Setting a `CosNull` object deletes the
object from the dictionary.

In case of `CosStream` objects the data is added to the extent dictionary.

# Example
```
julia> set!(catalog, cn"Version", cn"1.4")

julia> <<
...
/Version /1.4
...
>>
```
"""
function set!(dict::CosDict, name::CosName, obj::CosObject)
    if (obj === CosNull)
        delete!(dict.val, name)
    else
        dict.val[name] = obj
    end
    return dict
end

"""
```
    CosStream
```
A stream object in a PDF. Stream objects have an `extends` disctionary, followed
by binary data.
"""
mutable struct CosStream <: CosObject
    extent::CosDict
    isInternal::Bool
    CosStream(d::CosDict, isInternal::Bool=true) = new(d, isInternal)
end

get(stm::CosStream, name::CosName, defval::T = CosNull) where T =
    get(stm.extent, name, defval)

set!(stm::CosStream, name::CosName, obj::CosObject)=
    (set!(stm.extent, name, obj); stm)

"""
```
    get(stm::CosStream) -> IO
```
Decodes the stream and provides output as an `IO`.

# Example
```
julia> stm

448 0 obj
<<
	/FFilter	/FlateDecode
	/F	(/tmp/tmpIyGPhL/tmp9hwwaG)
	/Length	437
>>
stream
...
endstream
endobj

julia> io = get(stm)
IOBuffer(data=UInt8[...], readable=true, writable=true, ...)

```
"""
get(stm::CosStream) = decode(stm)

mutable struct CosObjectStream <: CosObject
    stm::CosStream
    n::Int
    first::Int
    oids::Vector{Int}
    oloc::Vector{Int}
    populated::Bool
    function CosObjectStream(s::CosStream)
        n = get(s, cn"N")
        @assert n != CosNull
        first = get(s, cn"First")
        @assert first !== CosNull
        n_n = get(n)
        first_n = get(first)
        oids = zeros(Int, n_n)
        oloc = zeros(Int, n_n)
        # cosStreamRemoveFilters(s)
        # read_object_info_from_stm(s, oids, oloc, n_n, first_n)
        new(s, n_n, first_n,oids, oloc, false)
    end
end

get(os::CosObjectStream, name::CosName, defval::T = CosNull) where T =
    get(os.stm, name, defval)

set!(os::CosObjectStream, name::CosName, obj::CosObject)=
    (set!(os.stm, name, obj); os)

get(os::CosObjectStream) = get(os.stm)

get(o::Union{CosIndirectObject{CosDict},
             CosIndirectObject{CosStream},
             CosIndirectObject{CosObjectStream}},
    name::CosName, defval::T = CosNull) where T =
        get(o.obj, name, defval)

set!(os::Union{CosIndirectObject{CosDict},
               CosIndirectObject{CosStream},
               CosIndirectObject{CosObjectStream}},
     name::CosName, obj::CosObject) =
         (set!(os.obj, name, obj); os)

#=
mutable struct CosXRefStream <: CosObject
  stm::CosStream
  isDecoded::Bool
  function CosXRefStream(s::CosStream, isDecoded::Bool=false)
      new(s, isDecoded)
  end
end

get(os::CosXRefStream, name::CosName) = get(os.stm, name)

get(os::CosIndirectObject{CosXRefStream}, name::CosName) = get(os.obj,name)

set!(os::CosXRefStream, name::CosName, obj::CosObject)=
    set!(os.stm, name, obj)

set!(os::CosIndirectObject{CosXRefStream}, name::CosName, obj::CosObject)=
    set!(os.obj,name,obj)

get(os::CosXRefStream) = get(os.stm)
=#

#=
Can be a Number Tree or a Name Tree.

`kids`: is `null` in case of a leaf node
`range`: is `null` in case of a root node
`values`: is `null` in case of an intermediate node

Intent: faster loookup without needing to load the complete tree structure.
Hence, the tree will not be loaded on full scan.
=#
mutable struct CosTreeNode{K <: Union{Int, String}}
    values::Union{Nothing, Vector{Tuple{K, CosObject}}}
    kids::Union{Nothing, Vector{CosIndirectObjectRef}}
    range::Union{Nothing, Tuple{K, K}}
    function CosTreeNode{K}() where {K <: Union{Int, String}}
        new(nothing, nothing, nothing)
    end
end

# If K is Int, it's a number tree else it's a String which is a name tree
function createTreeNode(::Type{K}, dict::IDD{CosDict}) where K
    range = get(dict, cn"Limits")
    kids  = get(dict, cn"Kids")
    node = CosTreeNode{K}()
    if (range !== CosNull)
        r = get(range, true)
        node.range = (K(r[1]), K(r[2]))
    end
    if (kids !== CosNull)
        node.kids = get(kids)
    end
    return populate_values(node, dict)
end

function populate_values(node::CosTreeNode{Int}, dict::IDD{CosDict})
    nums = get(dict, cn"Nums")
    if (nums !== CosNull)
        v = get(nums)
        values = [(get(v[2i-1]), v[2i]) for i=1:div(length(v), 2)]
        node.values = values
    end
    return node
end

function populate_values(node::CosTreeNode{String}, dict::IDD{CosDict})
    names = get(dict, cn"Names")
    if (names !== CosNull)
        v = get(names, false) #true => v[2i] is Tuple{Int} => error
        values = [(String(v[2i-1]), v[2i]) for i=1:div(length(v), 2)]
        node.values = values
    end
    return node
end

# All show methods

show(io::IO, o::CosObject) = print(io, o.val)

showref(io::IO, o::CosObject) = show(io, o)

show(io::IO, o::CosNullType) = print(io, "null")

show(io::IO, o::CosName) = print(io, "/", String(o))

show(io::IO, o::CosXString) =  print(io, "<", String(copy(o.val)), ">")

function show(io::IO, o::CosLiteralString)
    print(io, '(')
    for b in o.val
        c = Char(b)
        if isprint(c)
            print(io, c)
        else
            print(io, '\\')
            print(io, string(b, base=8, pad=3))
        end
    end
    print(io, ')')
end

function show(io::IO, o::CosArray)
  print(io, '[')
  for obj in o.val
    showref(io, obj)
    print(io, ' ')
  end
  print(io, ']')
end

function show(io::IO, o::CosDict)
    print(io, "<<\n")
    for (k,v) in o.val
        print(io, '\t')
        show(io, k)
        print(io, '\t')
        showref(io, v)
        println(io, "")
    end
    print(io, ">>")
end

show(io::IO, stm::CosStream) =
  (show(io, stm.extent); print(io, "\nstream\n...\nendstream"))

show(io::IO, os::CosObjectStream) = show(io, os.stm)

show(io::IO, o::CosIndirectObjectRef) = print(io, o.val[1], ' ', o.val[2], " R")

showref(io::IO, o::CosIndirectObject) = print(io, o.num, ' ', o.gen, " R")

function show(io::IO, o::CosIndirectObject)
    println(io, "")
    println(io, o.num, ' ', o.gen, " obj")
    print(io, o.obj)
    println(io, "\nendobj\n")
end

#=
```
    CosComment
```
A comment object in PDF which is normally ignored as a whitespace.
=#
struct CosComment <: CosObject
    val::String
end

CosComment(barr::Vector{UInt8}) = CosComment(String(Char.(barr)))

show(io::IO, os::CosComment) = print(io, '%', os.val)


text:
#------------------------------------------------------------------------------
# mechanic julia script
# based on mechanical example @ Ramos Jr. & Vieira [ICECE2000]
#
# ricardoaf@lccv.ufal.br (version 21/01/2021)
#------------------------------------------------------------------------------
include("DiscreteModel.jl")

function Mechanical()

	node, elem, α = init(
		connectivity=[1 2; 2 3; 2 3; 2 4; 3 4],
		stiffness=[300, 200, 300, 200, 300],
		displacement=[0, 0, 0, 0],
		force=[0, 0, 0, 10],
		fixnode=[1],
		type="mechanic")

	node, elem = run!(node, elem, α)

	println("  Displacements: ", node.u)
	println("         Forces: ", node.f)
	println("Internal forces: ", elem.f)
end

max_stars_repo_name: ricardoaf/discretemodel

max_issues_repo_name: ricardoaf/discretemodel

max_forks_repo_name: ricardoaf/discretemodel


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

using Markdown
using InteractiveUtils

# ╔═╡ cf365570-054c-11eb-2e15-d117ada9248c
try using AddPackage catch; using Pkg; Pkg.add("AddPackage"); using AddPackage end

# ╔═╡ 1ff5f7d0-0553-11eb-13f4-4fa1c8cfa104
begin
	using Base.Iterators: partition
	@add using CUDA
	@add using Images
	@add using Parameters
	@add using Flux
	using Flux.Data.MNIST
	using Flux: @epochs, mse, throttle
	using Statistics
end

# ╔═╡ 5771c7e0-054c-11eb-3cb4-075289b865c0
md"# Autoencoder for MNIST"

# ╔═╡ 06ab3f50-0554-11eb-1eb0-d536ab6204cf
md"## Training arguments"

# ╔═╡ 57ede040-054d-11eb-1cf9-0b53979535b0
@with_kw mutable struct Args
	α::Float64 = 1e-3      # learning rate
	epochs::Int = 20       # number of epochs
	N::Int = 64            # size of the encoding (i.e. hidden layer)
	batchsize::Int = 1000  # batch size for training
	num::Int = 20          # number of random digits in the sample image (UNUSED)
	throttle::Int = 1      # throttle timeout (called once every X seconds)
end

# ╔═╡ 0f009420-0554-11eb-25c9-4be26cf068ae
md"## Dataset processing"

# ╔═╡ ce88f610-055a-11eb-1d29-b937c6b45dd6
global X = MNIST.images()

# ╔═╡ 91295d80-054d-11eb-2c53-99b31d4390eb
function get_processed_data(args)
	# load images and convert image of type RBG to Float
	imgs = channelview.(X)
	# partition into batches of size `batchsize`
	batches = partition(imgs, args.batchsize)
	traindata = [float(hcat(vec.(imgs)...)) for imgs in batches]
	return gpu.(traindata)
end

# ╔═╡ 158854e0-0554-11eb-15d3-93d4ab118d01
md"## Model training"

# ╔═╡ cad32cf0-054d-11eb-0da8-4bc1c7f2bffc
function train(; kwargs...)
	args = Args(; kwargs...)

	traindata = get_processed_data(args)

	@info "Constructing model..."

	# You can try to make the encoder/decoder network larger
	# Also, the output of encder is a coding of the given input.
	# In this case, the input dimension is 28^2 and the output dimension of the
	# encoder is 32. This implies that the coding is a compressed representation.
	# We can make lossy compression via this `encoder`.
	encoder = Dense(28^2, args.N, leakyrelu) |> gpu
	decoder = Dense(args.N, 28^2, leakyrelu) |> gpu

	# define main model as a Chain of encoder and decoder models
	m = Chain(encoder, decoder)

	@info "Training model..."
	loss(x) = mse(m(x), x)

	# callback, optimizer, and training
	callback = throttle(() -> @show(mean(loss.(traindata))), args.throttle)
	opt = ADAM(args.α)

	@epochs args.epochs Flux.train!(loss, params(m), zip(traindata), opt, cb=callback)

	return m, args, callback
end

# ╔═╡ 1ae3caa0-0554-11eb-372e-edb775d1dc94
md"## Sample from model"

# ╔═╡ 5e5b45c0-054e-11eb-053a-7718f1043426
img(x::Vector) = Gray.(reshape(clamp.(x, 0, 1), 28, 28))

# ╔═╡ 6df36df0-054e-11eb-0b56-0d21d0f36305
function sample(m, num=1)
	# convert image of type RGB to Float
	imgs = channelview.(X)
	# number of random digits (truth)
	before = [imgs[i] for i in rand(1:length(imgs), num)]
	# after applying autoencoder to `before` input image
	after = img.(map(x->cpu(m)(float(vec(x))), before))
	# stack `before` and `after` images them all together
	Gray.(hcat(vcat.(before, after)...))
end

# ╔═╡ a9873310-054e-11eb-3d3e-6d39e2f065d2
m, args, cost = train();

# ╔═╡ 664c8f50-0553-11eb-0951-8b2516f026c1
Markdown.parse(string("\$\$\\frac{1}{m}\\sum_{i=1}^m\\operatorname{Loss}(\\mathbf{\\hat{x}}, \\mathbf{x}) = ", cost(), "\$\$"))

# ╔═╡ 5f088040-054f-11eb-0970-2f1c50b98cb5
sample(m, 20)

# ╔═╡ Cell order:
# ╟─5771c7e0-054c-11eb-3cb4-075289b865c0
# ╠═cf365570-054c-11eb-2e15-d117ada9248c
# ╠═1ff5f7d0-0553-11eb-13f4-4fa1c8cfa104
# ╟─06ab3f50-0554-11eb-1eb0-d536ab6204cf
# ╠═57ede040-054d-11eb-1cf9-0b53979535b0
# ╟─0f009420-0554-11eb-25c9-4be26cf068ae
# ╠═ce88f610-055a-11eb-1d29-b937c6b45dd6
# ╠═91295d80-054d-11eb-2c53-99b31d4390eb
# ╟─158854e0-0554-11eb-15d3-93d4ab118d01
# ╠═cad32cf0-054d-11eb-0da8-4bc1c7f2bffc
# ╟─1ae3caa0-0554-11eb-372e-edb775d1dc94
# ╠═5e5b45c0-054e-11eb-053a-7718f1043426
# ╠═6df36df0-054e-11eb-0b56-0d21d0f36305
# ╠═a9873310-054e-11eb-3d3e-6d39e2f065d2
# ╟─664c8f50-0553-11eb-0951-8b2516f026c1
# ╠═5f088040-054f-11eb-0970-2f1c50b98cb5


text:
using SeqMaker
using OpenGene
using Base.Test

function test_ngs()
    ngs(depth = 0.01)
    return true
end

# write your own tests here
info("Test starts, this may take minutes to hours if it need to download the genome assembly")
@test test_ngs()


text:
"""
    SimTools

Simulation toolbox used in DGE models.
The table of content is:
1. Private Types/Data Structures
2. Samplings
    1. `lhsnorm()` - Ordinary Latin Hypercube Sampling on normal distribution
    2. `lhsmnorm()` - Ordinary Latin Hypercube Sampling on multinvariate normal distribution
    3. `lhsempirical()` - Ordinary Latin Hypercube Sampling on an one-way empirical distribution
3. Algorithms
    1. `cdf_empirical()` - build empirical distribution (cdf) when given samples from a univariate continuous distribution
4. Sample/Testing Data

This module depends on:
1. `Distributions.jl`
2. `Random.jl` (std lib)
3. `LinearAlgebra.jl` (std lib)
4. `StatsBase.jl`
5. `StatsFuns.jl`
6. `Dierckx.jl` - wraper of classical Dierckx interpolation library in Fortran

@author: Clpr @ GitHub
"""
module SimTools
    import StatsFuns, Distributions, StatsBase, Dierckx
    import Random, LinearAlgebra  # std lib

    export lhsnorm, lhsmnorm, lhsempirical
    export cdf_empirical
# =============================

# ------------------
"""
    lhsnorm(N::Int ; μ::Real = 0.0, σ::Real = 1.0)

Ordinary Latin Hypercube Sampling (LHS) on Normal distribution.
Returns a `Vector{Float64}` with `N` elements
"""
function lhsnorm(N::Int ; μ::Real = 0.0, σ::Real = 1.0)
    @assert((N>0)&(σ>0), "invalid N or std")
    local Res::Vector{Float64} = rand(N)
    for x in 1:N
        Res[x] = StatsFuns.norminvcdf( μ, σ, (Res[x] + x - 1)/N )
    end # for
    Random.shuffle!(Res)  # shuffle the sample
    return Res::Vector{Float64}
end # lhsnorm

# ------------------
"""
    lhsmnorm(N::Int ; μ::Vector = [1.0, 1.0], Σ::Matrix = [1.0 0; 0 1.0] )

A wraper of `lhsmnorm(::Int, Distributions.AbstractMvNormal)`.
Receives custom mean vector `μ` and covariance matrix `Σ`.
"""
function lhsmnorm(N::Int ; μ::Vector = [1.0, 1.0], Σ::Matrix = [1.0 0; 0 1.0] )
    return (lhsmnorm(N, Distributions.MvNormal(μ,Σ)))::Matrix{Float64}
end # lhsmnorm
# -------------------
"""
    lhsmnorm(N::Int, D::Distributions.AbstractMvNormal)

Ordinary Latin Hypercube Sampling (LHS) on multinvariate normal distribution.
Returns a `Matrix{Float64}` with `N` rows/samples in which each row of the matrix is a sample.
The dimension of samples depends on the input mean value vector `μ` and covariance matrix `Σ`.

This algorithm follows:
1. Iman, R. L., and W. J. Conover. 1982. A Distribution-free Approach to Inducing Rank Correlation Among Input Variables. Communications in Statistics B 11:311-334
2. Zhang, Y. , & Pinder, G. . (2003). Latin hypercube lattice sample selection strategy for correlated random hydraulic conductivity fields. Water Resources Research, 39(8).
3. My blog: https://clpr.github.io/pages/blogs/190420_LHSonMvNormal.html
"""
function lhsmnorm(N::Int, D::Distributions.AbstractMvNormal)
    # init
    local Ddim::Int = Distributions.dim(D.Σ)
    local Res::Matrix{Float64} = rand(Float64, N, Ddim)
    local Rnk1::Matrix{Int} = zeros(Int, N, Ddim) # (row) ranks of each column for Res with correlation
    local EachDimStd = sqrt.(LinearAlgebra.diag(D.Σ)) # std of each marginal distribution
    # loop to re-scale & inverse to a MvNormal without correlation (each dim is independent)
    # NOTE: just use our lhsnorm() for univariate normal distribution which auto shuffles the results
    for z in 1:Ddim
        Res[:,z] = lhsnorm(N, μ=D.μ[z], σ=EachDimStd[z])
    end # for
    # compute R (Ddim * Ddim size), the correlation matrix (not cov) of the Res without correlation
    # NOTE: when N*P is large, it is a good idea to use a eye(Ddim) to approximate R
    local Rmat::Matrix{Float64} = (N > 1000) ? LinearAlgebra.diagm(0 => ones(Ddim)) : StatsBase.cor(Res)
    # get a new sample matrix Res1
    # NOTE: we do not specially declare Q,P but integrate them to the computing of X1
    # the complier is smart enough to optimize it for least memeory & time costs
    local Res1 = Res * transpose( LinearAlgebra.cholesky( Distributions.cor( D ) ).L * inv(LinearAlgebra.cholesky( Rmat ).L) )
    # record ranks of each column of Res1, then rearrange Res's columns (one by one) according to Res1's col-ranks
    for z in 1:Ddim
        # first, record X1's rank values
        Rnk1[:,z] = StatsBase.ordinalrank(Res1[:,z])
        # then, sort X0's column, the index is the rank value
        Res[:,z] = sort( Res[:,z] )
        # finally, rearrange indices according to Rnk1
        Res[:,z] = Res[Rnk1[:,z], z]
    end # for
    # return (we have already shuffled columns)
    return Res::Matrix{Float64}
end # lhsmnorm

# --------------------
"""
    lhsempirical(N::Int, X::Vector)

Uses Ordinary Latin Hypercube Sampling to sample on a continuous empirical distribution when given a sample vector `X`.
`N` is the number of samples to generate.
Using one-way (cubic) spline interpolation to interpolate the constructed empirical cdf (from sample vector `X`).
Returns a `Vector{Float64}`.

Note: usually, this function is more efficient when using relatively small size of samples to generate large size of samples.
"""
function lhsempirical(N::Int, X::Vector)
    # NOTE: I write the whole proc to estimate the empirical cdf to avoid unnecessary dependency
    # first, rebuild the distribution (empirical cdf)
    local Demp::Tuple = cdf_empirical(X) # this function is defined below
    # then, sampling on Uniform(0,1)
    local Res::Vector{Float64} = rand(N)
    # then, construct a spline interpolation instance
    # NOTE: we use Dierckx.jl, a wraper of Dierckx of Fortran; it is high-performance
    # NOTE: because we will use the inverse cdf, then we use cdf as X but realized values as Y here
    local CdfEmp::Dierckx.Spline1D = Dierckx.Spline1D(Demp[2], Demp[1])
    # then, do ordinary LHS (inverse cdf)
    for x in 1:N
        # rescale Uniform(0,1) to Uniform((x-1)/N, x/N)
        Res[x] = (Res[x] + x - 1)/N
    end # for
    # NOTE: "=" costs similar memory as ".=" but faster
    Res = CdfEmp.(Res)
    Random.shuffle!(Res)  # shuffle the sample
    return Res::Vector{Float64}
end # lhsempirical

# ---------------------
"""
    cdf_empirical(X::Vector)

Rebuilds an empirical cdf when given a sample vector `X`
from a univariate distribution.
Returns a tuple (image) `(X,F)` in which
the 1st member is the input `X` (ascendingly sorted) vector, and `F` is cdf values vector.
"""
function cdf_empirical(X::Vector)
    local Xdim::Int = length(X)
    @assert(Xdim>1, "too few samples (only one!)")
    local Xsort::Vector = sort(X)
    local CdfEmp = Float64[0.0]
    # then, build empirical cdf by estimating empirical quantiles
    # considering that in practice, there will not be so many replicated samples from a continuous distribution,
    # then, using a nested forward-looking searching at every `z` (as follows) is still much more economic than applying conditional filters on `X`
    for z in 2:Xdim
        # this handles the end of loop and the case that Xdim==2
        if z == Xdim
            push!(CdfEmp, 1.0)
            break # out of loop
        end # if
        # else, go on:
        if Xsort[z] != Xsort[z+1] # if there is no same-value in the next ascendingly-ordered sample value
            push!(CdfEmp, z/Xdim) # using formula F(x)= (number of samples which are less than or equal to current x)/N
        else # in this case, we need to consider same-value cases
            local m::Int = 1
            while z+m <= Xdim
                # if no more same-value sample, end loop
                (Xsort[z] == Xsort[z+m]) ? (m += 1) : break
            end # while
            # because we added an extra one in the above while loop, then we need to minus it now
            push!(CdfEmp, (z+m-1)/Xdim)
        end # if
    end # for
    return (Xsort,CdfEmp)::Tuple{Vector{Float64},Vector{Float64}}
end # cdf_empirical

# =============================
end # SimTools

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
function calc_temp(p, v, d, tt, annual_year)
    # for every year, do the same calculations, but then with the new annual_year variables
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    # linear interpolation between the years of analysis, forward-stepping, up to gettime(tt)
    ## interpolation here is also linear for 2015-2020
    if is_first(tt)
        frac = annual_year - 2015
        fraction_timestep = frac / ((gettime(tt)) - 2015)

        v.pt_g_preliminarygmst_ann[yr] = (v.pt_g_preliminarygmst[tt]) * (fraction_timestep) + p.rt_g0_baseglobaltemp * (1 - fraction_timestep)  # difference_to_next_tt * ratio_to_there
    else
        frac = annual_year - gettime(tt - 1)
        fraction_timestep = frac / ((gettime(tt)) - (gettime(tt - 1))) # check if +1 might also need to feature here.

        v.pt_g_preliminarygmst_ann[yr] = (v.pt_g_preliminarygmst[tt] ) * (fraction_timestep) + v.pt_g_preliminarygmst[tt - 1] * (1 - fraction_timestep) # difference_to_next_tt * ratio_to_there
    end
    # Without surface albedo, just equal
    v.rt_g_globaltemperature_ann[yr] = v.pt_g_preliminarygmst_ann[yr]

    # Setting regional temperature
    for r in d.region
        if use_variability
            r_variationtemp = rand(Normal(0.0, max(0., p.rvarsd_regionalvariabilitystandarddeviation[r]))) # NEW - with variation
            # r_variationtemp = r_variationtemp * p.var_multiplier # for sensitivity analysis for variability
        else
            r_variationtemp = 0 # NEW - no variation
        end
        v.rtl_realizedtemperature_ann[yr, r] = v.rt_g_globaltemperature_ann[yr] * p.ampf_amplification[r] + r_variationtemp # NEW - add interannual variability
    end

    # Setting global temperature
    if use_variability
        g_variationtemp = rand(Normal(0.0, max(0., p.gvarsd_globalvariabilitystandarddeviation)))
        # g_variationtemp = g_variationtemp * p.var_multiplier # for sensitivity analysis for variability
    else
        g_variationtemp = 0
    end
    v.rt_g_globaltemperature_ann[yr] = v.pt_g_preliminarygmst_ann[yr] + g_variationtemp

    # Land average temperature
    v.rtl_g_landtemperature_ann[yr] = sum(v.rtl_realizedtemperature_ann[yr, :]' .* p.area') / sum(p.area)

    # Ocean average temperature
    v.rto_g_oceantemperature_ann[yr] = (p.area_e_eartharea * v.rt_g_globaltemperature_ann[yr] - sum(p.area) * v.rtl_g_landtemperature_ann[yr]) / (p.area_e_eartharea - sum(p.area))

end

@defcomp ClimateTemperature begin

    region = Index()
    year = Index()

    # Basic parameters
    area = Parameter(index=[region], unit="km2")
    y_year_0 = Parameter(unit="year")
    y_year = Parameter(index=[time], unit="year")
    y_year_ann = Parameter(index=[year], unit="year")
    area_e_eartharea = Parameter(unit="km2", default=5.1e8)
    use_seaice = Parameter{Bool}()

    # Initial temperature outputs
    rt_g0_baseglobaltemp = Parameter(unit="degreeC", default=0.9461666666666667) # needed for feedback in CO2 cycle component
    rtl_0_baselandtemp = Variable(index=[region], unit="degreeC")
    rtl_g0_baselandtemp = Variable(unit="degreeC") # needed for feedback in CH4 and N2O cycles

    # variability parameters
    gvarsd_globalvariabilitystandarddeviation = Parameter(unit="degreeC", default=0.11294)  # default from https://github.com/jkikstra/climvar
    rvarsd_regionalvariabilitystandarddeviation = Parameter(index=[region], unit="degreeC") # provided in data folder

    # Rate of change of forcing
    fant_anthroforcing = Parameter(index=[time], unit="W/m2")
    ft_0_totalforcing0 = Parameter(unit="W/m2", default=3.202)
    fsd_g_0_directsulphate0 = Parameter(unit="W/m2", default=-0.46666666666666673)
    fsi_g_0_indirectsulphate0 = Parameter(unit="W/m2", default=-0.17529148701061475)

    # Climate sensitivity calculations
    tcr_transientresponse = Parameter(unit="degreeC", default=1.7666666666666668)
    frt_warminghalflife = Parameter(unit="year", default=28.333333333333332)

    ecs_climatesensitivity = Variable(unit="degreeC")

    # Unadjusted temperature calculations
    fslope_CO2forcingslope = Parameter(unit="W/m2", default=5.5)
    pt_g_preliminarygmst = Variable(index=[time], unit="degreeC")
    pt_g_preliminarygmst_ann = Variable(index=[year], unit="degreeC")

    # Surface albedo parameters
    alb_t_switch = Parameter(default=10.0)
    alb_saf_quadr_mean_t2_coeff = Parameter(default=-0.001955963673212411)
    alb_saf_quadr_mean_t1_coeff = Parameter(default=-0.006407351323331157)
    alb_saf_quadr_mean_t0_coeff = Parameter(default=0.36372124028281877)
    alb_saf_quadr_std = Parameter(default=0.10908771841232413)
    alb_saf_lin_mean = Parameter(default=0.0636051293985324)
    alb_saf_lin_std = Parameter(default=0.02907749205894582)
    alb_emulator_rand = Parameter(default=0)

    # Surface albedo internal variables
    alb_fsaf_y0 = Variable(unit="W/m2")
    alb_saf_y0 = Variable(unit="W/m2/degC")
    alb_fsaf_ecs = Variable(unit="W/m2")
    alb_saf_ecs = Variable(unit="W/m2/degC")
    alb_fsaf_t_switch = Variable(unit="W/m2")

    # Global outputs
    rt_g_globaltemperature = Variable(index=[time], unit="degreeC")
    rt_g_globaltemperature_ann = Variable(index=[year], unit="degreeC")
    rto_g_oceantemperature = Variable(index=[time], unit="degreeC")
    rto_g_oceantemperature_ann = Variable(index=[year], unit="degreeC")
    rtl_g_landtemperature = Variable(index=[time], unit="degreeC")
    rtl_g_landtemperature_ann = Variable(index=[year], unit="degreeC")

    # Regional outputs
    ampf_amplification = Parameter(index=[region])

    rtl_realizedtemperature = Variable(index=[time, region], unit="degreeC")
    rtl_realizedtemperature_ann = Variable(index=[year, region], unit="degreeC")

    function init(p, v, d)
        for rr in d.region
            v.rtl_0_baselandtemp[rr] = p.rt_g0_baseglobaltemp * p.ampf_amplification[rr]
        end

        # Equation 21 from Hope (2006): initial global land temperature
        v.rtl_g0_baselandtemp = sum(v.rtl_0_baselandtemp .* p.area) / sum(p.area)

        # Inclusion of transient climate response from Hope (2009)
        v.ecs_climatesensitivity = p.tcr_transientresponse / (1. - (p.frt_warminghalflife / 70.) * (1. - exp(-70. / p.frt_warminghalflife)))

        ## Surface albedo internal variables
        v.alb_fsaf_y0 = ((p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.rt_g0_baseglobaltemp) + (p.alb_saf_quadr_std * p.rt_g0_baseglobaltemp) * p.alb_emulator_rand)
        v.alb_saf_y0 = (p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t1_coeff * p.rt_g0_baseglobaltemp + p.alb_saf_quadr_mean_t0_coeff) + p.alb_saf_quadr_std * p.alb_emulator_rand

        if v.ecs_climatesensitivity <= p.alb_t_switch
            v.alb_fsaf_ecs = (p.alb_saf_quadr_mean_t2_coeff * (v.ecs_climatesensitivity^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.ecs_climatesensitivity^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.ecs_climatesensitivity) + (p.alb_saf_quadr_std * v.ecs_climatesensitivity) * p.alb_emulator_rand
        else
            v.alb_fsaf_ecs = (p.alb_saf_quadr_mean_t2_coeff * ( p.alb_t_switch^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * ( p.alb_t_switch^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.alb_t_switch + p.alb_saf_lin_mean * (v.ecs_climatesensitivity - p.alb_t_switch)) + (p.alb_saf_quadr_std * p.alb_t_switch + p.alb_saf_lin_std * (v.ecs_climatesensitivity - p.alb_t_switch)) * p.alb_emulator_rand
        end
        v.alb_saf_ecs = v.alb_fsaf_ecs / v.ecs_climatesensitivity
        v.alb_fsaf_t_switch = (p.alb_saf_quadr_mean_t2_coeff * (p.alb_t_switch^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.alb_t_switch^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.alb_t_switch) + (p.alb_saf_quadr_std * p.alb_t_switch) * p.alb_emulator_rand
    end

    function run_timestep(p, v, d, tt)
        # Rate of change of grand total forcing
        if is_first(tt)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            rate_fant = 0 # inferred from spreadsheet
            deltat = p.y_year[tt] - p.y_year_0
        elseif tt == TimestepIndex(2)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            rate_fant = (p.fant_anthroforcing[tt - 1] - fant0) / (p.y_year[tt - 1] - p.y_year_0)
            deltat = p.y_year[tt] - p.y_year[tt - 1]
        else
            rate_fant = (p.fant_anthroforcing[tt - 1] - p.fant_anthroforcing[tt - 2]) / (p.y_year[tt - 1] - p.y_year[tt - 2])
            deltat = p.y_year[tt] - p.y_year[tt - 1]
        end

        BB = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * rate_fant
        EXPT = exp(-deltat / p.frt_warminghalflife)

        if is_first(tt)
            fant0 = p.ft_0_totalforcing0 + p.fsd_g_0_directsulphate0 + p.fsi_g_0_indirectsulphate0
            AA = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * fant0
            v.pt_g_preliminarygmst[tt] = p.rt_g0_baseglobaltemp + (AA - p.rt_g0_baseglobaltemp) * (1 - EXPT) # Drop BB because 0
        else
            AA = v.ecs_climatesensitivity / (p.fslope_CO2forcingslope * log(2.0)) * p.fant_anthroforcing[tt - 1]
            v.pt_g_preliminarygmst[tt] = v.pt_g_preliminarygmst[tt - 1] + (AA - p.frt_warminghalflife * BB - v.pt_g_preliminarygmst[tt - 1]) * (1 - EXPT) + deltat * BB
        end

        if !p.use_seaice
            # Without surface albedo, just equal
            v.rt_g_globaltemperature[tt] = v.pt_g_preliminarygmst[tt]
        else
            alb_saf_approx = (((p.alb_saf_quadr_mean_t2_coeff * (v.pt_g_preliminarygmst[tt]^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.pt_g_preliminarygmst[tt]^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.pt_g_preliminarygmst[tt]) + (p.alb_saf_quadr_std * v.pt_g_preliminarygmst[tt]) * p.alb_emulator_rand) - v.alb_fsaf_y0) / (v.pt_g_preliminarygmst[tt] - p.rt_g0_baseglobaltemp)
            alb_saf_adjust = alb_saf_approx - v.alb_saf_ecs
            ecs_alb_adj = v.ecs_climatesensitivity / (1 - v.ecs_climatesensitivity * alb_saf_adjust / (log(2) * p.fslope_CO2forcingslope))
            frt_alb_adj = p.frt_warminghalflife / (1 - v.ecs_climatesensitivity * alb_saf_adjust / (log(2) * p.fslope_CO2forcingslope))
            if is_first(tt)
                alb_fsaf_adjust = ((p.alb_saf_quadr_mean_t2_coeff * (p.rt_g0_baseglobaltemp^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (p.rt_g0_baseglobaltemp^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * p.rt_g0_baseglobaltemp) + (p.alb_saf_quadr_std * p.rt_g0_baseglobaltemp) * p.alb_emulator_rand) - alb_saf_approx * p.rt_g0_baseglobaltemp
                v.rt_g_globaltemperature[tt] = p.rt_g0_baseglobaltemp + (ecs_alb_adj * (fant0 + alb_fsaf_adjust) / (log(2) * p.fslope_CO2forcingslope) - p.rt_g0_baseglobaltemp) * (1 - exp(-deltat / frt_alb_adj))
            else
                if v.rt_g_globaltemperature[tt - 1] <= p.alb_t_switch
                    alb_fsaf_adjust = ((p.alb_saf_quadr_mean_t2_coeff * (v.rt_g_globaltemperature[tt - 1]^3) / 3 + p.alb_saf_quadr_mean_t1_coeff * (v.rt_g_globaltemperature[tt - 1]^2) / 2 + p.alb_saf_quadr_mean_t0_coeff * v.rt_g_globaltemperature[tt - 1]) + (p.alb_saf_quadr_std * v.rt_g_globaltemperature[tt - 1]) * p.alb_emulator_rand) - alb_saf_approx * v.rt_g_globaltemperature[tt - 1]
                else
                    alb_fsaf_adjust = v.alb_fsaf_t_switch + p.alb_saf_lin_mean * (v.rt_g_globaltemperature[tt - 1] - p.alb_t_switch) + (p.alb_saf_lin_std * (v.rt_g_globaltemperature[tt - 1] - p.alb_t_switch)) * p.alb_emulator_rand - alb_saf_approx * v.rt_g_globaltemperature[tt - 1]
                end
                v.rt_g_globaltemperature[tt] = v.rt_g_globaltemperature[tt - 1] + (ecs_alb_adj * (p.fant_anthroforcing[tt - 1] - rate_fant * frt_alb_adj + alb_fsaf_adjust) / (log(2) * p.fslope_CO2forcingslope) - v.rt_g_globaltemperature[tt - 1]) * (1 - exp(-deltat / frt_alb_adj)) + ecs_alb_adj * (rate_fant * deltat) / (log(2) * p.fslope_CO2forcingslope)
            end
        end

        # Adding adjustment, from Hope (2009)
        for rr in d.region
            v.rtl_realizedtemperature[tt, rr] = v.rt_g_globaltemperature[tt] * p.ampf_amplification[rr]
        end

        # Land average temperature
        v.rtl_g_landtemperature[tt] = sum(v.rtl_realizedtemperature[tt, :]' .* p.area') / sum(p.area)

        # Ocean average temperature
        v.rto_g_oceantemperature[tt] = (p.area_e_eartharea * v.rt_g_globaltemperature[tt] - sum(p.area) * v.rtl_g_landtemperature[tt]) / (p.area_e_eartharea - sum(p.area))

        # annual interpolation of the timestep-calculated variables.
        if is_first(tt)
            for annual_year = 2015:gettime(tt) # because p["y_year_0"] does not exist (globally)
                calc_temp(p, v, d, tt, annual_year) # NEW -- introduce function for calulating temperature
            end
        else
            for annual_year = (gettime(tt - 1) + 1):gettime(tt)
                calc_temp(p, v, d, tt, annual_year) # NEW -- use  function for calulating temperature
            end
        end
    end
end

function addclimatetemperature(model::Model, use_seaice::Bool)
    climtemp = add_comp!(model, ClimateTemperature)

    climtemp[:use_seaice] = use_seaice

    return climtemp
end


text:
# [SOREN] TRUST REGION MODIFICATION OF NEWTON'S METHOD
# [N&W] Numerical optimization
# [Yuan] A review of trust region algorithms for optimization
abstract type TRSPSolver end
abstract type NearlyExactTRSP <: TRSPSolver end

include("subproblemsolvers/NWI.jl")
function minimize!(objective, x0, scheme::Newton, B0=nothing, options=OptOptions())
    minimize!(objective, x0, (scheme, NWI()), B0, options)
end
function minimize!(objective, x0, approach::Tuple{<:Any, <:TRSPSolver}, B0=nothing, options=OptOptions())
    tr_minimize!(objective, x0, approach, B0, options)
end
include("tr_optimize_inplace.jl")


text:
function verticalA( w :: Array{Float64,3}, T :: Array{Float64,3}, wintercept :: OnePointR = OnePointR())
  w1 = dimAverage(w,3, wintercept)
  return T.* w1
end

export verticalA


text:
using DataFrames, CSV, Dates, Statistics
data = dropmissing(CSV.read("../data/purchaseData.csv", copycols=true))

data[!,:Date] = Date.(data[!,:Date], "d/m/y")
println(first(sort(data, :Date), 3),"\n")

println(first(filter(row -> row[:Price] > 50000, data),3 ),"\n")

categorical!(data, :Grade)
println(first(data, 3), "\n")

println(
   by(data, :Grade, :Price =>
       x -> ( NumSold=length(x), AvgPrice=mean(x)) )
   )