text:
module LifeValReserve

using LifeValStructures

export Calc, Prem

# probability of surviving t years from age x given yearly survival probabilities
function tapxFunction(apx)
    tapx = [1.0]
    n = length(apx)
    for i in 2:n
        append!(tapx, tapx[i-1] .* apx[i])
    end
    return(tapx)
end

function escFunction(initVal, incr, rate)
    esc = [initVal]
    n = length(incr)
    for i in 2:n
        append!(esc, esc[i-1] .* (1 .+ rate .* incr[i]))
    end
    return(esc)
end

# reusable function to set expenses per policy when pricing, because P is unknown.
# common code for any product tyee
function perpolExp(policy, runtype, basisData)
    eBasis = basisData[7][runtype.runno, policy.spc].datfram

    exp_init = eBasis.INIT_EXP_PREM * (1 + runtype.exp_inc)

    #renewal expense also for increasing premiums need to recalc with prem=1
    j = 1:policy.projMax
    m = (runtype.valmonth .+ j) .% 12
    incr = map(m -> ifelse.(m == policy.esc_m, 1, 0), m)
    exp_ren_prem = eBasis.REN_EXP_PREM .* escFunction(1, incr, policy.prem_esc)

    return [exp_init, exp_ren_prem]
end

# reusable function to set basis for calculating reserves
# common code for any product tyee
function setBasis(policy, runtype, basisData, profit)

    # set mort basis
    mort = ifelse(runtype.morta==1, ifelse(policy.gender==0, :UXM_D, :UXF_D),
                            ifelse(policy.gender==0, :UXM_DS, :UXF_DS))

    mortBasis = ifelse(profit==0, basisData[1][runtype.runno, policy.spc].datfram[:, mort], #qxtables
                                   basisData[2][runtype.runno, policy.spc].datfram[:, mort]) #qxptables
    #qxData[runno].datfram[:, [:AGE, :mort]]
    # set wdl basis
    wdl = ifelse(runtype.wdlsa == 1, :uxw, :uxws)
    wdlBasis = ifelse(profit==0, basisData[3][runtype.runno, policy.spc].datfram[:, wdl],
                                basisData[4][runtype.runno, policy.spc].datfram[:, wdl])

    # set v basis
    v = ifelse(runtype.valR == 1, :vt, :vts)
    vBasis = ifelse(profit==0, basisData[5][runtype.runno, policy.spc].datfram[:, v],
                            basisData[6][runtype.runno, policy.spc].datfram[:, v])

    j = 1:policy.projMax
    age = policy.age .+ (j .- 1) /12
    term_wx = policy.termIF .+ j

    #read select mortality
    #for k in termwx:(termwx+policy.projMax)
    #    l += l
    #    if ceil.(Int, k/12) = 1

    #    elseif ceil.(Int, k/12) = 2

    #    else
    #        uxd[k] = mortBasis[floor.(Int, age)]
    #    end
    #end

    uxd = mortBasis[floor.(Int, age)]
    uxw = wdlBasis[ceil.(Int, term_wx/12)]
    aqx = 1 .- exp.((-uxd .-uxw)/12)
    aqxd = uxd ./(uxd + uxw) .*aqx
    awxd = uxw ./(uxd + uxw) .*aqx
    pop!(aqx)
    apx = append!([1.0], 1 .- aqx)
    tapx = tapxFunction(apx)
    vt = vBasis[j]

    #premium and SumAssured escalations
    m = (runtype.valmonth .+ j) .% 12
    incr = map(m -> ifelse.(m == policy.esc_m, 1, 0), m)
    prem = escFunction(policy.premium, incr, policy.prem_esc)
    sa =  escFunction(policy.SumAssured, incr, policy.sa_esc)

    #expenses
    eBasis = ifelse(profit==0, basisData[7][runtype.runno, policy.spc].datfram, basisData[8][runtype.runno, policy.spc].datfram)

    exp_init = (eBasis.INIT_EXP_AMT + eBasis.INIT_EXP_SA * policy.SumAssured +
                eBasis.INIT_EXP_PREM * policy.premium) *
                (1 + runtype.exp_inc)

    exp_ren = eBasis.REN_EXP_AMT * (1 + runtype.exp_inc)
    exp_ren = exp_ren[1]

    exp_ren_prem = eBasis.REN_EXP_PREM .* prem
    exp_ren_prem = exp_ren_prem[1]

    exp_claim = (eBasis.CLM_EXP_AMT .+ eBasis.CLM_EXP_SA .* sa) .*
                 (1 + runtype.exp_inc)

    #sur_claim = (eBasis[runtype.runno, policy.spc].datfram[:SUR_EXP_AMT]  .+
    #             eBasis[runtype.runno, policy.spc].datfram[:SUR_EXP_SA] .* sa) .*
    #             (1 + runtype.exp_inc)

    #mat_claim = (eBasis[runtype.runno, policy.spc].datfram[:MAT_EXP_AMT]  .+
    #              eBasis[runtype.runno, policy.spc].datfram[:MAT_EXP_SA] .* sa) .*
    #              (1 + runtype.exp_inc)

    exp_infl = (1 .+ eBasis.REN_EXP_INFL) * (1 + runtype.expinfl_inc) .- 1
    exp_infl = exp_infl[1]

    expenses = exp_ren_prem .+ (exp_ren/12) * (1 + exp_infl) .^((j .- 1)/12)

    return [aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa]
end

function setAlloc(policy, runno, prem, fmc_sens, basisData, profit)
    global chData

    alloc_P = basisData[9][runno, policy.spc].datfram.ALOC_P
    BOS = basisData[9][runno, policy.spc].datfram.BOS
    fmc = basisData[9][runno, policy.spc].datfram.FMC_PERC

    # set ch basis
    j = 1:policy.projMax
    term_group = policy.termIF .+ j
    rate_alloc = alloc_P[ceil.(Int, term_group/12)]
    rate_BOS = BOS[ceil.(Int, term_group/12)]
    rate_fmc = fmc[ceil.(Int, term_group/12)] * fmc_sens

    prem_alloc = prem .* rate_alloc .* rate_BOS

    return [prem_alloc, rate_fmc]
end

function calcUF(policy, prem_alloc, fmc, unitPrice)
    fundbef = Array{Float64}(undef, policy.projMax)
    fundbeg = Array{Float64}(undef, policy.projMax)
    fundmc = Array{Float64}(undef, policy.projMax)
    units = Array{Float64}(undef, policy.projMax)
    fv = Array{Float64}(undef, policy.projMax)

    units_purc = prem_alloc ./ unitPrice[1:policy.projMax]
    units_start = policy.units_1 #1 fund
    units[1] = (units_start + units_purc[1])
    fundbeg[1] = units[1] .* unitPrice[1]
    fundbef[1] = units[1] .* unitPrice[2]
    fundmc[1] = fundbef[1] * fmc[1]
    fv[1] = fundbef[1] - fundmc[1]

    for j in 2:policy.projMax
        units[j] = fv[j - 1] / unitPrice[j] + units_purc[j]
        fundbeg[j] = units[j] * unitPrice[j]
        fundbef[j] = units[j] * unitPrice[j + 1]
        fundmc[j] = fundbef[j] * fmc[j]
        fv[j] = fundbef[j] - fundmc[j]
    end

    return [units, fv, fundmc, fundbef, fundbeg]#, fundch]
end

function calcNUF(policy, netPrem, fundch, initExp, renExp, aqxd, apx, unitf, unitPrice)
    units_initexp = Array{Float64}(undef, policy.projMax)
    units_renexp = Array{Float64}(undef, policy.projMax)
    units = Array{Float64}(undef, policy.projMax)
    non_u_int = Array{Float64}(undef, policy.projMax)
    mort_rate = Array{Float64}(undef, policy.projMax)
    death_cost = Array{Float64}(undef, policy.projMax)
    cf = Array{Float64}(undef, policy.projMax)
    nonexp = Array{Float64}(undef, policy.projMax)
    reserve = Array{Float64}(undef, policy.projMax)
    fundmc = Array{Float64}(undef, policy.projMax)

    units_purc = netPrem ./ unitPrice[1:policy.projMax]
    units_start = policy.units_99 #NU FUND
    units_initexp = initExp ./ unitPrice[1]
    units_renexp = renExp ./ unitPrice[1:policy.projMax]

    units[1] = units_start + units_purc[1] - units_initexp[1]
    nonexp[1] = units_initexp[1] * unitPrice[1]
    non_u_int[1] = units[1] .* (unitPrice[2] .- unitPrice[1])
    fundmc[1] = fundch[1]
    mort_rate[1] = aqxd[1]
    death_cost[1] = max((policy.SumAssured - unitf[1]),0) * mort_rate[1]
    cf[1] = units[1] * unitPrice[1] + non_u_int[1] + fundmc[1] - death_cost[1]

    for j in 2:policy.projMax
        units[j] = units_purc[j] - units_renexp[j] # + fv[j - 1] / unitPrice[j] ?
        nonexp[j] = units_renexp[j] * unitPrice[j]
        non_u_int[j] = units[j] .* (unitPrice[j+1] .- unitPrice[j])
        fundmc[j] = fundch[j]
        mort_rate[j] = aqxd[j]
        death_cost[j] = max((policy.SumAssured - unitf[j]),0) * mort_rate[j]
        cf[j] = units[j] * unitPrice[j] + non_u_int[j] + fundmc[j] - death_cost[j]
    end

    #calc reserves
    if cf[policy.projMax]< 0
        reserve[policy.projMax] = -cf[policy.projMax] / unitPrice[policy.projMax] * unitPrice[policy.projMax - 1]
    end

    for j in (policy.projMax-1):-1:1
        if cf[j] < 0
            reserve[j]= (-cf[j] + apx[j]*reserve[j + 1]) / unitPrice[j + 1] * unitPrice[j]
        else
            reserve[j]= 0
        end
    end

    return [units, cf, non_u_int, death_cost, nonexp, reserve]#, fundch]
end

function Calc(policy::LifeValStructures.LifeAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy,runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    cf = deathClaims + expenses .- prem
    reserve = sum(cf .* tapx .* vt)
    return(reserve)
end

function Calc(policy::LifeValStructures.Annuity, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa/12 .+ exp_claim) .* apx
    cf = deathClaims + expenses .- prem
    reserve = sum(cf .* tapx .* vt)
    return(reserve)
end

function Calc(policy::LifeValStructures.TermAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    cf = deathClaims + expenses .- prem
    reserve = sum(cf .* tapx .* vt)
    return(reserve)
end

function Calc(policy::LifeValStructures.EndowmentAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    deathClaims[policy.projMax]=deathClaims[policy.projMax]/aqxd[policy.projMax]
    cf = deathClaims + expenses .- prem
    reserve = sum(cf .* tapx .* vt)
    return(reserve)
end

function Calc(policy::LifeValStructures.EndowmentAssuranceUL, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    prem_alloc, fmc = setAlloc(policy, runtype.runno, prem, runtype.fmc_sens, basisData, 0)

    #currently only 1 fund, 100% allocation
    fmcp = (1 .+ fmc) .^ (1/12) .- 1 #fmc is annual, need to make it monthly

    #unit fund
    units, fundv, fundch, fundb, fundbeg = calcUF(policy, prem_alloc, fmcp, unitPrices[runtype.runno,1,1,:])
    #output=table((time = 1:policy.projMax, alloc = prem_alloc, fundvalbeg = fundbeg, fundvalbef = fundb, fmc = fundch, fundvalend = fundv); pkey = [:time])
    #writedlm("bu2.csv", output, ',')

    #non-unit fund
    non_units, noncf, nint, death_cost, nuexp, non_reserve = calcNUF(policy, (prem - prem_alloc), fundch, exp_init, expenses, aqxd, apx, fundv, unitPrices[runtype.runno,2,1,:])
    #output=table((time = 1:policy.projMax, alloc = (prem - prem_alloc), expenses = nuexp, interest = nint, fundmc = fundch, d_cost = death_cost, cf = noncf, v = non_reserve); pkey = [:time])
    #writedlm("nbu1.csv", output, ',')

    #profit test BASIS
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 1)
    prem_alloc, fmc = setAlloc(policy, runtype.runno, prem, runtype.fmc_sens, basisData, 1)
    #currently only 1 fund, 100% allocation
    fmcp = (1 .+ fmc) .^ (1/12) .- 1 #fmc is annual, need to make it monthly

    #unit fund
    units_p, fundv_p, fundch_p, fundb_p, fundbeg_p = calcUF(policy, prem_alloc, fmcp, unitPrices[runtype.runno,1,2,:])
    #output=table((time = 1:policy.projMax, alloc = prem_alloc, fundvalbeg = fundbeg_p, fundvalbef = fundb_p, fmc = fundch_p, fundvalend = fundv_p); pkey = [:time])
    #writedlm("bup2.csv", output, ',')

    non_units_p, noncf_p, nint_p, death_cost_p, nuexp_p, non_reserve_p = calcNUF(policy, (prem - prem_alloc), fundch_p, exp_init, expenses, aqxd, apx, fundv_p, unitPrices[runtype.runno,2,2,:])
    #output=table((time = 1:policy.projMax, alloc = (prem - prem_alloc), expenses = nuexp_p, interest = nint_p, fundmc = fundch_p, d_cost = death_cost_p, cf = noncf_p, v = non_reserve_p); pkey = [:time])
    #writedlm("nbup1.csv", output, ',')

    #TODO Profitvectors and signatures

    reserve = fundv[1]+non_reserve[2]  #2 because it is at start of year 1.
    return(reserve)
end

function Prem(policy::LifeValStructures.LifeAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy,runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    cf = deathClaims + expenses

    fixed, ren = perpolExp(policy,runtype, basisData)

    premium = (sum(cf .* tapx .* vt) .+ exp_init) ./ (sum((1 .- ren) .* tapx .* vt) .- fixed)
    return(premium[1])
end

function Prem(policy::LifeValStructures.Annuity, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa/12 .+ exp_claim) .* apx
    cf = deathClaims + expenses

    fixed, ren = perpolExp(policy,runtype, basisData)

    premium = (sum(cf .* tapx .* vt) .+ exp_init) ./ (1 .- fixed) #lump sum
    return(premium[1])
end

function Prem(policy::LifeValStructures.TermAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    cf = deathClaims + expenses

    fixed, ren = perpolExp(policy,runtype, basisData)

    premium = (sum(cf .* tapx .* vt) .+ exp_init) ./ (sum((1 .- ren) .* tapx .* vt) .- fixed)
    return(premium[1])
end

function Prem(policy::LifeValStructures.EndowmentAssurance, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    deathClaims = (sa .+ exp_claim) .* aqxd
    deathClaims[policy.projMax]=deathClaims[policy.projMax]/aqxd[policy.projMax]
    cf = deathClaims + expenses

    fixed, ren = perpolExp(policy,runtype, basisData)

    premium = (sum(cf .* tapx .* vt) .+ exp_init) ./ (sum((1 .- ren) .* tapx .* vt) .- fixed)

    return(premium[1])
end

function Prem(policy::LifeValStructures.EndowmentAssuranceUL, runtype, basisData, unitPrices)
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 0)

    prem_alloc, fmc = setAlloc(policy, runtype.runno, prem, runtype.fmc_sens, basisData, 0)

    #currently only 1 fund, 100% allocation
    fmcp = (1 .+ fmc) .^ (1/12) .- 1 #fmc is annual, need to make it monthly

    #unit fund
    units, fundv, fundch, fundb, fundbeg = calcUF(policy, prem_alloc, fmcp, unitPrices[runtype.runno,1,1,:])
    #output=table((time = 1:policy.projMax, alloc = prem_alloc, fundvalbeg = fundbeg, fundvalbef = fundb, fmc = fundch, fundvalend = fundv); pkey = [:time])
    #writedlm("bu2.csv", output, ',')

    #non-unit fund
    non_units, noncf, nint, death_cost, nuexp, non_reserve = calcNUF(policy, (prem - prem_alloc), fundch, exp_init, expenses, aqxd, apx, fundv, unitPrices[runtype.runno,2,1,:])
    #output=table((time = 1:policy.projMax, alloc = (prem - prem_alloc), expenses = nuexp, interest = nint, fundmc = fundch, d_cost = death_cost, cf = noncf, v = non_reserve); pkey = [:time])
    #writedlm("nbu1.csv", output, ',')

    #profit test BASIS
    aqxd, tapx, vt, apx, expenses, exp_init, exp_claim, prem, sa = setBasis(policy, runtype, basisData, 1)
    prem_alloc, fmc = setAlloc(policy, runtype.runno, prem, runtype.fmc_sens, basisData, 1)
    #currently only 1 fund, 100% allocation
    fmcp = (1 .+ fmc) .^ (1/12) .- 1 #fmc is annual, need to make it monthly

    #unit fund
    units_p, fundv_p, fundch_p, fundb_p, fundbeg_p = calcUF(policy, prem_alloc, fmcp, unitPrices[runtype.runno,1,2,:])
    #output=table((time = 1:policy.projMax, alloc = prem_alloc, fundvalbeg = fundbeg_p, fundvalbef = fundb_p, fmc = fundch_p, fundvalend = fundv_p); pkey = [:time])
    #writedlm("bup2.csv", output, ',')

    non_units_p, noncf_p, nint_p, death_cost_p, nuexp_p, non_reserve_p = calcNUF(policy, (prem - prem_alloc), fundch_p, exp_init, expenses, aqxd, apx, fundv_p, unitPrices[runtype.runno,2,2,:])
    #output=table((time = 1:policy.projMax, alloc = (prem - prem_alloc), expenses = nuexp_p, interest = nint_p, fundmc = fundch_p, d_cost = death_cost_p, cf = noncf_p, v = non_reserve_p); pkey = [:time])
    #writedlm("nbup1.csv", output, ',')

    #TODO Profitvectors and signatures

    #TODO calc prem
    premium = fundv[1]+non_reserve[2]  #2 because it is at end of year 1.
    return(premium)
end

end


text:
arraytype(A) = Tullio.storage_type(A) <: CuArray ? CuArray : Array
arraytype(::Type{T}) where {T} = Array
arraytype(::Type{<:CuArray}) = CuArray
arraytype(::Type{<:CUDA.Adaptor}) = CuArray

components(a::AbstractArray) = (a,)

function numbercontiguous(A; by=identity)
    p = sortperm(A; by=by)
    notequalprevious = fill!(similar(p, Bool), false)
    @tullio notequalprevious[i] =
        @inbounds(begin by(A[p[i]]) != by(A[p[i-1]]) end) (i in 2:length(p))

    B = similar(p)
    B[p] .= cumsum(notequalprevious) .+ 1

    return B
end

_numfields(::Type) = 1
_numfields(::Type{S}) where {S <: SArray} = length(S)
_numfields(T::Tuple) = sum(map(_numfields, T))
_numfields(N::NamedTuple) = sum(map(_numfields, N))

_fieldtype(::Type{S}) where {S} = eltype(S)
_fieldtype(T::Tuple) = promote_type(map(_fieldtype, T)...)
_fieldtype(N::NamedTuple) = promote_type(map(_fieldtype, N)...)

function fieldarray(::UndefInitializer, S, ::Type{A}, dims::Dims) where {A}
    N = _numfields(S)
    T = _fieldtype(S)

    if N == 1 && S isa Type && S <: Number
        return A{S}(undef, dims)
    end

    if length(dims) == 0
        d = (N, )
    elseif length(dims) == 1
        d = (dims[1], N)
    else
        d = (dims[1:end-1]..., N, dims[end])
    end

    data = A{T}(undef, d)
    dataviews = ntuple(N) do i
        offset = length(dims) > 1 ? 0 : 1
        viewtuple = ntuple(j->(j==length(dims)+offset ? i : :), length(d))
        return view(data, viewtuple...)
    end

    return fieldarray(S, dataviews)
end

function _ckfieldargs(S, data::Tuple)
    if _numfields(S) != length(data)
        throw(ArgumentError("Number of data fields is incorrect."))
    end
    T = _fieldtype(S)
    dims = size(first(data))
    for d in data
        if eltype(d) != T
            throw(ArgumentError("Data arrays do not have the same eltype."))
        end
        if size(d) != dims
            throw(ArgumentError("Data arrays do not have the same size."))
        end
    end
end

function fieldarray(::Type{S}, data::Tuple) where {S}
    d = only(data)
    if S != eltype(d)
        throw(ArgumentError("Data array does not have the correct eltype."))
    end

    return d
end

function fieldarray(::Type{S}, data::Tuple) where {S <: SArray}
    _ckfieldargs(S, data)
    return StructArray{S}(data)
end

function fieldarray(S::NamedTuple, data::Tuple)
    _ckfieldargs(S, data)

    offsets = cumsum((1, map(_numfields, S)...))
    fields = ntuple(length(S)) do i
        fieldarray(S[i], data[offsets[i]:offsets[i+1]-1])
    end

    return StructArray(NamedTuple{keys(S)}(fields))
end

ext: jl


text:
function main()
	N = 20  +1 # num grid points = num squares +1

	# at border, only one path possible
	# init grid with what is known
	grid = zeros(eltype(N),N,N)
	for i = 1:N
		grid[i,N] = 1
		grid[N,i] = 1
	end

	# fill middle of grid
	for i = N-1:-1:1, j = N-1:-1:1
		grid[i,j] = grid[i+1,j] + grid[i,j+1]
	end

	println(grid)
	println(grid[1,1])
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
mutable struct WallTimer
  starttime_ns::typeof(Base.time_ns())
  paused_elapsed_ns::typeof(Base.time_ns())

  WallTimer() = new(0,0)
end

function start!(timer::WallTimer)
  timer.starttime_ns = (Base.time_ns)()
  return nothing
end

started(timer::WallTimer) = (timer.starttime_ns ≠ 0)

""" Return nanoseconds since timer was started or 0 if not yet started. """
function elapsed(timer::WallTimer)
  local elapsedtime_ns = (Base.time_ns)() - timer.starttime_ns
  return started(timer) * elapsedtime_ns
end

function pause!(timer::WallTimer)
  timer.paused_elapsed_ns = (Base.time_ns)() - timer.starttime_ns
  return nothing
end

function unpause!(timer::WallTimer)
  timer.starttime_ns = (Base.time_ns)()
  timer.starttime_ns -= timer.paused_elapsed_ns;
  return nothing
end


text:
# Obtained from https://ccrma.stanford.edu/~jos/filters/Matlab_listing_tmps_m.html
spec = [1 1 1 0 0 0 1 1]'; # Lowpass cutting off at fs*3/8
@show mps(spec)
@show abs(mps(spec))
@show ifft(spec)
@show ifft(mps(spec))


text:
function define_arithmetic(machine)
    define(machine, "+",        gen_op_nn_n(+))
    define(machine, "-",        gen_op_nn_n(-))
    define(machine, "*",        gen_op_nn_n(*))
    define(machine, "/",        gen_op_nn_n(/))
    define(machine, "<",        gen_op_nn_n(less))
    define(machine, ">",        gen_op_nn_n(more))
    define(machine, "=",        gen_op_nn_n(equal))
end

less(a, b) = a < b ? -1 : 0
more(a, b) = less(b, a)
equal(a, b) = a == b ? -1 : 0

function gen_op_nn_n(operation)
    return (machine, parent, myidx) -> begin
        arg2 = pop!(machine.stack)
        arg1 = pop!(machine.stack)
        push!(machine.stack, operation(arg1, arg2))
        return 1
    end
end


text:
export lambert, benchmark_lambert

function benchmark_lambert(times::Int)
    r0 = [5000.0, 10000.0, 2100.0]
    r = [-14600.0, 2500.0, 7000.0]
    tof = 3600.0
    mu = 3.986004418e5
    worst = -Inf
    best = Inf
    total = 0.0
    lambert(mu, r0, r, tof)
    for i = 1:times
        gc_enable(false)
        t = @elapsed lambert(mu, r0, r, tof)
        gc_enable(true)
        if t > worst
            worst = t
        end
        if t < best
            best = t
        end
        total += t
    end
    println("[$(total/times),$best,$worst]")
end

function lambert(k::Float64, r0::Vector{Float64}, r::Vector{Float64}, tof::Float64, short::Bool=true, numiter::Int=35, rtol::Float64=1e-8)
    if short
        t_m = 1
    else
        t_m = -1
    end

    norm_r0 = norm(r0)
    norm_r = norm(r)
    cos_dnu = dot(r0, r) / (norm_r0 * norm_r)

    A = t_m * sqrt(norm_r * norm_r0 * (1 + cos_dnu))

    if A == 0.0
        error("Cannot compute orbit, phase angle is 180 degrees")
    end

    psi = 0.0
    psi_low = -4*pi
    psi_up = 4*pi

    count = 0
    converged = false
    y = 0.0
    while count < numiter
        y = norm_r0 + norm_r + A * (psi * c3(psi) - 1) / sqrt(c2(psi))
        if A > 0.0 && y < 0.0
            while y < 0.0
                psi_low = psi
                psi = (0.8 * (1.0 / c3(psi)) *
                    (1.0 - (norm_r0 + norm_r) * sqrt(c2(psi)) / A))
                y = norm_r0 + norm_r + A * (psi * c3(psi) - 1) / sqrt(c2(psi))
            end
        end
        xi = sqrt(y / c2(psi))
        tof_new = (xi^3 * c3(psi) + A * sqrt(y)) / sqrt(k)

        if abs((tof_new - tof) / tof) < rtol
            converged = true
            break
        else
            count += 1
            if tof_new <= tof
                psi_low = psi
            else
                psi_up = psi
            end
            psi = (psi_up + psi_low) / 2
        end
    end

    if !converged
        error("Maximum number of iterations reached")
    end

    f = 1 - y / norm_r0
    g = A * sqrt(y / k)
    gdot = 1 - y / norm_r

    v0 = (r - f * r0) / g
    v = (gdot * r - r0) / g
    return v0, v
end

function c2(psi::Float64)
    eps = 1.0
    if psi > eps
        res = (1 - cos(sqrt(psi))) / psi
    elseif psi < -eps
        res = (cosh(sqrt(-psi)) - 1) / (-psi)
    else
        res = 1.0 / 2.0
        delta = (-psi) / gamma(2 + 2 + 1)
        k = 1
        while res + delta != res
            res += delta
            k += 1
            delta = (-psi)^k / gamma(2*k + 2 + 1)
        end
    end
    return res
end

function c3(psi::Float64)
    eps = 1.0
    if psi > eps
        res = (sqrt(psi) - sin(sqrt(psi))) / (psi * sqrt(psi))
    elseif psi < -eps
        res = (sinh(sqrt(-psi)) - sqrt(-psi)) / (-psi * sqrt(-psi))
    else
        res = 1.0 / 6.0
        delta = (-psi) / gamma(2 + 3 + 1)
        k = 1
        while res + delta != res
            res += delta
            k += 1
            delta = (-psi)^k / gamma(2*k + 3 + 1)
        end
    end
    return res
end

max_issues_repo_name: OpenAstrodynamics/benchmarks

max_forks_repo_name: OpenAstrodynamics/benchmarks


text:
#module test_sparse

#export SparseMatrix

type SparseMatrix{T}
    index_l::Vector{Int}
    index_r::Vector{Int}
    values::Vector{T}
end

function SparseMatrix{T}(A::Matrix{T})
    index_l = Int[]
    index_r = Int[]
    values = T[]
    for j=1:size(A,2), i=1:size(A,1)
        if abs2(A[i,j]) > 1e-15
            push!(index_l, i)
            push!(index_r, j)
            push!(values, A[i,j])
        end
    end
    return SparseMatrix(index_l, index_r, values)
end

function perm!(S::SparseMatrix, perm)
    S.index_l = S.index_l[perm]
    S.index_r = S.index_r[perm]
    S.values = S.values[perm]
    return nothing
end

optimize_left!(S::SparseMatrix) = perm!(S, sortperm(S.index_l))
optimize_right!(S::SparseMatrix) = perm!(S, sortperm(S.index_r))

function inner_SB_js{T}(j::Int, index_l::Vector{Int}, index_r::Vector{Int}, values::Vector{T}, B::Matrix{T}, result::Matrix{T})
    @inbounds for s=1:length(index_l)
        result[index_l[s], j] += values[s] * B[index_r[s], j]
    end
    return nothing
end

function mul_SB_js{T}(S::SparseMatrix{T}, B::Matrix{T}, result::Matrix{T})
    @inbounds for j=1:size(B,1)
        inner_SB_js(j, S.index_l, S.index_r, S.values, B, result)
    end
    return nothing
end

function inner_SB_sj{T}(s::Int, index_l::Int, index_r::Int, value::T, B::Matrix{T}, result::Matrix{T})
    @inbounds for j=1:size(B,1)
        result[index_l, j] += value * B[index_r, j]
    end
    return nothing
end

function mul_SB_sj{T}(S::SparseMatrix{T}, B::Matrix{T}, result::Matrix{T})
    @inbounds for s=1:length(S.index_l)
        inner_SB_sj(s, S.index_l[s], S.index_r[s], S.values[s], B, result)
    end
    return nothing
end

function inner_BS_js{T}(i::Int, index_l::Vector{Int}, index_r::Vector{Int}, values::Vector{T}, B::Matrix{T}, result::Matrix{T})
    @inbounds for s=1:length(index_l)
        result[i, index_r[s]] += B[i, index_l[s]] * values[s]
    end
    return nothing
end

function mul_BS_js{T}(B::Matrix{T}, S::SparseMatrix{T}, result::Matrix{T})
    @inbounds for i=1:size(B,1)
        inner_BS_js(i, S.index_l, S.index_r, S.values, B, result)
    end
    return nothing
end

function inner_BS_sj{T}(s::Int, index_l::Int, index_r::Int, value::T, B::Matrix{T}, result::Matrix{T})
    @inbounds for i=1:size(B,1)
        result[i, index_r] += B[i, index_l] * value
    end
    return nothing
end

function mul_BS_sj{T}(B::Matrix{T}, S::SparseMatrix{T}, result::Matrix{T})
    @inbounds for s=1:length(S.index_l)
        inner_BS_sj(s, S.index_l[s], S.index_r[s], S.values[s], B, result)
    end
    return nothing
end

function doprofile(N, fill=0.01, seed=0)
    srand(seed)
    Sfull = full(sprand(N, N, fill))
    S = SparseMatrix(Sfull)
    B = rand(N,N)

    #BS = B*Sfull
    #SB = Sfull*B

    result = zeros(eltype(B), N, N)
    # Run once for jit
    mul_SB_js(S, B, result)
    mul_SB_sj(S, B, result)
    mul_BS_js(B, S, result)
    mul_BS_sj(B, S, result)

    T = Float64[]

    println("=========SxB=========")
    optimize_left!(S)
    fill!(result, 0)
    push!(T, @elapsed(mul_SB_js(S, B, result)))
    #println(norm(SB-result))
    fill!(result, 0)
    push!(T, @elapsed(mul_SB_sj(S, B, result)))
    #println(norm(SB-result))

    optimize_right!(S)
    fill!(result, 0)
    push!(T, @elapsed(mul_SB_js(S, B, result)))
    #println(norm(SB-result))
    fill!(result, 0)
    push!(T, @elapsed(mul_SB_sj(S, B, result)))
    #println(norm(SB-result))

    println("=========BxS=========")
    optimize_left!(S)
    fill!(result, 0)
    push!(T, @elapsed(mul_BS_js(B, S, result)))
    #println(norm(BS-result))
    fill!(result, 0)
    push!(T, @elapsed(mul_BS_sj(B, S, result)))
    #println(norm(BS-result))

    optimize_right!(S)
    fill!(result, 0)
    push!(T, @elapsed(mul_BS_js(B, S, result)))
    #println(norm(BS-result))
    fill!(result, 0)
    push!(T, @elapsed(mul_BS_sj(B, S, result)))
    #println(norm(BS-result))
    return T
end

N = [10:10:1000]
T = Array{Float64,2}[]
for n=N
    push!(T, doprofile(n, 0.1)')
end

T = vcat(T...)

using Winston

figure(height=2000,width=2000)
hold(true)

colors = ["y", "m", "c", "r", "g", "b", "w", "k"]

for i=1:size(T,2)
    println(T[:,i])
    semilogy(N, T[:,i], colors[i])
end
savefig("t.png")
#=
error()

function _f11{T}(j::Int, index_l::Vector{Int}, index_r::Vector{Int}, values::Vector{T}, alpha::T, B::Matrix{T}, result::Matrix{T})
    @inbounds for i=1:length(values)
        result[index_l[i], j] += alpha * values[i] * B[index_r[i], j]
    end
end

function gemm!{T}(alpha::T, M::SparseMatrix{T}, B::Matrix{T}, beta::T, result::Matrix{T})
    for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for j=1:size(B,1)
        _f11(j, M.index_l, M.index_r, M.values, alpha, B, result)
    end
end

function _f22{T}(j::Int, index_l::Vector{Int}, index_r::Vector{Int}, values::Vector{T}, alpha::T, B::Matrix{T}, result::Matrix{T})
    @inbounds for i=1:length(values)
        result[j, index_r[i]] += alpha * values[i] * B[j, index_l[i]]
    end
end

function gemm!{T}(alpha::T, B::Matrix{T}, M::SparseMatrix{T}, beta::T, result::Matrix{T})
    for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for j=1:size(B,1)
        _f22(j, M.index_l, M.index_r, M.values, alpha, B, result)
    end
end

# function gemm!{T}(alpha::T, B::Matrix{T}, M::SparseMatrix2{T}, beta::T, result::Matrix{T})
#     for j=1:size(result,2), i=1:size(result,1)
#         result[i,j] *= beta
#     end
#     @inbounds for i=1:length(M.values)
#         m = alpha*M.values[i]
#         i_l = M.index_l[i]
#         i_r = M.index_r[i]
#         for j=1:size(B,1)
#             result[i_l, j] += m * B[i_r, j]
#         end
#     end
# end

type SparseMatrix3{T}
    indices::Vector{Int}
    values::Vector{T}
end

function SparseMatrix3{T}(a::Matrix{T})
    indices = Int[]
    values = T[]
    for i=1:size(a,1), j=1:size(a,2)
        if abs2(a[i,j]) > 1e-15
            push!(indices, i)
            push!(indices, j)
            push!(values, a[i,j])
        end
    end

    return SparseMatrix3(indices, values)
end

function _f2r{T}(j::Int, indices::Vector{Int}, values::Vector{T}, alpha::T, B::Matrix{T}, result::Matrix{T})
    counter = 1
    @inbounds for i=1:length(values)
        result[indices[counter], j] += alpha * values[i] * B[indices[counter+1], j]
        counter += 2
    end
end

function gemm!{T}(alpha::T, M::SparseMatrix3{T}, B::Matrix{T}, beta::T, result::Matrix{T})
    for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for j=1:size(B,1)
        _f2r(j, M.indices, M.values, alpha, B, result)
    end
end

N = 2000
r = full(sprand(N, N, 0.1))
A = rand(N,N)

M1 = SparseMatrix(r)
M3 = SparseMatrix3(r)
M2 = SparseMatrix2(r)

println("version 1 reversed")
result = zeros(eltype(A), N, N)
@time gemm!(1., A, M1, 0., result)
result = zeros(eltype(A), N, N)
@time gemm!(1., A, M1, 0., result)
#println(norm(result-A*r))

println("version 2 reversed")
result = zeros(eltype(A), N, N)
@time gemm!(1., A, M2, 0., result)
result = zeros(eltype(A), N, N)
@time gemm!(1., A, M2, 0., result)
#println(norm(result-A*r))

println("version 1")
result = zeros(eltype(A), N, N)
@time gemm!(1., M1, A, 0., result)
result = zeros(eltype(A), N, N)
@time gemm!(1., M1, A, 0., result)
# println(norm(result-r*A))

println("version 2")
result = zeros(eltype(A), N, N)
@time gemm!(1., M2, A, 0., result)
result = zeros(eltype(A), N, N)
@time gemm!(1., M2, A, 0., result)
# println(norm(result-r*A))

println("version 3")
result = zeros(eltype(A), N, N)
@time gemm!(1., M3, A, 0., result)
result = zeros(eltype(A), N, N)
@time gemm!(1., M3, A, 0., result)
# println(norm(result-r*A))
=#

#end


text:
# ---
# title: 990. Satisfiability of Equality Equations
# id: problem990
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Union Find, Graph
# link: <https://leetcode.com/problems/satisfiability-of-equality-equations/description/>
# hidden: true
# ---
#
# Given an array equations of strings that represent relationships between
# variables, each string `equations[i]` has length `4` and takes one of two
# different forms: `"a==b"` or `"a!=b"`.  Here, `a` and `b` are lowercase
# letters (not necessarily different) that represent one-letter variable names.
#
# Return `true` if and only if it is possible to assign integers to variable
# names so as to satisfy all the given equations.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["a==b","b!=a"]
#     Output: false
#     Explanation: If we assign say, a = 1 and b = 1, then the first equation is satisfied, but not the second.  There is no way to assign the variables to satisfy both equations.
#
#
# **Example 2:**
#
#
#
#     Input: ["b==a","a==b"]
#     Output: true
#     Explanation: We could assign a = 1 and b = 1 to satisfy both equations.
#
#
# **Example 3:**
#
#
#
#     Input: ["a==b","b==c","a==c"]
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: ["a==b","b!=c","c==a"]
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: ["c==c","b==d","x!=z"]
#     Output: true
#
#
#
#
# **Note:**
#
#   1. `1 <= equations.length <= 500`
#   2. `equations[i].length == 4`
#   3. `equations[i][0]` and `equations[i][3]` are lowercase letters
#   4. `equations[i][1]` is either `'='` or `'!'`
#   5. `equations[i][2]` is `'='`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using GNSSTools
using SatelliteToolbox
using ProgressMeter
using Makie, FileIO, Colors, GeometryBasics, GLMakie
using AbstractPlotting, AbstractPlotting.MakieLayout

GM = 3.986004418e14  # m³s⁻²
eop = get_eop();
user_lla = (40.01, -105.2437, 1655);

###########ONE-WEB############
a = Rₑ+1200*1000;  # meters
plane_num = 18
sat_per_plane = 36;
incl = 87.9;  # degrees
ΔΩ = 360/plane_num/2;
###########IRIDIUM############
a = Rₑ+780*1000;  # meters
plane_num = 6;
sat_per_plane = 11;
incl = 86;
ΔΩ = 360/plane_num/2;
###########STARLINK###########
# a = Rₑ+1150*1000;  # meters
# plane_num = 32;
# sat_per_plane = 50;
# incl = 53;
# ΔΩ = 360/plane_num;
a_starlink = Rₑ+1150*1000;  # meters
#############GPS##############
# a = Rₑ+20000*1000;
# plane_num = 6;
# sat_per_plane = 4;
# incl = 56;
# ΔΩ = 360/plane_num;
##############################
orbital_period = 2π*sqrt(a^3 / GM)  # seconds
orbital_period_starlink = 2π*sqrt(a_starlink^3 / GM)  # seconds
# orb_period_fraction = 2
# orb_period_fraction = (60*60*24) / orbital_period
orb_period_fraction = 1
frames_per_second = 25  # frames/second
total_orbit_time = orb_period_fraction*orbital_period
total_orbit_time_starlink = orb_period_fraction*orbital_period_starlink
total_vid_time = 10  # seconds
# Δt = total_orbit_time/total_vid_time/frames_per_second  # seconds
Δt = total_orbit_time_starlink/total_vid_time/frames_per_second  # seconds
t_range = Array(0:Δt:total_orbit_time);
freq = 9650e6;
constellation = define_constellation(a, plane_num,
                                     sat_per_plane, incl,
                                     t_range; obs_lla=user_lla,
                                     print_steps=true,
                                     eop=eop,
                                     ΔΩ=ΔΩ,
                                     show_plot=false);
# Loop through each time step and produce a constellation plot to be combined
# into a video later
obs_ecef = GeodetictoECEF(deg2rad(user_lla[1]), deg2rad(user_lla[2]), user_lla[3])
N = length(constellation.t_range)
sat_num = length(constellation.satellites)
# Earth
u = Array(range(0, 2π, length=200)) .+ (180*pi/180)
v = Array(range(0, 1π, length=200))
u, v = meshgrid(u, v)
E_x = Rₑ.*cos.(u).*sin.(v)
E_y = Rₑ.*sin.(u).*sin.(v)
E_z = Rₑ.*cos.(v)

# earth = load(download("https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"));
earth = load("bluemarble-2048.png");

# function make_plot(i, scene=Scene(resolution=(1000, 1000), color=:black);
#                    sat2user_vectors=true)
#     pt_size = 9e8
#     strokewidth = 4
#     linewidth = 3
#     ϕ = deg2rad(270)
#     θ = deg2rad(65)
#     scaler = 3
#     camera_pos = (scaler*constellation.a).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]
#     light = Vec{3,Float32}[[1.0, 1.0, 1.0], [0.1, 0.1, 0.1], [0.9, 0.9, 0.9], [0.0, -20.0, 0.0]]
#     to_eci = rECEFtoECI(ITRF(), TEME(), constellation.epoch+constellation.t_range[i], eop)
#     # Transform Earth and CO from ECEF to ECI frames
#     E_xt = zeros(size(E_x))
#     E_yt = zeros(size(E_y))
#     E_zt = zeros(size(E_z))
#     for j in 1:size(E_x)[1]
#         for k in 1:size(E_x)[2]
#             x, y, z = to_eci*[E_x[j,k], E_y[j,k], E_z[j,k]]
#             E_xt[j,k] = x
#             E_yt[j,k] = y
#             E_zt[j,k] = z
#         end
#     end
#     surface!(E_xt, E_yt, E_zt, color=earth, shading=false, show_axis=false)
#     obs_eci = to_eci*[obs_ecef[1], obs_ecef[2], obs_ecef[3]]
#     for j in 1:sat_num
#         r_ecef = constellation.satellites[j].r_ecef[i,:]
#         r_eci = kepler_to_sv(constellation.satellites[j].orbit[i]).r
#         sat_range, az, el = calcelevation(r_ecef, user_lla)
#         scatter!([r_eci[1]], [r_eci[2]], [r_eci[3]],
#                 strokewidth=strokewidth, strokecolor=:grey)
#         if el > 10
#             lines!([obs_eci[1], r_eci[1]],
#                    [obs_eci[2], r_eci[2]],
#                    [obs_eci[3], r_eci[3]], color=:red, linewidth=linewidth)
#         end
#     end
#     for j in 1:sat_num
#         if (j%sat_per_plane) == (sat_per_plane-1)
#             r_ecis = zeros(size(constellation.satellites[j].r_ecef))
#             for k in 1:size(r_ecis)[1]
#                 # r_ecis[k,:] = to_eci*constellation.satellites[j].r_ecef[k,:]
#                 r_ecis[k,:] =  kepler_to_sv(constellation.satellites[j].orbit[k]).r
#             end
#             lines!(r_ecis[:,1], r_ecis[:,2], r_ecis[:,3], color=:grey, linewidth=linewidth)
#         end
#     end
#     update_cam!(scene, camera_pos, 0)
#     return scene
# end

function rotate_earth(i)
    to_eci = rECEFtoECI(ITRF(), TEME(), constellation.epoch+constellation.t_range[i], eop)
    # Transform Earth and CO from ECEF to ECI frames
    E_xt = zeros(size(E_x))
    E_yt = zeros(size(E_y))
    E_zt = zeros(size(E_z))
    for j in 1:size(E_x)[1]
        for k in 1:size(E_x)[2]
            x, y, z = to_eci*[E_x[j,k], E_y[j,k], E_z[j,k]]
            E_xt[j,k] = x
            E_yt[j,k] = y
            E_zt[j,k] = z
        end
    end
    return (E_xt, E_yt, E_zt)
end

function get_obj_positions(i)
    to_eci = rECEFtoECI(ITRF(), TEME(), constellation.epoch+constellation.t_range[i], eop)
    obs_eci = to_eci*obs_ecef
    sat_pos = zeros(length(constellation.satellites), 3)
    planes = []
    visible = []
    visible_sat = []
    for j in 1:length(constellation.satellites)
        r_ecef = constellation.satellites[j].r_ecef[i,:]
        r_eci = kepler_to_sv(constellation.satellites[j].orbit[i]).r
        sat_pos[j,:] = r_eci
        sat_range, az, el = calcelevation(r_ecef, user_lla)
        if el > 10
            push!(visible, [[obs_eci[1], r_eci[1]], [obs_eci[2], r_eci[2]], [obs_eci[3], r_eci[3]]])
            push!(visible_sat, r_eci)
        end
        if (constellation.satellites[j].id%sat_per_plane) == 1
            r_ecis = zeros(size(constellation.satellites[j].r_ecef))
            for k in 1:size(r_ecis)[1]
                r_ecis[k,:] =  kepler_to_sv(constellation.satellites[j].orbit[k]).r
            end
            push!(planes, r_ecis)
        end
    end
    visible_sats = zeros(length(sat_pos),3)
    for j in 1:size(sat_pos)[1]
        if j > length(visible_sat)
            visible_sats[j,:] = [NaN, NaN, NaN]
        else
            visible_sats[j,:] = visible_sat[j]
        end
    end
    return (sat_pos, planes, visible, visible_sats)
end

function initialize_figure(i; show_vectors=true, visible_color=missing)
    ϕ = deg2rad(180)
    θ = deg2rad(70)
    scaler = 3
    camera_pos = (scaler*constellation.a).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]
    light_ϕ = deg2rad(270)
    light_θ = deg2rad(90)
    light_scaler = 149597870700/Rₑ
    light_camera_pos = (light_scaler*constellation.a).*[cos(light_ϕ)*sin(light_θ), sin(light_ϕ)*sin(light_θ), cos(light_θ)]
    lightpos = Vec3f0(light_camera_pos...)
    pt_size = 9e8
    strokewidth = 3
    linewidth = 2
    linestyle = nothing
    pt_color = :lightgrey
    line_color = :lightgrey
    vector_color = :lightblue
    background_color = :white
    if ismissing(visible_color)
        visible_color = pt_color
    end
    E_xt, E_yt, E_zt = rotate_earth(i)
    sat_pos, planes, visible, visible_sats = get_obj_positions(i)
    scene = Scene(resolution=(1000, 1000), backgroundcolor=background_color)
    surf = surface!(scene, E_xt, E_yt, E_zt, shading=true, color=earth,
                    show_axis=false, lightposition=lightpos,
                    colormap=(:white, :white),
                    ambient = Vec3f0(0.5, 0.5, 0.5),
                    diffuse = Vec3f0(0.4, 0.4, 0.2),
                    specular = Vec3f0(0.5, 0.5, 0.5),
                    shininess = 0.5f0)
    scatter_sat = scatter!(scene, sat_pos[:,1], sat_pos[:,2],
    sat_pos[:,3],
                           strokewidth=strokewidth, strokecolor=pt_color)
    scatter_visible_sat = scatter!(scene,
                                   visible_sats[:,1],
                                   visible_sats[:,2],
                                   visible_sats[:,3],
                                   strokewidth=strokewidth,
                                   strokecolor=visible_color)
    for j in 1:length(planes)
        line_planes = lines!(scene, planes[j][:,1], planes[j][:,2],
                             planes[j][:,3],
                             color=line_color, linewidth=linewidth,
                             linestyle=linestyle)
    end
    if show_vectors
        val = [NaN, NaN]
        for j in 1:ceil(Int, length(constellation.satellites)/2)
            line_visible = lines!(scene, val, val,  val,
                                  color=vector_color, linewidth=linewidth)
        end
        for j in 1:length(visible)
            scene.plots[j+length(planes)+3][1] = visible[j][1]
            scene.plots[j+length(planes)+3][2] = visible[j][2]
            scene.plots[j+length(planes)+3][3] = visible[j][3]
        end
    end
    update_cam!(scene, camera_pos, 0)
    return scene
end

function modify_figure!(scene, i; show_vectors=true)
    ϕ = deg2rad(180)
    θ = deg2rad(70)
    scaler = 3
    camera_pos = (scaler*constellation.a).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]
    E_xt, E_yt, E_zt = rotate_earth(i)
    sat_pos, planes, visible, visible_sats = get_obj_positions(i)
    scene.plots[1][1] = E_xt
    scene.plots[1][2] = E_yt
    scene.plots[1][3] = E_zt
    scene.plots[2][1] = sat_pos[:,1]
    scene.plots[2][2] = sat_pos[:,2]
    scene.plots[2][3] = sat_pos[:,3]
    scene.plots[3][1] = visible_sats[:,1]
    scene.plots[3][2] = visible_sats[:,2]
    scene.plots[3][3] = visible_sats[:,3]
    for j in 1:length(planes)
        scene.plots[j+3][1] = planes[j][:,1]
        scene.plots[j+3][2] = planes[j][:,2]
        scene.plots[j+3][3] = planes[j][:,3]
    end
    if show_vectors
        val = [NaN, NaN]
        for j in 1:(length(scene.plots)-3-length(planes))
            scene.plots[j+length(planes)+3][1] = val
            scene.plots[j+length(planes)+3][2] = val
            scene.plots[j+length(planes)+3][3] = val

        end
        for j in 1:length(visible)
            scene.plots[j+length(planes)+3][1] = visible[j][1]
            scene.plots[j+length(planes)+3][2] = visible[j][2]
            scene.plots[j+length(planes)+3][3] = visible[j][3]
        end
    end
    update_cam!(scene, camera_pos, 0)
    return scene
end

function make_plots(N; show_vectors=true, visible_color=missing)
    scene = initialize_figure(1; show_vectors=show_vectors,
                              visible_color=visible_color)
    p = Progress(N, 1, "Generating $N frames...")
    compression = 10
    record(scene, "gps_constellation_black_compressed.mp4", 1:N; sleep=false,
           compression=compression, framerate=frames_per_second) do i
        modify_figure!(scene, i; show_vectors=show_vectors)
        next!(p)
    end
end

function zoom_earth(N=11*frames_per_second, start_i=9*frames_per_second,
                    end_i=10*frames_per_second;
                    show_vectors=true, visible_color=missing)
    a1 = Rₑ+20000*1000
    a2 = Rₑ+1150*1000;
    ϕ = deg2rad(180)
    θ = deg2rad(70)
    scaler = 3
    camera_pos = (scaler*a1).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]
    light_ϕ = deg2rad(270)
    light_θ = deg2rad(90)
    light_scaler = 149597870700/Rₑ
    light_camera_pos = (light_scaler*a1).*[cos(light_ϕ)*sin(light_θ), sin(light_ϕ)*sin(light_θ), cos(light_θ)]
    lightpos = Vec3f0(light_camera_pos...)
    E_xt, E_yt, E_zt = rotate_earth(1)
    scene = Scene(resolution=(1000, 1000), backgroundcolor=:black)
    surf = surface!(scene, E_xt, E_yt, E_zt, shading=true, color=earth,
                    show_axis=false, lightposition=lightpos,
                    colormap=(:white, :white),
                    ambient = Vec3f0(0.5, 0.5, 0.5),
                    diffuse = Vec3f0(0.4, 0.4, 0.2),
                    specular = Vec3f0(0.5, 0.5, 0.5),
                    shininess = 0.5f0)
    as = fill(a1, N)
    Δas = (a2 - a1)
    for i in start_i:end_i
        as[i] = a1 - Δas*(cos((i - start_i) * π / (end_i - start_i))/2 - 0.5)
    end
    as[end_i+1:end] .= a2
    p = Progress(N, 1, "Generating $N frames...")
    compression = 0
    record(scene, "zoom_earth.mp4", 1:N; sleep=false,
           compression=compression, framerate=frames_per_second) do i
        E_xt, E_yt, E_zt = rotate_earth(i)
        scene.plots[1][1] = E_xt
        scene.plots[1][2] = E_yt
        scene.plots[1][3] = E_zt
        camera_pos = (scaler*as[i]).*[cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]
        update_cam!(scene, camera_pos, 0)
        next!(p)
    end
end

function make_hist_plots()
    # GPS
    a = Rₑ+20000*1000;
    plane_num = 6;
    sat_per_plane = 4;
    incl = 56;
    ΔΩ = 360/plane_num;
    orbital_period = 2π*sqrt(a^3 / GM)  # seconds
    # orb_period_fraction = 2
    # orb_period_fraction = (60*60*24) / orbital_period
    orb_period_fraction = 10
    frames_per_second = 1  # frames/second
    total_orbit_time = orb_period_fraction*orbital_period
    total_vid_time = total_orbit_time
    # Δt = total_orbit_time/total_vid_time/frames_per_second  # seconds
    Δt = total_orbit_time/total_vid_time/frames_per_second  # seconds
    t_range = Array(0:Δt:total_orbit_time);
    # freq = 9650e6;
    freq = L1_freq
    doppler, doppler_rate, doppler_bounds,
    doppler_rate_bounds = doppler_distribution(a, plane_num,
                                               sat_per_plane, incl,
                                               t_range, user_lla,
                                               freq, show_plot=false,
                                               print_steps=true,
                                               eop=eop,
                                               ΔΩ=ΔΩ);
    # Make plot
    scene, layout = layoutscene(resolution=(840, 360))
    ax1 = layout[1, 1] = LAxis(scene, xlabel="Doppler Frequency (kHz)")
    hideydecorations!(ax1, label=true, ticklabels=true, ticks=true, grid=true)
    hidexdecorations!(ax1, label=false, ticklabels=false, ticks=false, grid=true)
    barplot!(ax1, doppler.edges[1][2:end]./1000, doppler.weights, color=:black,
            strokewidth=1, strokecolor=:black, show_axis=true)
    save("gps_hist_doppler.png", scene)
    scene, layout = layoutscene(resolution=(840, 360))
    ax1 = layout[1, 1] = LAxis(scene, xlabel = "Doppler Frequency Rate (Hz/s)")
    hideydecorations!(ax1, label=true, ticklabels=true, ticks=true, grid=true)
    hidexdecorations!(ax1, label=false, ticklabels=false, ticks=false, grid=true)
    barplot!(ax1, doppler_rate.edges[1][2:end], doppler_rate.weights, color=:black,
            strokewidth=1, strokecolor=:black, show_axis=true)
    save("gps_hist_doppler_rate.png", scene)
    # Starlink
    a = Rₑ+1150*1000;  # meters
    plane_num = 32;
    sat_per_plane = 50;
    incl = 53;
    ΔΩ = 360/plane_num;
    orbital_period = 2π*sqrt(a^3 / GM)  # seconds
    # orb_period_fraction = 2
    # orb_period_fraction = (60*60*24) / orbital_period
    orb_period_fraction = 10
    frames_per_second = 1  # frames/second
    total_orbit_time = orb_period_fraction*orbital_period
    total_vid_time = total_orbit_time
    # Δt = total_orbit_time/total_vid_time/frames_per_second  # seconds
    Δt = total_orbit_time/total_vid_time/frames_per_second  # seconds
    t_range = Array(0:Δt:total_orbit_time);
    # freq = 9650e6;
    freq = L1_freq
    doppler, doppler_rate, doppler_bounds,
    doppler_rate_bounds = doppler_distribution(a, plane_num,
                                               sat_per_plane, incl,
                                               t_range, user_lla,
                                               freq, show_plot=false,
                                               print_steps=true,
                                               eop=eop,
                                               ΔΩ=ΔΩ);
    # Make plot
    scene, layout = layoutscene(resolution=(840, 360))
    ax1 = layout[1, 1] = LAxis(scene, xlabel="Doppler Frequency (kHz)")
    hideydecorations!(ax1, label=true, ticklabels=true, ticks=true, grid=true)
    hidexdecorations!(ax1, label=false, ticklabels=false, ticks=false, grid=true)
    barplot!(ax1, doppler.edges[1][2:end]./1000, doppler.weights, color=:black,
            strokewidth=1, strokecolor=:black, show_axis=true)
    save("starlink_hist_doppler.png", scene)
    scene, layout = layoutscene(resolution=(840, 360))
    ax1 = layout[1, 1] = LAxis(scene, xlabel = "Doppler Frequency Rate (Hz/s)")
    hideydecorations!(ax1, label=true, ticklabels=true, ticks=true, grid=true)
    hidexdecorations!(ax1, label=false, ticklabels=false, ticks=false, grid=true)
    barplot!(ax1, doppler_rate.edges[1][2:end], doppler_rate.weights, color=:black,
            strokewidth=1, strokecolor=:black, show_axis=true)
    ax1.xticks = -200:50:0
    save("starlink_hist_doppler_rate.png", scene)
end
# run(`ffmpeg -r $frames_per_second -f image2 -i constellation_figures/constellation_%0$(pad_num)d.png -y constellation.mp4`)
run(`ffmpeg -i LabVideo_Sergei.mp4 -filter_complex "[0:v] palettegen" -y palette.png`)
run(`ffmpeg -i LabVideo_Sergei.mp4 -i palette.png -filter_complex "[0:v][1:v] paletteuse" -y LabVideo_Sergei.gif`)

# L1 C/A signal undergoing extreme doppler
using PyPlot
pygui(false)
t_length = 2
sig_freq = 9650e6;
sigtype = define_l1ca_code_type(t_length; sig_freq=sig_freq);
doppler_t = Array(0:0.001:t_length+0.1);
f_d = 150*1e3
fd_rate = -1000
fdd_rate = 175
fddd_rate = 0
doppler_curve = f_d .+ fd_rate .* doppler_t .+ fdd_rate .* doppler_t.^2 .+ fddd_rate .* doppler_t.^3;
signal = definesignal(sigtype, 5e6, t_length);
generatesignal!(signal; doppler_curve=doppler_curve, doppler_t=doppler_t);
acqresults, trackresults = process(signal, sigtype, 1; fd_center=f_d, show_plot=false);
# Iridium constellation @ X-band
t_range = Array(0:1:12*60*60);
a = Rₑ+780*1000;
plane_num = 6;
sat_per_plane = 11;
incl = 86;
freq = 9650e6;
ΔΩ = 30;
doppler, doppler_rate, doppler_bounds,
doppler_rate_bounds = doppler_distribution(a, plane_num,
                                           sat_per_plane, incl,
                                           t_range, user_lla,
                                           freq,
                                           print_steps=true,
                                           eop=eop,
                                           ΔΩ=ΔΩ;
                                           show_plot=false);
# Make animated plots
figsize = (17,8)
frames_per_second = 25
total_length = 10
N = floor(Int, frames_per_second*t_length)
pad_num = length(string(N))
ΔN = floor(Int, length(trackresults.t)/N)
p = Progress(N, 1, "Generating $N frames...")
for i in 0:N
    idx = i*ΔN
    fig = figure(figsize=figsize)
    ax1 = subplot(221)
    ax1.plot(trackresults.t[1:idx], trackresults.n0[1:idx].%1023, "k")
    xlim([trackresults.t[1], trackresults.t[end]])
    ylim([minimum(trackresults.n0.%1023), maximum(trackresults.n0.%1023)])
    title("Code Phase", fontsize=20)
    xticks([])
    yticks([])
    ax2 = subplot(222)
    ax2.plot(trackresults.t[1:idx], trackresults.dphi_meas[1:idx], "k")
    xlim([trackresults.t[1], trackresults.t[end]])
    ylim([minimum(trackresults.dphi_meas), maximum(trackresults.dphi_meas)])
    title("Phase Residuals", fontsize=20)
    xticks([])
    yticks([])
    ax3 = subplot(223)
    ax3.plot(trackresults.t[1:idx], trackresults.fds[1:idx], "k")
    xlim([trackresults.t[1], trackresults.t[end]])
    ylim([minimum(trackresults.fds), maximum(trackresults.fds)])
    title("Doppler", fontsize=20)
    xticks([])
    yticks([])
    ax4 = subplot(224)
    ax4.plot(trackresults.t[1:idx], real.(trackresults.ZP[1:idx]), "k")
    xlim([trackresults.t[1], trackresults.t[end]])
    ylim([minimum(real.(trackresults.ZP)), maximum(real.(trackresults.ZP))])
    # ax4.plot(trackresults.t, trackresults.data_bits, "k", linestyle="-")
    title("Navigation Data", fontsize=20)
    xticks([])
    yticks([])
    subplots_adjust(top=0.96, bottom=0.004, right=0.998, left=0.001,
                    hspace=0.2, wspace=0.1)
    # margins(0,0,0)
    # gca().xaxis.set_major_locator(PyPlot.plt.NullLocator())
    # gca().yaxis.set_major_locator(PyPlot.plt.NullLocator())
    if i == N
        num = string(0, pad=pad_num)
    else
        num = string(i+1, pad=pad_num)
    end
    savefig("data_figures/data_$num.png")
    close(fig)
    next!(p)
end
run(`ffmpeg -r $frames_per_second -f image2 -i data_figures/data_%0$(pad_num)d.png -y data.mp4`)
run(`ffmpeg -i data.mp4 -filter_complex "[0:v] palettegen" -y palette1.png`)
run(`ffmpeg -i data.mp4 -i palette1.png -filter_complex "[0:v][1:v] paletteuse" -y data.gif`)


text:
immutable StencilsnN{n_inner, n_outer} <: PartialDImplicitStrategy
  inner :: Int64
  outer :: Int64
  function StencilsnN()
    @assert eltype(n_inner) <: Int64
    @assert eltype(n_outer) <: Int64
    @assert n_inner <= n_outer
    new(n_inner, n_outer)
  end
end
export StencilsnN

# Defaultly set \Delta x = 1.0
# Defaultly set 3/3 Stencils, for simplify LU decompositon
function eval{n1}(
  pf :: PartialD_OP{PartialD{n1}, StencilsnN{3,3}, NaiveBoundary},
  data :: Array{Float64,1}; delta_x :: Float64 = 1.0
  )
  datax = initR(pf, data)
  (a,b,c)= initL(pf, data)
  (l,d,u) = naiveLU(a,b,c)
  return BackwardThomas(d,u, ForwardThomas(l,datax)) / ((delta_x)^(pf.pfunctor.n))
end

function eval{n1}(
  pf :: PartialD_OP{PartialD{n1}, StencilsnN{3,3}, NaivePeriodicShared},
  data_unshared :: Array{Float64,1}; delta_x :: Float64 = 1.0
  )
  n_shared = pf.bfunctor.n_shared
  data = cat((1,),data_unshared[end-n_shared+(1:n_shared)], data_unshared, data_unshared[1:n_shared] )
  datax = initR(pf, data)
  (a,b,c)= initL(pf, data)
  (l,d,u) = naiveLU(a,b,c)
  res =  BackwardThomas(d,u, ForwardThomas(l,datax)) / ((delta_x)^(pf.pfunctor.n))
  res[1+n_shared:end-n_shared]
end

include("initRL.jl")
include("naiveLU.jl")


text:
# deprecations should also include a comment stating the version
# for which the deprecation should be removed.


text:
#=
  Julia tokenizer for use with the Syntax tool

  https:#www.npmjs.com/package/syntax-cli

  See `--custom-tokenizer` to skip this generation, and use a custom one.

=#

#  Token: encapsulates token type, the matched value, and also location data.
Base.@kwdef mutable struct Token
    type::Int
    value::String
    startoffset = 1
    endoffset = 1
    startline = 1
    endline = 1
    startcolumn = 1
    endcolumn = 1
end

# TokenizerData: state of the tokenizer process
Base.@kwdef mutable struct TokenizerData
    tokensDict = {{{TOKENS}}}
    EOF_TOKEN = Token(type = tokensDict[EOF], value = EOF)
    lexRules = {{{LEX_RULES}}}
    lexRulesByConditionsDict = {{{LEX_RULES_BY_START_CONDITIONS}}}
    initstring
    states::Stack{String}
    cursor = 1
    tokensQueue::Queue{String}
    currentLine = 1
    currentColumn = 1
    currentLineBeginOffset = 1
    tokenStartoffset = 1
    tokenEndoffset = 1
    tokenStartline = 1
    tokenEndline = 1
    tokenStartcolumn = 1
    tokenEndcolumn = 1
end

# injected by parser
{{{LEX_RULE_HANDLERS}}}

#=
  Throws default "Unexpected token" exception, showing the actual
  line from the source, pointing with the ^ marker to the bad token.
  In addition, shows line:column location.
=#
function throwunexpectedtoken(tokenizerData::TokenizerData, symbol, line::Integer, column::Integer)
    throw(SyntaxError(string("\n\n", split(tokenizerData.initstring, "\n")[line], "\n", " "^(column - 1), "^\nUnexpected Token: \"", symbol, "\" at ", line, ":", column, ".")))
end

function inittokenizer(tokenizingString::AbstractString)
    mydata = TokenizerData(
        initstring = tokenizingString,
        states = Stack{String}(),
        tokensQueue = Queue{String}()
    )
    push!(mydata.states, "INITIAL")
    return mydata
end

#=
  All functions take the TokenizerData structure as a first value,
  similar to "self" pattern in other languages.

  Note that by convention in Julia, functions that change their
  arguments have ! which provides clarity to the user.
=#
function getcurrentstate(tokenizerData::TokenizerData)
    return first(tokenizerData.states)
end

function pushstate!(tokenizerData::TokenizerData, newState::AbstractString)
    push!(tokenizerData.states, newState)
end

function begin!(tokenizerData::TokenizerData, beginState::AbstractString)
    pushstate!(tokenizerData, beginState)
end

function popstate!(tokenizerData::TokenizerData)
    return pop!(tokenizerData.states)
end

function getnexttoken!(tokenizerData::TokenizerData)
    if !isempty(tokenizerData.tokensQueue)
        # process tokens waiting in the queue
        return totoken(tokenizerData, dequeue(tokenizerData.tokensQueue), "")
    end
    if !hasmoretokens(tokenizerData)
        return tokenizerData.EOF_TOKEN
    end

    ss = SubString(tokenizerData.initstring, tokenizerData.cursor)
    lexrulesforstate = tokenizerData.lexRulesByConditionsDict[getcurrentstate(tokenizerData)]
    # loop through all the lexer rules for this state to see what we can match
    for rulenum ∈ lexrulesforstate
        rule = tokenizerData.lexRules[rulenum+1]
        regexmatch = match(rule[1], ss)
        matchstr = ""
        if !isnothing(regexmatch)
            matchstr = regexmatch.match
            capturelocation(tokenizerData, matchstr)
            tokenizerData.cursor += length(matchstr)
        end
        # EOF token
        if length(ss) == 0 && !isnothing(regexmatch) && length(matchstr) == 0
            tokenizerData.cursor += 1
        end
        if !isnothing(regexmatch)
            global yytext = matchstr
            global yylength = length(matchstr)

            # the rules have strings that represent the names of functions to call
            rulefunction = getfield(SyntaxParser, Symbol(rule[2]))
            tokens = rulefunction()
            local token
            if isnothing(tokens)
                return getnexttoken!(tokenizerData)
            end
            if tokens isa AbstractVector
                token = tokens[1]
                for i in 2:length(tokens)
                    enqueue!(tokenizerData.tokensQueue, tokens[i])
                end
            else
                token = tokens
            end
            return totoken(tokenizerData, token, yytext)
        end
    end

    # If we are at the end of the file, push the cursor past the end and return that we have eaten EOF
    if isEOF(tokenizerData)
        tokenizerData.cursor += 1
        return tokenizerData.EOF_TOKEN
    end
    # we should not have reached here
    throwunexpectedtoken(tokenizerData, ss[1], tokenizerData.currentLine, tokenizerData.currentColumn)
end

# Given a string that matches a token, captures the location and start/end offsets
function capturelocation(tokenizerData::TokenizerData, matched::AbstractString)
    newline = r"\n"

    # absolute offsets
    tokenizerData.tokenStartoffset = tokenizerData.cursor

    # token start line-based locations
    tokenizerData.tokenStartline = tokenizerData.currentLine
    tokenizerData.tokenStartcolumn = tokenizerData.tokenStartoffset - tokenizerData.currentLineBeginOffset + 1

    # extract new line in the matched token
    for i in [x.offset for x in eachmatch(newline, matched)]
        tokenizerData.currentLine += 1
        tokenizerData.currentLineBeginOffset = tokenizerData.tokenStartoffset + i
    end
    tokenizerData.tokenEndoffset = tokenizerData.cursor + length(matched)

    # token end line-based locations
    tokenizerData.tokenEndline = tokenizerData.currentLine
    tokenizerData.tokenEndcolumn = tokenizerData.currentColumn = tokenizerData.tokenEndoffset - tokenizerData.currentLineBeginOffset
end

function totoken(tokenizerData::TokenizerData, tokenType::AbstractString, yytext::AbstractString)
    return Token(
        type = tokenizerData.tokensDict[tokenType],
        value = yytext,
        startoffset = tokenizerData.tokenStartoffset,
        endoffset = tokenizerData.tokenEndoffset,
        startline = tokenizerData.tokenStartline,
        endline = tokenizerData.tokenEndline,
        startcolumn = tokenizerData.tokenStartcolumn,
        endcolumn = tokenizerData.tokenEndcolumn
    )
end

function hasmoretokens(tokenizerData::TokenizerData)
    return tokenizerData.cursor <= length(tokenizerData.initstring)
end

function isEOF(tokenizerData::TokenizerData)
    return tokenizerData.cursor == (length(tokenizerData.initstring) + 1)
end

max_stars_repo_name: davidmclifton/syntax

max_issues_repo_name: davidmclifton/syntax

max_forks_repo_name: davidmclifton/syntax


text:
using Printf
using ShaperOriginDesignLib
using Unitful
using UnitfulUS
using Test

@testset "ShaperOriginDesignLib.jl" begin
    @test svgval(1.5) == 1.5
    @test svgval(4.5u"inch") == 4.5
    d = pathd(
        ["M", 0, 0],
        ["h", 3u"cm"],
        ["v", 4u"cm"],
        ["z"]
    )
    @test d == @sprintf("M 0 0 h %3f v %3f z",
                        3 / 2.54, 4 / 2.54)
    @test shaper_style_string(:on_line_cut) ==
        "fill: none; stroke: gray; stroke-width: 0.01; opacity: 1.0;"
end

@testset "viewport_attributes" begin
    left, top, right, bottom = (0.0u"inch", 0.0u"cm", 5.0u"inch", 2.54u"cm")
    attrs = viewport_attributes(left, top, right, bottom, u"inch")
    @test attrs.width == "5.0in"
    @test attrs.height == "1.0in"
    @test attrs.viewBox == "0.0 0.0 5.0 1.0"
end

max_stars_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_issues_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_forks_repo_name: MarkNahabedian/ShaperOriginDesignLib


text:
function ChainRulesCore.rrule(::typeof(eachrow), X::VecOrMat)
    eachrow_pullback(ΔΩ::Tangent) = (NoTangent(), ΔΩ.f.A)
    return eachrow(X), eachrow_pullback
end

function ChainRulesCore.rrule(::typeof(inv), D::Diagonal{<:Real})
    Ω = inv(D)
    function inv_Diagonal_pullback(ΔΩ::NamedTuple{(:diag,)})
        return (NoTangent(), (diag = .-ΔΩ.diag .* Ω.diag .^2,))
    end
    function inv_Diagonal_pullback(ΔΩ::Diagonal)
        return (NoTangent(), Diagonal(.-ΔΩ.diag .* Ω.diag .^2))
    end
    return Ω, inv_Diagonal_pullback
end

function ChainRulesCore.rrule(::typeof(inv), C::Cholesky{<:BLAS.BlasFloat, <:StridedMatrix})
    Ω = inv(C)
    function inv_Cholesky_pullback(ΔΩ::StridedMatrix{<:BLAS.BlasFloat})
        return (NoTangent(), (factors = -C.U' \ (ΔΩ + ΔΩ') * Ω, ))
    end
    function inv_Cholesky_pullback(ΔΩ::AbstractMatrix{<:BLAS.BlasFloat})
        return inv_Cholesky_pullback(collect(ΔΩ))
    end
    return Ω, inv_Cholesky_pullback
end


text:
using StatsBase
using LinearAlgebra, Random, Test

@testset "StatsBase.Histogram" begin

@testset "Histogram binindex and binvolume" begin
    edg1 = -2:0.5:9
    edg1f0 = -2:0.5f0:9
    edg2 = [-2, -1, 2, 7, 19]
    h1 = Histogram(edg1)
    h2 = Histogram((edg1, edg2))
    h3 = Histogram((edg1f0, edg2))

    @test h1 == Histogram(edg1, :left, false)

    @test @inferred StatsBase.binindex(h1, -0.5) == 4
    @test @inferred StatsBase.binindex(h2, (1.5, 2)) == (8, 3)

    @test [StatsBase.binvolume(h1, i) for i in axes(h1.weights, 1)] ≈ diff(edg1)
    @test [StatsBase.binvolume(h2, (i,j)) for i in axes(h2.weights, 1), j in axes(h2.weights, 2)] ≈ diff(edg1) * diff(edg2)'

    @test typeof(@inferred(StatsBase.binvolume(h2, (1,1)))) == Float64
    @test typeof(@inferred(StatsBase.binvolume(h3, (1,1)))) == Float32
    @test typeof(@inferred(StatsBase.binvolume(Float64, h3, (1,1)))) == Float64
end

@testset "Histogram append" begin
    h = Histogram(0:20:100, Float64, :left, false)
    @test @inferred(append!(h, 0:0.5:99.99)) == h
    @test append!(Histogram(0:20:100, Float64, :left, false), 0:0.5:99.99).weights ≈ [40,40,40,40,40]
    @test append!(Histogram(0:20:100, Float64, :left, true), 0:0.5:99.99).weights ≈ [2,2,2,2,2]
    @test append!(Histogram(0:20:100, Float64, :left, false), 0:0.5:99.99, fill(2, 200)).weights ≈ [80,80,80,80,80]
    @test append!(Histogram(0:20:100, Float64, :left, true), 0:0.5:99.99, fill(2, 200)).weights ≈ [4,4,4,4,4]
end

@testset "Histogram fit" begin
    @test sum(fit(Histogram,[1,2,3]).weights) == 3
    @test fit(Histogram,Int[]).weights == Int[]
    @test fit(Histogram,[1]).weights == [1]
    @test fit(Histogram,[1,2,3],[0,2,4]) == Histogram([0,2,4],[1,2], :left)
    @test fit(Histogram,[1,2,3],[0,2,4]) != Histogram([0,2,4],[1,1], :left)
    @test fit(Histogram,[1,2,3],0:2:4) == Histogram(0:2:4,[1,2], :left)
    @test all(fit(Histogram,[0:99;]/100,0.0:0.01:1.0).weights .==1)
    @test fit(Histogram,[1,1,1,1,1]).weights[1] == 5
    @test sum(fit(Histogram,(rand(100),rand(100))).weights) == 100
    @test fit(Histogram,1:100,nbins=5,closed=:right).weights == [20,20,20,20,20]
    @test fit(Histogram,1:100,nbins=5,closed=:left).weights == [19,20,20,20,20,1]
    @test fit(Histogram,0:99,nbins=5,closed=:right).weights == [1,20,20,20,20,19]
    @test fit(Histogram,0:99,nbins=5,closed=:left).weights == [20,20,20,20,20]

    @test fit(Histogram,(0:99,0:99),nbins=5).weights == Matrix(Diagonal([20,20,20,20,20]))
    @test fit(Histogram,(0:99,0:99),nbins=(5,5)).weights == Matrix(Diagonal([20,20,20,20,20]))

    @test fit(Histogram,0:99,weights(ones(100)),nbins=5).weights == [20,20,20,20,20]
    @test fit(Histogram,0:99,weights(2*ones(100)),nbins=5).weights == [40,40,40,40,40]
    @test fit(Histogram{Int32},0:99,weights(2*ones(100)),nbins=5).weights == [40,40,40,40,40]
    @test fit(Histogram{Float32},0:99,weights(2*ones(100)),nbins=5).weights == [40,40,40,40,40]
end

@testset "Histogram element type" begin
    @test eltype(@inferred(fit(Histogram,1:100,weights(ones(Int,100)),nbins=5)).weights) == Int
    @test eltype(@inferred(fit(Histogram{Float32},1:100,weights(ones(Int,100)),nbins=5)).weights) == Float32
    @test eltype(@inferred(fit(Histogram,1:100,weights(ones(Float64,100)),nbins=5)).weights) == Float64
    @test eltype(@inferred(fit(Histogram{Float32},1:100,weights(ones(Float64,100)),nbins=5)).weights) == Float32
end

@testset "histrange" begin
    # Note: atm histrange must be qualified
    @test @inferred(StatsBase.histrange(Float64[], 0, :left)) == 0.0:1.0:0.0
    @test StatsBase.histrange(Float64[1:5;], 1, :left) == 0.0:5.0:10.0
    @test StatsBase.histrange(Float64[1:10;], 1, :left) == 0.0:10.0:20.0
    @test StatsBase.histrange(1.0, 10.0, 1, :left) == 0.0:10.0:20.0

    @test StatsBase.histrange([0.201,0.299], 10, :left) == 0.2:0.01:0.3
    @test StatsBase.histrange([0.2,0.299], 10, :left) == 0.2:0.01:0.3
    @test StatsBase.histrange([0.2,0.3], 10, :left)  == 0.2:0.01:0.31
    @test StatsBase.histrange(0.2, 0.3,  10, :left)  == 0.2:0.01:0.31
    @test StatsBase.histrange([0.2,0.3], 10, :right) == 0.19:0.01:0.3
    @test StatsBase.histrange(0.2, 0.3,  10, :right) == 0.19:0.01:0.3

    @test StatsBase.histrange([200.1,299.9], 10, :left) == 200.0:10.0:300.0
    @test StatsBase.histrange([200.0,299.9], 10, :left) == 200.0:10.0:300.0
    @test StatsBase.histrange([200.0,300.0], 10, :left) == 200.0:10.0:310.0
    @test StatsBase.histrange([200.0,300.0], 10, :right) == 190.0:10.0:300.0

    @test @inferred(StatsBase.histrange(Int64[1:5;], 1, :left)) == 0:5:10
    @test StatsBase.histrange(Int64[1:10;], 1, :left) == 0:10:20

    @test StatsBase.histrange([0, 1, 2, 3], 4, :left) == 0.0:1.0:4.0
    @test StatsBase.histrange([0, 1, 1, 3], 4, :left) == 0.0:1.0:4.0
    @test StatsBase.histrange([0, 9], 4, :left) == 0.0:5.0:10.0
    @test StatsBase.histrange([0, 19], 4, :left) == 0.0:5.0:20.0
    @test StatsBase.histrange([0, 599], 4, :left) == 0.0:200.0:600.0
    @test StatsBase.histrange([-1, -1000], 4, :left) == -1000.0:500.0:0.0

    # Base issue #13326
    l,h = extrema(StatsBase.histrange([typemin(Int),typemax(Int)], 10, :left))
    @test l <= typemin(Int)
    @test h >= typemax(Int)

    @test_throws ArgumentError StatsBase.histrange([1, 10], 0, :left)
    @test_throws ArgumentError StatsBase.histrange([1, 10], -1, :left)
    @test_throws ArgumentError StatsBase.histrange([1.0, 10.0], 0, :left)
    @test_throws ArgumentError StatsBase.histrange([1.0, 10.0], -1, :left)
    @test_throws ArgumentError StatsBase.histrange(Float64[],-1, :left)
    @test_throws ArgumentError StatsBase.histrange([0.], 0, :left)
end

@testset "Histogram show" begin
    # hist show
    show_h = sprint(show, fit(Histogram,[0,1,2]))
    @test occursin("edges:\n  0.0:1.0:3.0", show_h)
    @test occursin("weights: $([1,1,1])", show_h)
    @test occursin("closed: left", show_h)
    @test occursin("isdensity: false", show_h)
end

@testset "Histogram norm and normalize" begin
    rng = MersenneTwister(345678)
    edges = (
        cumsum(rand(rng) * rand(rng, 9)),
        cumsum(rand(rng, 1:10) * rand(rng, 1:100, 11)),
        cumsum(5 * rand(rng) * rand(rng, 14))
    )

    n = 100000

    data = (
        maximum(edges[1]) .* (randn(rng, n) ./ 6 .+ 0.5),
        rand(rng, 1:maximum(edges[2]), n),
        maximum(edges[3]) .* rand(rng, n)
    )

    h = fit(Histogram, data, edges, closed = :left)

    weight_sum = sum(h.weights)
    bin_vols = [ x * y * z for x in diff(edges[1]), y in diff(edges[2]), z in diff(edges[3])]

    @test norm(h) ≈ sum(h.weights .* bin_vols)

    @test @inferred(normalize(h, mode = :none)) == h

    h_pdf = normalize(h, mode = :pdf)
    @test h_pdf.weights ≈ h.weights ./ bin_vols ./ weight_sum
    @test h_pdf.isdensity == true
    @test @inferred(norm(h_pdf)) ≈ 1
#    @test @inferred(normalize(h_pdf, mode = :pdf)) == h_pdf
    @test @inferred(normalize(h_pdf, mode = :density)) == h_pdf
#    @test @inferred(normalize(h_pdf, mode = :probability)) == h_pdf

    h_density = normalize(h, mode = :density)
    @test h_density.weights ≈ h.weights ./ bin_vols
    @test h_density.isdensity == true
    @test @inferred(norm(h_density)) ≈ weight_sum
    @test @inferred(normalize(h_density, mode = :pdf)) ==
        Histogram(h_density.edges, h_density.weights .* (1/norm(h_density)), h_density.closed, true)
    @test normalize(h_density, mode = :pdf).weights ≈ h_pdf.weights
    @test normalize(h_density, mode = :density) == h_density
    @test normalize(h_density, mode = :probability).weights ≈ h_pdf.weights

    h_fraction = normalize(h, mode = :probability)
    @test sum(h_fraction.weights) ≈ 1
    @test h_fraction.isdensity == false
    @test normalize(h_fraction, mode = :pdf).weights ≈ h_pdf.weights
    @test normalize(h_fraction, mode = :density).weights ≈ h_pdf.weights
    @test normalize(h_fraction, mode = :probability).weights ≈ h_fraction.weights

    h_copy = deepcopy(float(h))
    @test @inferred(normalize!(h_copy, mode = :density)) == h_copy

    h2 = deepcopy(float(h))
    mod_h2 = normalize!(h2, mode = :density)
    @test mod_h2 === h2 && mod_h2.weights === h2.weights
    @test h2.weights == h_density.weights

    aux_weights = sqrt.(h.weights)
    divor0 = (a,b) -> (a == 0 && b == 0) ? 0 : a/b
    divor0_cmp = (a_n, a_d, b_n, b_d) -> maximum(abs.(map(divor0, a_n, a_d) - map(divor0, b_n, b_d))) < 1e-10

    h_pdf2, h_pdf2_aux = normalize(float(h), aux_weights, mode = :pdf)
    @test divor0_cmp(h_pdf2_aux, aux_weights, h_pdf2.weights, h.weights)

    h_density2, h_density2_aux = normalize(float(h), aux_weights, mode = :density)
    @test divor0_cmp(h_density2_aux, aux_weights, h_density2.weights, h.weights)

    h_density3, h_density3_aux = normalize(h_density2, h_density2_aux, mode = :pdf)
    @test divor0_cmp(h_density3_aux, h_density2_aux, h_density3.weights, h_density2.weights)
end

@testset "Histogram zero" begin
    h = fit(Histogram, (rand(100), rand(100)))
    h2 = @inferred zero(h)
    @test all(x -> x≈0, h2.weights)
    @test !(h.weights === h2.weights)
    @test h.edges == h2.edges
    @test h.closed == h2.closed
    @test h.isdensity == h2.isdensity
end

@testset "Histogram merge" begin
    histograms = [fit(Histogram, (rand(100), 10 * rand(100)), (0:0.1:1, 0:1:10)) for _ in 1:10]
    h = zero(histograms[1])
    merge!(h, histograms ...)
    @test h.weights == (+).((x->x.weights).(histograms)...)
    @test (@inferred merge(histograms...)) == h
end

@testset "midpoints" begin
    @test StatsBase.midpoints([1, 2, 4]) == [1.5, 3.0]
    @test StatsBase.midpoints(range(0, stop = 1, length = 5)) == 0.125:0.25:0.875
end

end # @testset "StatsBase.Histogram"


text:
using Test

import MutableArithmetics
const MA = MutableArithmetics

using DynamicPolynomials

@testset "Mutable Arithmetics with $T" for T in [Int, BigInt]
    @polyvar x y z
    p = T(2)x^2 + T(4)y^2
    q = T(3)x * y
    @test p === MA.add!(p, q)
    @test q == T(3)x * y
    @test nterms(p) == 3
    @test issorted(monomials(p), rev=true)
    @test p == T(2)x^2 + T(3)x * y + T(4)y^2
    @test p === MA.add!(p, x * y)
    @test nterms(p) == 3
    @test issorted(monomials(p), rev=true)
    @test p == T(2)x^2 + T(4)x * y + T(4)y^2
    @test p === MA.add!(p, T(2)x^3 + T(3)x^2 + T(5)x + T(2))
    @test nterms(p) == 6
    @test issorted(monomials(p), rev=true)
    @test p == T(2)x^3 + T(5)x^2 + T(4)x * y + T(4)y^2 + T(5)x + T(2)
    q = T(5)z^3 + T(1)x*z + T(4)z^2 + T(3)z
    @test p === MA.add!(p, q)
    @test nterms(p) == 10
    @test issorted(monomials(p), rev=true)
    @test p == T(2)x^3 + T(5)x^2 + T(4)x * y + T(4)y^2 + T(5)x + T(2) + q

    @testset "Issue #62" begin
        @polyvar x y
        p = x^2 + x + 1
        q = rem(p, [x^2-y])
        @test q == x + y + 1
    end
end


text:
include("../src/modeling/util.jl")

function model(t :: Trace, data)
    loc = sample(t, "loc", Normal(0.0, 10.0))
    scale = sample(t, "scale", Gamma(3.0, 3.0))
    for (i, d) in enumerate(data)
        observe(t, "data $i", Normal(loc, scale), d)
    end
end

@testset "nested sampling 1" begin
    true_loc = 8.0
    true_scale = 2.0
    data = rand(4) .* true_scale .+ true_loc
    @info "True loc = $true_loc"
    @info "True scale = $true_scale"

    for n in [3, 4, 8, 15, 25, 50]
        @info "Nested sampling inference with $n points"
        @time results = nested(model; params = (data,), num_points = n)
        mean_loc = mean(results, "loc")
        mean_std = mean(results, "scale")
        @info "Posterior mean loc = $mean_loc"
        @info "Posterior mean scale = $mean_std"
        log_Z = logsumexp(results.log_weights)
        @info "log p(x) ≈ $log_Z"
    end
end