text:
using Maker
using Base.Test

writecsv("in1.csv", rand(3,3))
writecsv("in2.csv", rand(3,3))

@desc "Combine input files"
Maker.file("x.csv", ["in1.csv", "in2.csv"]) do
    println("Reading input data.")
    x = readcsv("in1.csv") + readcsv("in2.csv")
    println("Writing 'x.csv'.")
    writecsv("x.csv", x)
end

@desc "Load variable `x`"
Maker.variable("x", "x.csv") do
    println("Reading `x`.")
    readcsv("x.csv")
end

function process_x2()
    println("Processing `x`.")
    x2 = copy(x)
    x2 = pi .* x
    println("Writing 'x2.csv'.")
    writecsv("x2.csv", x2)
end

@desc "Make x2 using second input file"
Maker.file(process_x2, "x2.csv", "x")

show(tasks("x2.csv"))

@desc "Read `x2`"
Maker.variable("x2", "x2.csv") do
    println("Reading `x2`.")
    readcsv("x2.csv")
end

Maker.task("default", "x2.csv")

@desc "Load all variables"
Maker.task("vars", ["x", "x2"])

@desc "Delete generated csv files"
Maker.task("clean") do
    println("Deleting generated csv files.")
    Maker.rm("in1.csv")
    Maker.rm("in2.csv")
    Maker.rm("x.csv")
    Maker.rm("x2.csv")
end

make()

make()

@desc "Make x2 using second input file"
Maker.file("x2.csv", "x") do
    println("Processing `x` after redefining the x2.csv task.")
    x2 = copy(x)
    x2 = pi .* x
    println("Writing 'x2.csv'.")
    writecsv("x2.csv", x2)
end

make()

println("== Touch'ing in1.csv...")

run(`touch in1.csv`)

make()

println("== Changing x...")

tmp = x[1,1]
x[1,1] = 3

make()  # This should change `x` back to what it was.

@test tmp == x[1,1]

make("clean")

@test length(tasks()) == 7
@test tasks("clean").name == "clean"
show(tasks())
show(tasks("x.csv"))

println("Done.")


text:
mutable struct Body{T} <: AbstractBody{T}
    id::Int64
    name::String

    m::T
    J::SMatrix{3,3,T,9}

    x::Vector{SVector{3,T}}
    q::Vector{Quaternion{T}}

    F::Vector{SVector{3,T}}
    τ::Vector{SVector{3,T}}

    s0::SVector{6,T}
    s1::SVector{6,T}
    f::SVector{6,T}

    function Body(m::T, J::AbstractArray{T,2}; name::String="") where T
        x = [zeros(T, 3)]
        q = [Quaternion{T}()]

        F = [zeros(T, 3)]
        τ = [zeros(T, 3)]

        s0 = zeros(T, 6)
        s1 = zeros(T, 6)
        f = zeros(T, 6)

        new{T}(getGlobalID(), name, m, J, x, q, F, τ, s0, s1, f)
    end

    function Body(shape::Shape; name::String="")
        body = Body(shape.m, shape.J; name=name)
        push!(shape.bodyids, body.id)
        return body
    end

    # Body(m::T, J::SMatrix{3,3,T,9}, x::Vector{SVector{3,T}}, q::Vector{Quaternion{T}}, F::Vector{SVector{3,T}}, τ::Vector{SVector{3,T}},
    #     s0::SVector{6,T}, s1::SVector{6,T}, f::SVector{6,T}) where T = new{T}(getGlobalID(), m, J, x, q, F, τ, s0, s1, f)
end

mutable struct Origin{T} <: AbstractBody{T}
    id::Int64
    name::String

    Origin{T}(; name::String="") where T = new{T}(getGlobalID(), name)

    function Origin(shape::Shape{T}; name::String="") where T
        origin = Origin{T}(name=name)
        push!(shape.bodyids, origin.id)
        return origin
    end
end

function Base.deepcopy(b::Body)
    contents = []
    for i = 2:10
        push!(contents, deepcopy(getfield(b, i)))
    end

    Body(contents...)
end

function Base.deepcopy(b::Body, shape::Shape)
    contents = []
    for i = 2:10
        push!(contents, deepcopy(getfield(b, i)))
    end

    body = Body(contents...)
    @assert (body.m == shape.m && body.J == shape.J)
    push!(shape.bodyids, body.id)
    return body
end

Base.length(::Body) = 6

function setInit!(body::Body{T};
    x::AbstractVector = zeros(T, 3),
    q::Quaternion = Quaternion{T}(),
    F::AbstractVector = zeros(T, 3),
    τ::AbstractVector = zeros(T, 3)) where T

    body.x[1] = convert(SVector{3,T}, x)
    body.q[1] = q
    body.F[1] = convert(SVector{3,T}, F)
    body.τ[1] = convert(SVector{3,T}, τ)
    return
end

function setInit!(body1::Body{T}, body2::Body{T}, p1::AbstractVector, p2::AbstractVector;
    q::Quaternion = Quaternion{T}(),
    F::AbstractVector = zeros(T, 3),
    τ::AbstractVector = zeros(T, 3)) where T

    p1 = convert(SVector{3,T}, p1)
    p2 = convert(SVector{3,T}, p2)

    x2 = body1.x[1] + vrotate(p1, body1.q[1]) - vrotate(p2, q)
    setInit!(body2; x = x2, q = q, F = F, τ = τ)
end

function setInit!(body1::Origin{T}, body2::Body{T}, p1::AbstractVector, p2::AbstractVector;
    q::Quaternion = Quaternion{T}(),
    F::AbstractVector = zeros(T, 3),
    τ::AbstractVector = zeros(T, 3)) where T

    p2 = convert(SVector{3,T}, p2)

    x2 = p1 - vrotate(p2, q)
    setInit!(body2; x = x2, q = q, F = F, τ = τ)
end

@inline getM(body::Body{T}) where T = [[SMatrix{3,3,T,9}(I*body.m);@SMatrix zeros(T,3,3)] [@SMatrix zeros(T,3,3);body.J]]

@inline getx3(body::Body, Δt) = getvnew(body) * Δt + body.x[2]
@inline getq3(body::Body, Δt) = Quaternion(Lmat(body.q[2]) * ωbar(body, Δt))
#TODO why not simple from s0?
@inline getv1(body::Body, Δt) = (body.x[2] - body.x[1]) / Δt
@inline function getω1(body::Body, Δt)
    q1 = body.q[1]
    q2 = body.q[2]
    2 / Δt * (q1.s * imag(q2) - q2.s * imag(q1) - cross(imag(q1), imag(q2)))
    # 2/Δt*(VLᵀmat(body.q[1])*body.q[2])
end

# TODO use SOneTo and SUnitRange once they are faster
# @inline getvnew(body::Body) = body.s1[SOneTo(3)]
# @inline getωnew(body::Body) = body.s1[SUnitRange(4,6)]
@inline getvnew(body::Body) = body.s1[SVector(1, 2, 3)]
@inline getωnew(body::Body) = body.s1[SVector(4, 5, 6)]

@inline function derivωbar(body::Body{T}, Δt) where T
    ωnew = getωnew(body)
    msq = -sqrt(4 / Δt^2 - dot(ωnew, ωnew))
    Δt / 2 * [ωnew' / msq; SMatrix{3,3,T,9}(I)]
end

@inline function ωbar(body::Body, Δt)
    ωnew = getωnew(body)
    Δt / 2 * Quaternion(sqrt(4 / Δt^2 - dot(ωnew, ωnew)), ωnew)
end

@inline function dynamics(mechanism, body::Body{T}) where T
    No = mechanism.No
    Δt = mechanism.Δt

    ezg = SVector{3,T}(0, 0, -mechanism.g)
    dynT = body.m * ((getvnew(body) - getv1(body, Δt)) / Δt + ezg) - body.F[No]

    J = body.J
    ω1 = getω1(body, Δt)
    ωnew = getωnew(body)
    sq1 = sqrt(4 / Δt^2 - ω1' * ω1)
    sq2 = sqrt(4 / Δt^2 - ωnew' * ωnew)
    dynR = skewplusdiag(ωnew, sq2) * (J * ωnew) - skewplusdiag(ω1, sq1) * (J * ω1) - 2 * body.τ[No]

    body.f = [dynT;dynR]

    for cid in connections(mechanism.graph, body.id)
        GtλTof!(mechanism, body, geteqconstraint(mechanism, cid))
    end

    for cid in ineqchildren(mechanism.graph, body.id)
        NtγTof!(mechanism, body, getineqconstraint(mechanism, cid))
    end

    return body.f
end

@inline function ∂dyn∂vel(body::Body{T}, Δt) where T
    J = body.J
    ωnew = getωnew(body)
    sq = sqrt(4 / Δt^2 - ωnew' * ωnew)

    dynT = SMatrix{3,3,T,9}(body.m / Δt * I)
    dynR = skewplusdiag(ωnew, sq) * J - J * ωnew * (ωnew' / sq) - skew(J * ωnew)

    Z = @SMatrix zeros(T, 3, 3)

    return [[dynT; Z] [Z; dynR]]
end

@inline torqueFromForce(F::AbstractVector{T}, r::AbstractVector{T}) where T = cross(r, F)
@inline function setForce!(body::Body, F, τ, No)
    for i = 1:No
        body.F[i] = F
        body.τ[i] = τ
    end
end

ext: jl


text:
using LightGraphs
using MetaGraphs

# Imports
include("entities/MetaGraphsDFG.jl")
include("services/MetaGraphsDFG.jl")
include("DFGPlots.jl")

# Exports
# Deprecated -
# export MetaGraphsDFG
export exists
export getLabelDict, getDescription, setDescription, getAddHistory, getSolverParams, setSolverParams
#
export getAddHistory, getDescription, getLabelDict
export addVariable!, addFactor!
export ls, lsf, getVariables, getFactors, getVariableIds, getFactorIds
export getVariable, getFactor
export updateVariable!, updateFactor!
export deleteVariable!, deleteFactor!
export getIncidenceMatrix, getIncidenceMatrixSparse
export getNeighbors
export getSubgraphAroundNode
export getSubgraph
export isFullyConnected, hasOrphans
export toDot, toDotFile


text:
macro R_str(s)
	s
end

y = b"""
"\xa8br\x8bjr\""
"nq"
"zjrfcpbktjmrzgsz\xcaqsc\x03n\"huqab"
"daz\\zyyxddpwk"
"draes\xa2n\\g\x27ek\"lj\"\\viqych"
"nnx\\krnrfomdnt\x2flbl\xd2xpo\"cp\"k"
"kwdaapalq"
"u\"ptk"
"ckhorczuiudfjmmcc\\u\"wozqxibsfjma"
"ydctdrxat\"pd\"lwi\"bjesevfw\xe8"
"v\"\xa8rrzep\"\"r"
"nbydghkfvmq\\\xe0\"lfsrsvlsj\"i\x61liif"
"jsas\"u"
"odipikxlo"
"\"rnubsgwltqkbsu\"pcpcs"
"eitk\\f\\mhcqqoym\\ji"
"vnedc"
"\"lhcaurdqzyjyu"
"haxzsa\"zcn\"y\"foclgtjfcnv\"m\x68krc"
"\"eoeggg\"tmiydvcay\"vfavc"
"snqvyqoncwxcvwbdktoywch"
"rnfgjsyr\xd5wacy"
"ik\"hebrpvsts"
"txw"
"\x15pxtdkogd\"urbm\"gevhh\"nxr\x3erxtk"
"cetqtcy"
"inleep\\mgl"
"uflwbxvww\x2cxzezqnaply\"yh\"qlllzk"
"eepak\"xqtedzt"
"na\x61qzfieafvyrsnwkssznohjmc"
"yceaonylz\xc1\\jrlbbkzwsidfi"
"ybqafngkcqpbp"
"\xaft"
"yidjpaobqydso"
"ju\\ldxig\\lrdrhjcmm\x77rc"
"tylacqeslnwj\x48ds\"tjxa"
"efbfm"
"\\fxkgoprgdcjgyajykg\\dtbrz"
"eujvva"
"h\x7acwfpikme\\vwthyvrqdnx\""
"rbpbrxm\\\"\"\"voxx"
"ykiw\"tkb\\lforu\"rsf\\tf\"x\"rqti"
"e\\wh\x77aqeugiq\\ihhfqfuaij"
"g\"t\\o"
"nxzo\"hf\\xp"
"dxiaqfo\xea"
"kali\\zczhiqkqzybjj\"fgdjnik"
"zdkgrqmdv"
"bimxim\xb6lrwsaj\"ui\"a"
"\"rrznitibgx\\olpsjmjqzctxaubdifsq"
"zb\"khzixaacmhuzmlymoformipdzml"
"qfwi"
"hjwsxfpphttjy\"\"zixais\xbblgnqfto"
"puj\\qmyu\"nqgaqfthbwjokbmrpbhpi"
"cyxdpkh\\\""
"q"
"m"
"tbxdzzllarlo"
"gbtys"
"gytilk\\vlqxvcuutjunrqc"
"uugkvcuzan\\eyhb"
"yaxr\"genlbgw\"\\uc"
"nrgecjeip\\sjdvgqaqxwsqactopu"
"pu\"r\"txpyrkfny\\zmwfneyvwmnkkdipv"
"jm\xa3bhwvq"
"qxojmnml\"w\x9airr"
"xbzsuihs\x4dcedy\xaclrhgii\\\""
"drgjirusrekrwmvxllwdm"
"\x28hfxnfpycmpnkku\"csuf\xaarxlqyg\"x"
"\"zvz\\rmg\"\\sxxoifffyqfyn\"iq\"ps"
"\"z"
"zbwkmk\"sgzos\x93gtc\""
"bvm\x28aa\\\\\"pywuhaniox\\z\\hbp\xd7mold"
"aszgvsyna"
"qf\"vdwuss"
"lnohni\"qwiacjsjegstlbfq\\kyjhyd"
"c\\naawulxlqplnacvytspry\xf5ytxxqq"
"razwqmsqgbaaxcd\\f"
"radggyrjrg\"zx"
"\"pu\x11t\\ajcjuieinlkvya"
"veggiskh"
"eglfhjxiet\"kouqfskwsy\"hpthsldel"
"mv\xc1b\"f\\shrssnjwcpmurepdxdlcj"
"dlayjd\"suvzotgdtc"
"\xa9pvxeopn"
"lpplsaxy\"oiwaq"
"hqwh\\lusv"
"hykykwlx\"\xa5atkgh\\d\x63dff"
"vfktanpjy\"xxetc"
"dnhwkgjnsmsswfuelvihvjl\"jtf"
"x\"dwvzra\"nbbsewftehczgbvfzd\"rau"
"csfi\"mzejnjqkqupwadrgti\"von"
"xckf\xf7xsm\\pgvlpetjndpyblais\\z"
"yecy\x6fuj\x58bwpgeuiw\"mdu"
"fgb"
"c\\lx\x3efthet\xfdelgvwvpem"
"kgyrmarvfwjinlowt"
"yzte"
"vc\"z"
"sxevqfzmmdwsuu\""
"fxbaercmcy\xb6md"
"f"
"m\x44gqbcppho\\b"
"gtafr\x57m\x11jy\"\"erwmmpiwjkbckuw"
"ufdjt\"kssprzxqixzxmq\x58q"
"yzbyo\"lfdbyaxexyfbnyv\\\xe8xmre"
"u\x43ntr\\\\byyfjr\"iveujvnwsqbnpuvrta"
"us\xf6bai"
"c\\edh"
"tzckolphexfq\\\x23\xfbdqv\\\"m"
"yjafhbvhhj\x1b\"bplb"
"\"o"
"rubahvmp\""
"qmkukrnrmqumh"
"wdpxyvyidhwjf\\nabbijwhr\xc5bksvy\"p"
"u\"prlpg\""
"nsvcquyxbwilsxxemf\xd9leq"
"y\xcetxuafl"
"it"
"kwdlysf\\xjpelae"
"viwh\x58wpjjlnvryuti\x2chngrx\\nhtkui"
"vhn\x9ehre\xc3ncsqbozms\"nl"
"ytc\xa3mgeeogjcqavmmmd"
"xzlexlitseozoxtpzzutfq"
"cish\x07lmovj"
"ekbflwqzaiivdr\"pq\\azrfbntqwkn"
"uc\"xdbegmlmhksofzohavtrnxf"
"xfdnrdqdrcjzbe"
"ndg\"ckgrpisib\"rg\"p\\lmpfzlssnvk"
"witfjwpbyyzlop"
"zonlww\"emrbcsgdtrg\"rjzy\x64zqntlw"
"dvgb\"zn\\vrbzema\"ckmd"
"\\vdlmxhlvldk\"pmzazeip"
"\"\"r"
"rsntinv"
"iy"
"lr\x20efh"
"csgexlb\"zqdavlxxhtdbh\"\"\x0fkpvhiphm"
"ouwhp\"ogbft"
"cm\\ckltng\"dw\x8brf\xf0eppgckd"
"zmnlsgalhpkejsizfsbtnfliu\"nhc"
"pnrkaayqvwpdjbhcrbb\"yfeq\"aq"
"ozh\\hoxow\x2csrtr\\r\""
"bqxabj\"u\\s"
"cpsjti\"gy"
"aa\"p\\nki\\ajijkqev"
"q\"\"lfdentjgd\\"
"bmokvpoebutfki"
"pielvcbne\xf6efvzxn"
"kx"
"zlgmqagcrbhrwtwtmmg"
"aiyhmntcqjbpv\xb5hhswxbryoedvos"
"tdpaxrb"
"fu\"\x7dttkyvhrlwko"
"oirc\\\"cqlnqffjqt\\k"
"edxlia\\tcyby"
"jpeybgwfayerfrfbvfog\"ol"
"ysr"
"bzwzilgwfugjk"
"tlcc\x75nukvwjgftetjcs\xaecwc"
"dsqssa\"vzrf\"sewbp\\ahhlmhbeihlh"
"qtgmjck\"n\"guki\"gmdivwqxismqj"
"\"f"
"wuorvlovucngbzdszqpikyk"
"dfrdsacoukmgvhbq\"\"iwto"
"\"ey\"ch\\wcgioe\\\"ouvligmsw"
"ciqlszzgs"
"\\tzyrkaoi\"sopjaq"
"lmtnv"
"ar\"fqoroigiertjjlm\"ymgi\\kkjewsxd"
"wehcimlvudpxtamdn\"rwy"
"hr\"zvrwthr\"vruzqfrldn\"b"
"sggekodkiwvym\"mhsco"
"ltlkfbrrdvk\\"
"uut\"sfjnz\"\\ef"
"hxilg\\"
"zsredsiwlzrpedibn"
"vtfi"
"\\h"
"qekfrc\xf6wduodbwrguqcng\\n"
"\"lljlfdrxftwidn\\pkv\xd9ij"
"mrvgqynpehkliuijlpp"
"gikjph"
"yoxcdrdt\"wbaurnyhoyxoihu"
"onmomwuxuammbzxe"
"rnrr\\twviz\x61gqaljr\x0dmtw"
"r\"vupaoi"
"l"
"sei"
"jwxtdtbkd\\kxd"
"\x22v\\"
"ahd"
"j\"bjqxs"
"\\i\x24gglxub\\nzsajokt"
"lviwpu\"uxdlh\\zuy\"xqy\"ytdzlx\"r"
"kptfmys"
"fwxzikfhczkjwyjszqdbkepaeellc"
"nlqpsvbrbd\\ns"
"qryuwkjiodw\"\"vaqyq\"dmyifm"
"tw\x15kdmaudjl\\zorhp\"alwh"
"aatrvczesykekkjfyb\"kb"
"usqcutbqbxxhucwxo\xc1ltb\"j\"bghjcvws"
"ilhsrnzxkz"
"bianqfdfdhvw"
"hqibqs\x7ax\"qoxqoaqtcsz"
"htxtoojbbauztwxuiq\\ngyfy\\obzc"
"rxn\\moxlj"
"mtus\x84erh\"dbe"
"asx\x50huvsitcxadt"
"\"bugggtnrc\"\"kl\"hmpu\x83hqrvhpo"
"ewisbp\"\"vuzf\\w\x5fvalszdhl"
"scusplpwxfnxu\x57\"zynpn\x99xerc\\ri"
"m\\kinmkke\x0cl"
"xhuzit\x7fd"
"kfbo\x04\x50ruqirn"
"t\"\"xpbdscmdoug"
"punvpsgnbgyxe\"sptmpz"
"bxukkazijr"
"nxyrcdaoo\"rjkk\"wntehcvcip\"vrd"
"rdpvqskmihqaw"
"p\\gwdhtqnpgthod"
"nwnuf\"\"yebycearom\"nqym\"\xd4sii\xccle"
"alda\"ptspo\"wkkv\"zoi\"hkb\"vnntyd"
"ixpgpfzbqv"
"znui\"\\fzn\x03qozabh\"rva\"pv\x67"
"e\"zswmwuk"
"hcccygwfa"
"ngmace\\rtyllolr\"\x68bw"
"\\c\"jyufbry\"ryo\"xpo\x26ecninfeckh\\s"
"hdnpngtuc\"dzbvvosn\x31fwtpzbrt"
"hesbpd\xd4"
"dsdbstuzrdfmrnyntufs\"dmv"
"d\xeeibcwhcvkt"
"fvzwrsfjdqdmy\"\"v"
"ns\"dqafz\\lkyoflnazv\"mn\x37\"o\"yj\"e"
"dypilgbwzccayxa\"bnmuernx"
"q\xa9ztqrhreb\"\"kxfeyodqb"
"iz\xa5qjxqulaawuwz\"rqmpcj\\yel"
"z\"\\pq\"\"y\x67zpjtn"
"ifxqvivp\"kiiftdoe"
"jxzebj\"\x35\"qr\"ecglcutuoyywqumcs\"kk"
"q"
"yob\x85qmpuwexptczbkrl"
"cjiavv\"uudpozvibyycnmxhxpxmpjoz"
"xro\\uiqyrcid"
"nod\\k"
"d\"neiec"
"tqyrqvwyvmz\\pzgzzcqsqsrgbqbtapoz"
"r\"xvocpeuxfxslgueb\x05kzyyie\"aoec"
"\"du\\uirlhcbgv\\cjqhfreqnvn"
"zp\x04\x15\"pbjwhrjtmiba"
"\\cv\""
"k\"rwnb\\hiu\"rqd\"rc\\nyakrhly"
"klrmafjzandiddodgz"
"xipzhqzhvlpykzcuppx"
"zdvrvn\xd0mtfvpylbn\\\\sxcznrzugwznl"
"ody\\pvm\"kpjiudzhxazirgxzvumeat\"o"
"kllvhdp\"prjikzrrc\"adgpegc\\\"gk"
"sqtpug\xbcaauxaamw"
"wegxxrrxdvpivrqievfeokmnojsk"
"\\bo"
"gijhz"
"ylowluvabwrigssdgtxdwsiorxev\xdd"
"\""
"ghnsrnsqtxpygikahkrl"
"\"rcfqkbjf\"sgxg\"vnd\\rotn"
"ap\"smgsuexjrbuqs\"mpbstogj\"x"
"koaunz\\sgt\"opv"
"yialiuzwix"
"yp\"ndxgwzml\"bt"
"lpcjxmggfsy\\szbxccarjkqzasqkb\xcfd\x0c"
"x"
"mgakc"
"vjieunoh\x73fjwx"
"erbvv\"qulsd"
"mimycrbfhqkarmz"
"tihfbgcszuej\"c\xfbvoqskkhbgpaddioo"
"mziavkwrmekriqghw"
"izk\\tnjd\\ed\\emokvjoc"
"c\"nhbqzndro\\g"
"usfngdo"
"aypljdftvptt"
"ym\"afvq\xbcc"
"zabi\"wjpvugwhl"
"ebvptcjqjhc\"n\"p\"dxrphegr\\"
"mzlqqxokhye\xd9\\rffhnzs"
"hnipqknwpsjakanuewe"
"rqgbfcjdrmiz\"h"
"kzzp\\z\\txmkwaouxictybwx"
"yzmspjkqrteiydswlvb"
"gjpxklgpzv\"txri\\hotpuiukzzzd"
"p\"rxergtbsxmjmkeeqwvoagnki\""
"santipvuiq"
"\"ihjqlhtwbuy\"hdkiv\"mtiqacnf\\"
"oliaggtqyyx"
"fwwnpmbb"
"yrtdrieazfxyyneo"
"nywbv\\"
"twc\\ehfqxhgomgrgwpxyzmnkioj"
"qludrkkvljljd\\xvdeum\x4e" """

x = R"""
"\xa8br\x8bjr\""
"nq"
"zjrfcpbktjmrzgsz\xcaqsc\x03n\"huqab"
"daz\\zyyxddpwk"
"draes\xa2n\\g\x27ek\"lj\"\\viqych"
"nnx\\krnrfomdnt\x2flbl\xd2xpo\"cp\"k"
"kwdaapalq"
"u\"ptk"
"ckhorczuiudfjmmcc\\u\"wozqxibsfjma"
"ydctdrxat\"pd\"lwi\"bjesevfw\xe8"
"v\"\xa8rrzep\"\"r"
"nbydghkfvmq\\\xe0\"lfsrsvlsj\"i\x61liif"
"jsas\"u"
"odipikxlo"
"\"rnubsgwltqkbsu\"pcpcs"
"eitk\\f\\mhcqqoym\\ji"
"vnedc"
"\"lhcaurdqzyjyu"
"haxzsa\"zcn\"y\"foclgtjfcnv\"m\x68krc"
"\"eoeggg\"tmiydvcay\"vfavc"
"snqvyqoncwxcvwbdktoywch"
"rnfgjsyr\xd5wacy"
"ik\"hebrpvsts"
"txw"
"\x15pxtdkogd\"urbm\"gevhh\"nxr\x3erxtk"
"cetqtcy"
"inleep\\mgl"
"uflwbxvww\x2cxzezqnaply\"yh\"qlllzk"
"eepak\"xqtedzt"
"na\x61qzfieafvyrsnwkssznohjmc"
"yceaonylz\xc1\\jrlbbkzwsidfi"
"ybqafngkcqpbp"
"\xaft"
"yidjpaobqydso"
"ju\\ldxig\\lrdrhjcmm\x77rc"
"tylacqeslnwj\x48ds\"tjxa"
"efbfm"
"\\fxkgoprgdcjgyajykg\\dtbrz"
"eujvva"
"h\x7acwfpikme\\vwthyvrqdnx\""
"rbpbrxm\\\"\"\"voxx"
"ykiw\"tkb\\lforu\"rsf\\tf\"x\"rqti"
"e\\wh\x77aqeugiq\\ihhfqfuaij"
"g\"t\\o"
"nxzo\"hf\\xp"
"dxiaqfo\xea"
"kali\\zczhiqkqzybjj\"fgdjnik"
"zdkgrqmdv"
"bimxim\xb6lrwsaj\"ui\"a"
"\"rrznitibgx\\olpsjmjqzctxaubdifsq"
"zb\"khzixaacmhuzmlymoformipdzml"
"qfwi"
"hjwsxfpphttjy\"\"zixais\xbblgnqfto"
"puj\\qmyu\"nqgaqfthbwjokbmrpbhpi"
"cyxdpkh\\\""
"q"
"m"
"tbxdzzllarlo"
"gbtys"
"gytilk\\vlqxvcuutjunrqc"
"uugkvcuzan\\eyhb"
"yaxr\"genlbgw\"\\uc"
"nrgecjeip\\sjdvgqaqxwsqactopu"
"pu\"r\"txpyrkfny\\zmwfneyvwmnkkdipv"
"jm\xa3bhwvq"
"qxojmnml\"w\x9airr"
"xbzsuihs\x4dcedy\xaclrhgii\\\""
"drgjirusrekrwmvxllwdm"
"\x28hfxnfpycmpnkku\"csuf\xaarxlqyg\"x"
"\"zvz\\rmg\"\\sxxoifffyqfyn\"iq\"ps"
"\"z"
"zbwkmk\"sgzos\x93gtc\""
"bvm\x28aa\\\\\"pywuhaniox\\z\\hbp\xd7mold"
"aszgvsyna"
"qf\"vdwuss"
"lnohni\"qwiacjsjegstlbfq\\kyjhyd"
"c\\naawulxlqplnacvytspry\xf5ytxxqq"
"razwqmsqgbaaxcd\\f"
"radggyrjrg\"zx"
"\"pu\x11t\\ajcjuieinlkvya"
"veggiskh"
"eglfhjxiet\"kouqfskwsy\"hpthsldel"
"mv\xc1b\"f\\shrssnjwcpmurepdxdlcj"
"dlayjd\"suvzotgdtc"
"\xa9pvxeopn"
"lpplsaxy\"oiwaq"
"hqwh\\lusv"
"hykykwlx\"\xa5atkgh\\d\x63dff"
"vfktanpjy\"xxetc"
"dnhwkgjnsmsswfuelvihvjl\"jtf"
"x\"dwvzra\"nbbsewftehczgbvfzd\"rau"
"csfi\"mzejnjqkqupwadrgti\"von"
"xckf\xf7xsm\\pgvlpetjndpyblais\\z"
"yecy\x6fuj\x58bwpgeuiw\"mdu"
"fgb"
"c\\lx\x3efthet\xfdelgvwvpem"
"kgyrmarvfwjinlowt"
"yzte"
"vc\"z"
"sxevqfzmmdwsuu\""
"fxbaercmcy\xb6md"
"f"
"m\x44gqbcppho\\b"
"gtafr\x57m\x11jy\"\"erwmmpiwjkbckuw"
"ufdjt\"kssprzxqixzxmq\x58q"
"yzbyo\"lfdbyaxexyfbnyv\\\xe8xmre"
"u\x43ntr\\\\byyfjr\"iveujvnwsqbnpuvrta"
"us\xf6bai"
"c\\edh"
"tzckolphexfq\\\x23\xfbdqv\\\"m"
"yjafhbvhhj\x1b\"bplb"
"\"o"
"rubahvmp\""
"qmkukrnrmqumh"
"wdpxyvyidhwjf\\nabbijwhr\xc5bksvy\"p"
"u\"prlpg\""
"nsvcquyxbwilsxxemf\xd9leq"
"y\xcetxuafl"
"it"
"kwdlysf\\xjpelae"
"viwh\x58wpjjlnvryuti\x2chngrx\\nhtkui"
"vhn\x9ehre\xc3ncsqbozms\"nl"
"ytc\xa3mgeeogjcqavmmmd"
"xzlexlitseozoxtpzzutfq"
"cish\x07lmovj"
"ekbflwqzaiivdr\"pq\\azrfbntqwkn"
"uc\"xdbegmlmhksofzohavtrnxf"
"xfdnrdqdrcjzbe"
"ndg\"ckgrpisib\"rg\"p\\lmpfzlssnvk"
"witfjwpbyyzlop"
"zonlww\"emrbcsgdtrg\"rjzy\x64zqntlw"
"dvgb\"zn\\vrbzema\"ckmd"
"\\vdlmxhlvldk\"pmzazeip"
"\"\"r"
"rsntinv"
"iy"
"lr\x20efh"
"csgexlb\"zqdavlxxhtdbh\"\"\x0fkpvhiphm"
"ouwhp\"ogbft"
"cm\\ckltng\"dw\x8brf\xf0eppgckd"
"zmnlsgalhpkejsizfsbtnfliu\"nhc"
"pnrkaayqvwpdjbhcrbb\"yfeq\"aq"
"ozh\\hoxow\x2csrtr\\r\""
"bqxabj\"u\\s"
"cpsjti\"gy"
"aa\"p\\nki\\ajijkqev"
"q\"\"lfdentjgd\\"
"bmokvpoebutfki"
"pielvcbne\xf6efvzxn"
"kx"
"zlgmqagcrbhrwtwtmmg"
"aiyhmntcqjbpv\xb5hhswxbryoedvos"
"tdpaxrb"
"fu\"\x7dttkyvhrlwko"
"oirc\\\"cqlnqffjqt\\k"
"edxlia\\tcyby"
"jpeybgwfayerfrfbvfog\"ol"
"ysr"
"bzwzilgwfugjk"
"tlcc\x75nukvwjgftetjcs\xaecwc"
"dsqssa\"vzrf\"sewbp\\ahhlmhbeihlh"
"qtgmjck\"n\"guki\"gmdivwqxismqj"
"\"f"
"wuorvlovucngbzdszqpikyk"
"dfrdsacoukmgvhbq\"\"iwto"
"\"ey\"ch\\wcgioe\\\"ouvligmsw"
"ciqlszzgs"
"\\tzyrkaoi\"sopjaq"
"lmtnv"
"ar\"fqoroigiertjjlm\"ymgi\\kkjewsxd"
"wehcimlvudpxtamdn\"rwy"
"hr\"zvrwthr\"vruzqfrldn\"b"
"sggekodkiwvym\"mhsco"
"ltlkfbrrdvk\\"
"uut\"sfjnz\"\\ef"
"hxilg\\"
"zsredsiwlzrpedibn"
"vtfi"
"\\h"
"qekfrc\xf6wduodbwrguqcng\\n"
"\"lljlfdrxftwidn\\pkv\xd9ij"
"mrvgqynpehkliuijlpp"
"gikjph"
"yoxcdrdt\"wbaurnyhoyxoihu"
"onmomwuxuammbzxe"
"rnrr\\twviz\x61gqaljr\x0dmtw"
"r\"vupaoi"
"l"
"sei"
"jwxtdtbkd\\kxd"
"\x22v\\"
"ahd"
"j\"bjqxs"
"\\i\x24gglxub\\nzsajokt"
"lviwpu\"uxdlh\\zuy\"xqy\"ytdzlx\"r"
"kptfmys"
"fwxzikfhczkjwyjszqdbkepaeellc"
"nlqpsvbrbd\\ns"
"qryuwkjiodw\"\"vaqyq\"dmyifm"
"tw\x15kdmaudjl\\zorhp\"alwh"
"aatrvczesykekkjfyb\"kb"
"usqcutbqbxxhucwxo\xc1ltb\"j\"bghjcvws"
"ilhsrnzxkz"
"bianqfdfdhvw"
"hqibqs\x7ax\"qoxqoaqtcsz"
"htxtoojbbauztwxuiq\\ngyfy\\obzc"
"rxn\\moxlj"
"mtus\x84erh\"dbe"
"asx\x50huvsitcxadt"
"\"bugggtnrc\"\"kl\"hmpu\x83hqrvhpo"
"ewisbp\"\"vuzf\\w\x5fvalszdhl"
"scusplpwxfnxu\x57\"zynpn\x99xerc\\ri"
"m\\kinmkke\x0cl"
"xhuzit\x7fd"
"kfbo\x04\x50ruqirn"
"t\"\"xpbdscmdoug"
"punvpsgnbgyxe\"sptmpz"
"bxukkazijr"
"nxyrcdaoo\"rjkk\"wntehcvcip\"vrd"
"rdpvqskmihqaw"
"p\\gwdhtqnpgthod"
"nwnuf\"\"yebycearom\"nqym\"\xd4sii\xccle"
"alda\"ptspo\"wkkv\"zoi\"hkb\"vnntyd"
"ixpgpfzbqv"
"znui\"\\fzn\x03qozabh\"rva\"pv\x67"
"e\"zswmwuk"
"hcccygwfa"
"ngmace\\rtyllolr\"\x68bw"
"\\c\"jyufbry\"ryo\"xpo\x26ecninfeckh\\s"
"hdnpngtuc\"dzbvvosn\x31fwtpzbrt"
"hesbpd\xd4"
"dsdbstuzrdfmrnyntufs\"dmv"
"d\xeeibcwhcvkt"
"fvzwrsfjdqdmy\"\"v"
"ns\"dqafz\\lkyoflnazv\"mn\x37\"o\"yj\"e"
"dypilgbwzccayxa\"bnmuernx"
"q\xa9ztqrhreb\"\"kxfeyodqb"
"iz\xa5qjxqulaawuwz\"rqmpcj\\yel"
"z\"\\pq\"\"y\x67zpjtn"
"ifxqvivp\"kiiftdoe"
"jxzebj\"\x35\"qr\"ecglcutuoyywqumcs\"kk"
"q"
"yob\x85qmpuwexptczbkrl"
"cjiavv\"uudpozvibyycnmxhxpxmpjoz"
"xro\\uiqyrcid"
"nod\\k"
"d\"neiec"
"tqyrqvwyvmz\\pzgzzcqsqsrgbqbtapoz"
"r\"xvocpeuxfxslgueb\x05kzyyie\"aoec"
"\"du\\uirlhcbgv\\cjqhfreqnvn"
"zp\x04\x15\"pbjwhrjtmiba"
"\\cv\""
"k\"rwnb\\hiu\"rqd\"rc\\nyakrhly"
"klrmafjzandiddodgz"
"xipzhqzhvlpykzcuppx"
"zdvrvn\xd0mtfvpylbn\\\\sxcznrzugwznl"
"ody\\pvm\"kpjiudzhxazirgxzvumeat\"o"
"kllvhdp\"prjikzrrc\"adgpegc\\\"gk"
"sqtpug\xbcaauxaamw"
"wegxxrrxdvpivrqievfeokmnojsk"
"\\bo"
"gijhz"
"ylowluvabwrigssdgtxdwsiorxev\xdd"
"\""
"ghnsrnsqtxpygikahkrl"
"\"rcfqkbjf\"sgxg\"vnd\\rotn"
"ap\"smgsuexjrbuqs\"mpbstogj\"x"
"koaunz\\sgt\"opv"
"yialiuzwix"
"yp\"ndxgwzml\"bt"
"lpcjxmggfsy\\szbxccarjkqzasqkb\xcfd\x0c"
"x"
"mgakc"
"vjieunoh\x73fjwx"
"erbvv\"qulsd"
"mimycrbfhqkarmz"
"tihfbgcszuej\"c\xfbvoqskkhbgpaddioo"
"mziavkwrmekriqghw"
"izk\\tnjd\\ed\\emokvjoc"
"c\"nhbqzndro\\g"
"usfngdo"
"aypljdftvptt"
"ym\"afvq\xbcc"
"zabi\"wjpvugwhl"
"ebvptcjqjhc\"n\"p\"dxrphegr\\"
"mzlqqxokhye\xd9\\rffhnzs"
"hnipqknwpsjakanuewe"
"rqgbfcjdrmiz\"h"
"kzzp\\z\\txmkwaouxictybwx"
"yzmspjkqrteiydswlvb"
"gjpxklgpzv\"txri\\hotpuiukzzzd"
"p\"rxergtbsxmjmkeeqwvoagnki\""
"santipvuiq"
"\"ihjqlhtwbuy\"hdkiv\"mtiqacnf\\"
"oliaggtqyyx"
"fwwnpmbb"
"yrtdrieazfxyyneo"
"nywbv\\"
"twc\\ehfqxhgomgrgwpxyzmnkioj"
"qludrkkvljljd\\xvdeum\x4e" """

inpt = R"""
""
"abc"
"aaa\"aaa"
"\x27" """

count = 0
for a in split(x, "\n")
	s = strip(a)[2:end-1]
	# print(s, " ")
	count += 2 # for the quotes
	count += 2*length(matchall(r"\\\\\\\\", s))
	s = replace(s, "\\\\\\\\", "z")
	count += length(matchall(r"\"", s))
	s = replace(s, "\"", "z")
	count += length(matchall(r"\\\\", s))
	s = replace(s, "\\\\", "z")
	count += 3 * length(matchall(r"\\x", s))
	s = replace(s, "\\x", "z")
	# println(s)
end

println(count)

count = 0
for a in split(x, "\n")
	st = strip(a)[2:end-1]
	print(st, " ")
	count += 4 # for the quotes
	for s in st
		if s == '\\'
			count += 1
			println("found \\")
		elseif s == '\"'
			count += 2
			println("found \"")
		end
	end
	println(st)
end

println(count)

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
module TestLyceumBase

include("preamble.jl")

@includetests "LyceumBase.jl"

end # module


text:
"""
    biconnectedComponents(network, ignoreTrivial=false)

Calculate biconnected components (aka "blobs") using Tarjan's algorithm:
the output is an array of arrays of edges.
These blobs are returned in post-order, but within a blob,
edges are *not* necessarily sorted in topological order.
If `ignoreTrivial` is true, trivial components (of a single edge)
are not returned.
The network is assumed to be connected.

**Warnings**: for nodes, fields `k`, `inCycle`, and `prev`
are modified during the algorithm. They are used to store the
node's "index" (time of visitation), "lowpoint", and the node's
"parent", as defined by the order in which nodes are visited.

References:
- p. 153 of [Tarjan (1972)](http://epubs.siam.org/doi/pdf/10.1137/0201010).
  Depth first search and linear graph algorithms,
  SIAM Journal on Computing, 1(2):146-160
- on [geeksforgeeks](https://www.geeksforgeeks.org/biconnected-components/),
  there is an error (as of 2018-01-30):
  `elif v != parent[u] and low[u] > disc[v]:` (python version)
  should be replaced by
  `elif v != parent[u] and disc[u] > disc[v]:`
- nice explanation at this
  [url](https://www.cs.cmu.edu/~avrim/451f12/lectures/biconnected.pdf)
"""
function biconnectedComponents(net, ignoreTrivial=false::Bool)
    for n in net.node
        n.inCycle = -1 # -1 for missing, node not visited yet
        n.k = -1       # inCycle = lowpoint
        n.prev = nothing # parent: will be a node later
    end
    S = Edge[] # temporary stack
    blobs = Vector{Edge}[] # will contain the blobs
    biconnectedComponents(net.node[net.root], [0], S, blobs, ignoreTrivial)
    # if stack not empty: create last connected component
    if length(S)>0
        #println("stack not empty at the end.")
        bb = Edge[]
        while length(S)>0
            e2 = pop!(S)
            push!(bb, e2)
        end
        if !ignoreTrivial || length(bb)>1
            push!(blobs, bb)
        end
    end
    return blobs
end

"""
    biconnectedComponents(node, index, S, blobs, ignoreTrivial)

Helper recursive function starting at a node (not a network).
`index` is an array containing a single integer, thus mutable:
order in which nodes are visited.
"""
function biconnectedComponents(node, index, S, blobs, ignoreTrivial)
    #println("\nentering biconnect, index=$(index[1])")
    children = 0
    # set depth index for v to the smallest unused index
    node.k = index[1]       # index / disc
    node.inCycle = index[1] # lowpoint / lowlink
    #println(" for node $(node.number): k=$(node.k), low=$(node.inCycle)")
    index[1] += 1
    #print(" stack: "); @show [e.number for e in S]

    for e in node.edge # each (v, w) in E do
        w = (e.node[1] == node) ? e.node[2] : e.node[1]
        #println(" w: $(w.number) along edge $(e.number)")
        if w.k == -1 # w not yet visited
            w.prev = node # set parent of w to current node
            #println(" w's parent = $(w.prev.number)")
            children += 1
            push!(S, e)
            #println(" edge $(e.number) on stack, intermediate step for node=$(node.number)")
            biconnectedComponents(w, index, S, blobs, ignoreTrivial)
            # check if subtree rooted at w has connection to
            # one of the ancestors of node
            # Case 1 -- per Strongly Connected Components Article
            node.inCycle = min(node.inCycle, w.inCycle) # lowpoint
            #print(" case 1, low=$(node.inCycle) for node $(node.number); ")
            #println("low=$(w.inCycle) for w $(w.number)")
            # if node is an articulation point: pop until e
            if (node.prev == nothing && children > 1) ||
               (node.prev != nothing && w.inCycle >= node.k)
                # node is either root or an articulation point
                # start a new strongly connected component
                bb = Edge[]
                while length(S)>0
                    e2 = pop!(S)
                    #println(" popping: edge $(e2.number)")
                    push!(bb, e2)
                    e2 != e || break
                end
                if !ignoreTrivial || length(bb)>1
                    push!(blobs, bb)
                end
            end
        elseif w != node.prev && node.k > w.k
            # e is a back edge, not cross edge. Case 2 in article.
            node.inCycle = min(node.inCycle, w.k)
            #println(" case 2, node $(node.number): low=$(node.inCycle)")
            push!(S, e)
            #println(" edge $(e.number) now on stack, final step for node $(node.number)")
        else
            #println(" nothing to do from node $(node.number) to w $(w.number) along edge $(e.number)")
        end
    end
    return nothing
end

"""
    blobInfo(network, ignoreTrivial=true)

Calculate the biconnected components (blobs) using function
[`biconnectedComponents`](@ref) then:
- set node field `isExtBadTriangle` to true at the root of each
  non-trivial blob (and at the network root), false otherwise.
  (a better name for the field would be something like "isBlobRoot".)
- output:
  1. array of nodes that are the roots of each non-trivial blob,
     and the network root. If the root of the full network is
     not part of a non-trivial blob, a corresponding blob is
     added to the list.
  2. array of arrays: for each non-trivial blob,
     array of major hybrid edges in that blob.
  3. array of arrays: same as #2 but for minor hybrid edges,
     with hybrids listed in the same order, for each blob.

Blobs are ordered in reverse topological ordering
(aka post order).
If `ignoreTrivial` is true, trivial components are ignored.

keyword argument: `checkPreorder`, true by default. If false,
the `isChild1` edge field and the `net.nodes_changed` network field
are supposed to be correct.
"""
function blobInfo(net, ignoreTrivial=true::Bool;
    checkPreorder=true::Bool)
    if checkPreorder
      directEdges!(net) # update isChild1, needed for preorder
      preorder!(net) # creates / updates net.nodes_changed
    end
    bcc = biconnectedComponents(net, ignoreTrivial)
    bccRoots = Node[]         # one root node for each blob
    bccMajor = Vector{Edge}[] # one array for each blob
    bccMinor = Vector{Edge}[]
    for bicomp in bcc
        bccMa = Edge[]
        jmin = length(net.node)
        for edge in bicomp
            n = getParent(edge)
            j = something(findfirst(x -> x===n, net.nodes_changed), 0)
            jmin = min(j, jmin)
            if edge.hybrid && edge.isMajor
                push!(bccMa, edge)
            end
        end
        push!(bccRoots, net.nodes_changed[jmin])
        push!(bccMajor, bccMa)
        bccmi = Edge[] # find minor hybrid edges, in same order
        for edge in bccMa
            e = getPartner(edge)
            !e.isMajor || @warn "major edge $(edge.number) has a major partner: edge $(e.number)"
            push!(bccmi, e)
        end
        push!(bccMinor, bccmi)
    end
    # add the network root, if it was in a trivial bi-component (no hybrids)
    # blobs in post-order, so if there the root is there, it's the last blob
    if length(bcc)==0 || bccRoots[end]!=net.node[net.root]
        push!(bccRoots,net.node[net.root])
        push!(bccMajor, Edge[])
        push!(bccMinor, Edge[])
    end
    # update `isExtBadTriangle` to mark nodes that are blob roots:
    # these nodes will serve as dummy leaves when reached from other blobs
    for n in net.node n.isExtBadTriangle = false; end
    for r in bccRoots r.isExtBadTriangle = true;  end
    return(bccRoots, bccMajor, bccMinor)
end

"""
    blobDecomposition!(network)
    blobDecomposition(network)

Find blobs using [`biconnectedComponents`](@ref); find their roots
using [`blobInfo`](@ref); create a forest in the form of a
disconnected network (for efficiency), by deconnecting the
root of each non-trivial blob from its parent.
The root of each blob corresponds to a new leaf
(in another tree of the forest):
the number of the blob's root is given to the newly created leaf.

The first (bang) version modifies the network and returns
the array of blob roots. The second version copies the network
then returns a tuple: the forest and the blob roots.
"""
function blobDecomposition(net)
    net2 = deepcopy(net)
    blobR = blobDecomposition!(net2)
    return net2, blobR
end
function blobDecomposition!(net)
    nextnumber = maximum([n.number for n in net.node])+1
    blobR, tmp, tmp = blobInfo(net, true) # true: ignore trivial single-edge blobs
    for r in blobR
        for e in r.edge
            r == e.node[e.isChild1 ? 1 : 2] || continue
            removeEdge!(r,e) # detach edge e from root r
            removeNode!(r,e) # detach root r from edge e
            dummyleaf = Node(nextnumber, true) # true: leaf
            nextnumber += 1
            dummyleaf.name = string("dummy ", r.number)
            setEdge!(dummyleaf, e) # attach e to new dummy leaf
            setNode!(e,dummyleaf)  # attach new leaf to edge e
            e.isChild1 = false
            pushNode!(net, dummyleaf)
            break
        end
    end
    return blobR
end


text:
using GNSSTools
using ProgressMeter
using JLD
using PyPlot
pygui(true)
include("kalman_filter_gain_error_functions.jl")

function kf_run_results(sigtype, f_s, t_length;
                        f_d=0, fd_rate=0, phi_init=0, f_if=0,
                        include_phase_noise=true, include_thermal_noise=true,
                        include_adc=true, include_carrier=true,
                        acquisition_T=1e-3, M=1, fine_acq_T=10e-3,
                        tracking_T=1e-3, fine_acq_method=:fft,
                        CN0=45, code_start_idx=1, fd_range=5000, prn=1,
                        σᵩ²=GNSSTools.phase_noise_variance(CN0, tracking_T),
                        h_parms=[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2,
                        fd_center=0, nADC=8, dynamickf=true, state_num=3,
                        Tsys=535, channel="I", q_a=1, use_fine_acq=true,
                        dll_b=0.1, p_d=0.9, p_fa=1e-7, σω=1,
                        show_process_result_plot=false, fd_rate_init=0,
                        figsize=(8, 7.5), window_size=10,
                        save_to_file=missing, top=0.9, bottom=0.1,
                        left=0.1, right=0.9, wspace=0.15, hspace=0.4)
    sig_freq = sigtype.sig_freq
    signal = definesignal(sigtype, f_s, t_length; f_d=f_d,
                          prn=prn, fd_rate=fd_rate, f_if=f_if,
                          include_phase_noise=include_phase_noise,
                          include_thermal_noise=include_thermal_noise,
                          include_adc=include_adc, receiver_h_parms=h_parms,
                          include_carrier=include_carrier, nADC=nADC,
                          CN0=CN0, phi=phi_init, Tsys=Tsys,
                          code_start_idx=code_start_idx)
    generatesignal!(signal)
    results, trackresults, corr_result, SNR_est,
    P_d, above_threshold = process(signal, sigtype, 1, channel;
                                   return_corrresult=true,
                                   return_Pd=true,
                                   show_plot=show_process_result_plot,
                                   acquisition_T=acquisition_T,
                                   M=M,
                                   fine_acq_method=fine_acq_method,
                                   use_fine_acq=use_fine_acq,
                                   fine_acq_T=fine_acq_T,
                                   tracking_T=tracking_T,
                                   h₋₂=h₋₂,
                                   h₀=h₀,
                                   σᵩ²=σᵩ²,
                                   q_a=q_a,
                                   state_num=state_num,
                                   dynamickf=dynamickf,
                                   fd_rate=fd_rate_init,
                                   dll_b=dll_b,
                                   p_d=p_d,
                                   p_fa=p_fa,
                                   σω=σω)

    # Get the matrices used in signal tracking
    A = trackresults.A
    C = trackresults.C
    Q = trackresults.Q
    R = trackresults.R
    # Get KF gain values
    t = trackresults.t
    ZP = trackresults.ZP
    K = trackresults.K
    k₁ = K[1,:]
    k₂ = K[2,:]
    if state_num == 3
        k₃ = K[3,:]
    end
    B_exp_numerical = zeros(length(t))
    B_exp_analytical = zeros(length(t))
    t_window = t[window_size:window_size:end]
    dphi_meas = trackresults.dphi_meas
    dphi_filt = trackresults.dphi_filt
    σᵩ_unfilt = zeros(length(t_window))
    σᵩ_filt = zeros(length(t_window))
    pᵩ = zeros(length(t))
    for i in 1:length(t)
        # Calculate experimental steady state gain using the numerical solution
        B_exp_numerical[i] = get_B(init_H(K[:,i]...))
        # Calculate experimental steady state gain using the analytical equation
        B_exp_analytical[i] = calc_B(K[:,i]...)
        # Calculate experimental pᵩ
        P = trackresults.P_full[i,:,:]
        pᵩ[i] = sqrt((C*P*C')[1])
        if i%window_size == 0
            # Calculate standard deviation of unfiltered carrier phase error
            σᵩ_unfilt[Int(i/window_size)] = std(dphi_meas[i-(window_size-1):i])
            # Calculate standard deviation of filtered carrier phase error
            σᵩ_filt[Int(i/window_size)] = std(dphi_filt[i-(window_size-1):i])
        end
    end

    # Calculate the theoretical steady state kalman gain
    K_steady_state = dkalman(A, C, Q, GNSSTools.Diagonal(R))
    k₁_steady_state = K_steady_state[1,:]
    k₂_steady_state = K_steady_state[2,:]
    if state_num == 3
        k₃_steady_state = K_steady_state[3,:]
    end
    # Calculate theoretical steady state gain using the numerical solution
    B_theo_numerical = get_B(init_H(K_steady_state...))
    # Calculate theoretical steady state gain using the analytical equation
    B_theo_analytical = calc_B(K_steady_state...)
    # Calculate theoretical pᵩ
    P_steady_state = GNSSTools.dare(A', C', Q, GNSSTools.Diagonal(R)')
    pᵩ_steady_state = sqrt((C*P_steady_state*C')[1])
    # Calculate the steady state unfiltered carrier phase error (σᵩ)
    σᵩ_unfilt_steady_state = sqrt(R[1])
    σᵩ_filt_steady_state_num = NaN
    σᵩ_filt_steady_state_ana = NaN
    if include_phase_noise  # Calculate σᵩ with oscillator phase noise
        # Calculate the steady state filtered carrier phase error using
        # numerical solution
        σᵩ_filt_steady_state_num = sqrt(phase_noise_variance(B_theo_numerical,
                                                            CN0, tracking_T,
                                                            h₀, h₋₁, h₋₂))
        # Calculate the steady state filtered carrier phase error using
        # analytical equation
        try
            σᵩ_filt_steady_state_ana = sqrt(phase_noise_variance(B_theo_analytical,
                                                                CN0, tracking_T,
                                                                h₀, h₋₁, h₋₂))
        catch
            σᵩ_filt_steady_state_ana = NaN
            @warn("Unable to calculate value for `σᵩ_filt_steady_state_ana` due to negative sqrt.")
        end
    else # Calculate σᵩ without oscillator phase noise
        # Calculate the steady state filtered carrier phase error using
        # numerical solution
        σᵩ_filt_steady_state_num = sqrt(phase_noise_variance(B_theo_numerical,
                                                            CN0, tracking_T))
        # Calculate the steady state filtered carrier phase error using
        # analytical equation
        try
            σᵩ_filt_steady_state_ana = sqrt(phase_noise_variance(B_theo_analytical,
                                                                CN0, tracking_T))
        catch
            σᵩ_filt_steady_state_ana = NaN
            @warn("Unable to calculate value for `σᵩ_filt_steady_state_ana` due to negative sqrt.")
        end
    end

    # Make plot
    fig = figure(figsize=figsize)
    rows = 3
    cols = 2

    # KF gains
    ax1 = fig.add_subplot(rows, cols, 1)
    ax1.plot(t, k₁, "-", color="blue", label="k₁")
    ax1.plot(t, k₂, "-", color="orange", label="k₂")
    ax1.plot([t[1], t[end]], fill(k₁_steady_state, 2), ":", color="blue")
    ax1.plot([t[1], t[end]], fill(k₂_steady_state, 2), ":", color="orange")
    if state_num == 3
        ax1.plot(t, k₃, "-", color="green", label="k₃")
        ax1.plot([t[1], t[end]], fill(k₃_steady_state, 2), ":", color="green")
    end
    ax1.plot([], [], "k:", label="Steady State kₙ")
    xlabel("Time (s)")
    ylabel("KF Gain")
    legend()
    title("(a)")

    # Filtered carrier phase errors based off P
    ax2 = fig.add_subplot(rows, cols, 2)
    ax2.plot(t, pᵩ.*(180/π), "k-", label="pᵩ")
    ax2.plot([t[1], t[end]], fill(pᵩ_steady_state, 2).*(180/π), "k:",
             label="Steady State pᵩ")
    xlabel("Time (s)")
    ylabel("pᵩ (degrees)")
    legend()
    title("(b)")

    # Unfiltered carrier phase error
    ax3 = fig.add_subplot(rows, cols, 3)
    ax3.plot(t_window, σᵩ_unfilt.*(180/π), "k-", label="σᵩ Meas.")
    ax3.plot([t[1], t[end]], fill(σᵩ_unfilt_steady_state, 2).*(180/π), "k:",
              label="σᵩ Theo.")
    xlabel("Time (s)")
    ylabel("Unfiltered σᵩ (degrees)")
    legend()
    title("(c)")

    # Filtered carrier phase error
    ax4 = fig.add_subplot(rows, cols, 4)
    ax4.plot(t_window, σᵩ_filt.*(180/π), "k-", label="σᵩ Meas.")
    ax4.plot([t[1], t[end]], fill(σᵩ_filt_steady_state_num, 2).*(180/π), "k:",
              label="σᵩ Theo. Num.")
    # ax4.plot([t[1], t[end]], fill(σᵩ_filt_steady_state_ana, 2).*(180/π), ":",
            #   label="σᵩ Theo. Ana.")
    xlabel("Time (s)")
    ylabel("Filtered σᵩ (degrees)")
    legend()
    title("(d)")

    # Estimated equivalent KF noise bandwidth
    ax5 = fig.add_subplot(rows, cols, 5)
    ax5.plot(t, B_exp_numerical, "k-", label="B Exp. Num.")
    # ax5.plot(t, B_exp_analytical, ":", color="blue", label="B Exp. Ana.")
    ax5.plot([t[1], t[end]], fill(B_theo_numerical, 2), "k:",
             label="B Theo. Num.")
    # ax5.plot([t[1], t[end]], fill(B_theo_analytical, 2), ":", color="orange",
            #  label="B Theo. Ana.")
    xlabel("Time (s)")
    ylabel("Bandwidth (Hz)")
    ylim([-2, 15])
    legend()
    title("(e)")

    # Prompt correlator outputs
    ax6 = fig.add_subplot(rows, cols, 6)
    ax6.plot(t, real.(ZP), label="real(Zₚ)")
    ax6.plot(t, imag.(ZP), label="imag(Zₚ)")
    xlabel("Time (s)")
    ylabel("Zₚ")
    legend()
    title("(f)")

    subplots_adjust(hspace=hspace, wspace=wspace,  top=top, left=left, right=right,
                    bottom=bottom)
    if ~ismissing(save_to_file)
        savefig(save_to_file, dpi=300)
    end
end

tcxo_h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
ocxo_h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO
h_parms = tcxo_h_parms

t_length = 1
f_s = 2.5e6
f_if = 0
f_d = 1000
fd_rate = 0
phi_init = 0
code_start_idx = 1
CN0 = 45
nADC = 8

sigtype = define_l1ca_code_type()
channel = "I"

dll_b = 1
acquisition_T = 40e-3
M = 5
fine_acq_T = 100e-3
fine_acq_method = :fft
tracking_T = 1e-3
q_a = 0.1
σω = 1
state_num = 2
window_size = 50
fd_rate_init = 0

include_thermal_noise = true
include_phase_noise = true
include_adc = true
include_carrier = true

show_process_result_plot = false

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/figures/"
save_to_file = string(directory, "ch4_kf_response_successful_track_$(state_num)-state.pdf")
# save_to_file = string(directory, "ch4_kf_response_unsuccessful_track_$(state_num)-state.pdf")

kf_run_results(sigtype, f_s, t_length;
                        f_d=f_d, fd_rate=fd_rate, phi_init=phi_init, f_if=f_if,
                        include_phase_noise=include_phase_noise,
                        include_thermal_noise=include_thermal_noise,
                        include_adc=include_adc, include_carrier=include_carrier,
                        acquisition_T=acquisition_T, M=M,
                        fine_acq_T=fine_acq_T,
                        tracking_T=tracking_T,
                        fine_acq_method=fine_acq_method,
                        CN0=CN0, code_start_idx=code_start_idx,
                        h_parms=h_parms, nADC=nADC,
                        state_num=state_num,
                        channel=channel, q_a=q_a,
                        dll_b=dll_b,
                        show_process_result_plot=show_process_result_plot,
                        figsize=(8, 7.5), window_size=window_size,
                        top=0.95, bottom=0.1, left=0.1, right=0.9,
                        wspace=0.25, hspace=0.52, fd_rate_init=fd_rate_init,
                        save_to_file=save_to_file)