text:
# define  struct to hold resulting data
struct ModelStats
    comparisonByAY::DataFrame
    comparisonByLOB::Dict
    absErrorsByAY::Array
    comparisons::Dict
end

function buildTriangle(data;rowIdentifier=:AY,columns=[Symbol(string("PayCum", lpad(string(i), 2, "0"))) for i = 0:11])
    res = DataFrames.aggregate(data[vcat(columns, rowIdentifier)], rowIdentifier, sum);
    res2 = convert(Array, res)
    return res2[:,2:end]
end

function buildTriangleOLD(data;rowIdentifier=:AY,columns=[Symbol(string("PayCum", lpad(string(i), 2, "0"))) for i = 0:11])
    minRow, MaxRow = extrema(data[rowIdentifier])
    nRows = MaxRow - minRow + 1
    nColumns = length(columns)
    res = zeros(Int, nRows, nColumns)
    have = names(data)
    @assert issubset(columns, have)
    for i = 1:size(data, 1)
        for j = 1:length(columns)
            @inbounds thisAY = data[rowIdentifier][i]
            @inbounds res[thisAY - minRow + 1,j] += data[columns[j]][i]
        end
    end
    return res
end

function chainLadder(x;tail=1.0)
    # x is the cumulative triangle
    k = 1
    height = size(x, 1)
    width = size(x, 2)
    LDFs = zeros(width)
    factorsToUltimate = ones(width)
    ultimate = zeros(height)
    paidToDate = zeros(height)
    for j = 1:size(x, 2) - 1
        LDFs[j] = sum(x[1:height - j,j + 1]) / sum(x[1:height - j,j])
    end
    LDFs[end] = tail
    factorsToUltimate[end] = LDFs[end]
    for j = size(x, 2):-1:2
        factorsToUltimate[j - 1] *= (factorsToUltimate[j] * LDFs[j - 1])
        ultimate[height - j + 2] = x[height - j + 2,j - 1] * factorsToUltimate[j - 1]
        paidToDate[height - j + 2] = x[height - j + 2,j - 1]
    end
    j = 1 + size(x, 2)
    ultimate[height - j + 2] = x[height - j + 2,j - 1] * factorsToUltimate[j - 1]
    paidToDate[height - j + 2] = x[height - j + 2,j - 1]
    reserves = ultimate .- paidToDate
    # construct output DataFrame
    reverse!(LDFs)
    reverse!(factorsToUltimate)
    df = DataFrame(LDFS=LDFs, factorsToUltimate=factorsToUltimate, paidToDate=paidToDate, reserves=reserves, ultimate=ultimate)
    return df
end

function aggregateReservesOfIndividualCLModels(fullData, leafNrs)
    reservesCombined = zeros(length(unique(fullData[:AY])))
    uqLeafNrs = sort(unique(leafNrs))
    for leaf in uqLeafNrs
        partialData = fullData[leafNrs .== leaf,:]
        triangle = buildTriangle(partialData)
        resDF = chainLadder(triangle)
        reservesCombined .+= resDF[:reserves]
    end
    return reservesCombined
end

function calculateEstimates(fullData, LDFArray, aggCL, paidToDatePerRow)
    ayears = unique(fullData[:AY])
    # apply Chainladder to each row
    ultimatePerRow = zeros(size(fullData, 1))
    cumulativefactors = zeros(size(fullData, 1))
    ultimatePerAY = zeros(length(ayears))
    truthPerAY = zeros(length(ayears))
    paycumcols = [Symbol(string("PayCum", lpad(string(i), 2, "0"))) for i = 0:11]
    minAY, maxAY = extrema(ayears)
    for iRow = 1:size(fullData, 1)
        @inbounds row = fullData[:AY][iRow] - minAY + 1
        @inbounds rowInversed = length(ayears) - row + 1
        @inbounds cumulativefactor = prod(1 ./ LDFArray[iRow,rowInversed:end])
        @inbounds cumulativefactors[iRow] = cumulativefactor
        @inbounds ultimatePerRow[iRow] = cumulativefactor * paidToDatePerRow[iRow]
        @inbounds ultimatePerAY[row] += ultimatePerRow[iRow]
        @inbounds truthPerAY[row] += fullData[:PayCum11][iRow]
    end
    reservesPerAY = ultimatePerAY .- aggCL[:paidToDate]
    errPerAY = ultimatePerAY .- truthPerAY
    dfPerRow = DataFrame(paidToDate=paidToDatePerRow, reserves=ultimatePerRow .- paidToDatePerRow, ultimate=ultimatePerRow, cumulativefactors=cumulativefactors)
    df = DataFrame(paidToDate=deepcopy(aggCL[:paidToDate]), reserves=reservesPerAY, ultimate=ultimatePerAY, truth=truthPerAY, error=errPerAY)
    return dfPerRow, df
end

function getPaidToDatePerRow(fullData)
    ayears = unique(fullData[:AY])
    paidToDatePerRow = zeros(size(fullData, 1))
    minAY, maxAY = extrema(ayears)
    paycumcols = [Symbol(string("PayCum", lpad(string(i), 2, "0"))) for i = 0:11]
    for iRow = 1:size(fullData, 1)
        @inbounds row = fullData[:AY][iRow] - minAY + 1
        @inbounds rowInversed = length(ayears) - row + 1
        @inbounds paidToDatePerRow[iRow] = fullData[paycumcols[rowInversed]][iRow]
    end
    return paidToDatePerRow
end

function getcorrectedQuantilesOfError(estPerRow, truthPerRow, qtl_range, CLTotalUltimate)
    errPerRow = estPerRow[:ultimate] - truthPerRow
    totalUltimate = sum(estPerRow[:ultimate])
    # what if we correct the ultimate total to match the CL total?
    corrFactor = 1 / totalUltimate * CLTotalUltimate
    @assert isapprox(0, sum(corrFactor .* estPerRow[:ultimate]) - CLTotalUltimate, atol=1e-6)
    errPerRowCorrected = corrFactor .* estPerRow[:ultimate] - truthPerRow
    return Statistics.quantile(errPerRowCorrected, qtl_range)
end

function checkIfPaymentsAreIncreasing(data)
    # consider payment columns
    paymCols = [Symbol(string("PayCum", lpad(ii, 2, "0"))) for ii = 0:11]
    isNonDecreasing = BitVector(undef, size(data, 1))
    data2 = data[paymCols]
    dataAsArray = transpose(convert(Array, data2))
    for ii = 1:size(data2, 1)
        @inbounds isNonDecreasing[ii] = issorted(dataAsArray[:,ii])
    end
    return isNonDecreasing
end

function customSummary(treeEstimateAgg::DataFrame, treeEstimatePerRow::DataFrame, trueReservesPerLOB, clPerLOB;writeResultToTemp=true)
    # NOTE: this function is using several global variables!
    # overall Result
    comparisonByAY = hcat(ayears, trueReserves, treeEstimateAgg[:reserves], clAllLOBs[:reserves])
    comparisonByAY = vcat(comparisonByAY, sum(comparisonByAY, dims=1))
    comparisonByAY[end,1] = 9999
    comparisonByAY2 = DataFrame(comparisonByAY, :auto)
    names!(comparisonByAY2, [:AY,Symbol("True Outstanding Amount"),Symbol("Estimated Reserves"),Symbol("CL Estimated Reserves")])

    absErrorsByAY = ceil.(Int, abs.(comparisonByAY[:,2] .- comparisonByAY[:,3]) / 1000)
    @show absErrorsByAY
    @show absErrorsByAY[end]

    # result by LOB versus four CL  models (one per LOB)
    treeEstimatePerRow[:AY] = dataKnownByYE2005[:AY]
    comparisonByLOB = Dict{String,DataFrame}()
    for lob in sort(unique(dataKnownByYE2005[:LoB]))
        sumRes = DataFrames.aggregate(treeEstimatePerRow[dataKnownByYE2005[:LoB] .== lob,[:reserves,:AY]], :AY, sum)
        sumRes2 = hcat(ayears, trueReservesPerLOB[lob][1:end - 1], sumRes[:reserves_sum], clPerLOB[lob][:reserves][1:end - 1])
        sumRes2 = vcat(sumRes2, sum(sumRes2, dims=1))
        sumRes2[end,1] = 9999
        sumRes3 = DataFrame(sumRes2, :auto)
        names!(sumRes3, [:AY,Symbol("True Outstanding Amount"),Symbol("Estimated Reserves"),Symbol("CL Estimated Reserves")])
        comparisonByLOB[lob] = deepcopy(sumRes3)
    end

    comparisons = Dict{String,DataFrame}()
    # byVars=["cc","inj_part","age","LoB"]
    byVars = ["tenGroupsSortedByUltimate", "twentyGroupsSortedByUltimate","cc","inj_part","age","LoB"]
    for vv in byVars
        # attach explanatory data to the tree estimates
        treeEstimatePerRow[Symbol(vv)] = dataKnownByYE2005[Symbol(vv)]
        # create comparison table
        comparison = DataFrames.aggregate(fullData[vcat(Symbol(vv), :PayCum11, :paidToDate)], Symbol(vv), sum)
        # rename table
        names!(comparison, vcat(Symbol(vv), :PayCum11, :paidToDate))
        # calculate reserves
        comparison[Symbol("True Outstanding Amount")] = comparison[:PayCum11] .- comparison[:paidToDate]
        sort!(comparison, Symbol(vv))
        # calculate estimated reserves (tree model)
        est = DataFrames.aggregate(treeEstimatePerRow[vcat(Symbol(vv), :ultimate, :reserves, :paidToDate)], Symbol(vv), sum)
        sort!(est, Symbol(vv))
        # rename table
        names!(est, vcat(Symbol(vv), :ultimate, :reserves, :paidToDate))
        # check
        @assert all(isapprox.(est[:paidToDate] .- comparison[:paidToDate], 0))
        # attach estimated data (this should be a join (but as long as all values occurr in the data, we are fine))
        comparison[Symbol("Estimated Reserves")] = est[:reserves]
        # reorder columns
        comparison = comparison[[1,4,5,3,2]]
        # sort the table
        if !isa(comparison[1,Symbol(vv)], Number)
            comparison[:sort] = Meta.parse.(comparison[Symbol(vv)])
        else
            comparison[:sort] = deepcopy(comparison[Symbol(vv)])
        end
        sort!(comparison, :sort)
        # add variable name column
        comparison[:variable] = vv
        comparison[:value] = string.(comparison[Symbol(vv)])
        delete!(comparison, Symbol(vv))
        # reorder columns
        comparison = comparison[[7,1,2,3,4,5,6]]

        comparisons[vv] = deepcopy(comparison)
        if writeResultToTemp
            if isdir("C:\\temp\\")
                CSV.write(string("C:\\temp\\", vv, ".csv"), comparison)
            end
        end
        println(vv, " : ", ceil(Int, sum(abs.(comparison[Symbol("True Outstanding Amount")] .- comparison[Symbol("Estimated Reserves")])) / 1000))
    end

    return comparisonByAY2, comparisonByLOB, absErrorsByAY, comparisons
end

function runModels!(dataKnownByYE2005, dtmKnownByYE2005, modelsWeightsPerLDF, treeResults, treeResultsAgg, selected_explanatory_vars, categoricalVars, folderForOutput, clAllLOBs, paidToDatePerRow, settOrig::ModelSettings)
    kk = 0
    Random.seed!(1240)
    ayears = sort(unique(dataKnownByYE2005[:AY]))
    maxAY = maximum(ayears)

    # array of LDFs
    LDFArray = zeros(size(dataKnownByYE2005, 1), length(ayears))
    LDFArray[:,end] .= 1.0 # tail factor 1.0
    LDFArraySmoothed = deepcopy(LDFArray)
    LDFArrayUnSmoothed = deepcopy(LDFArray)
    some_stars = "******************************************************************************"

    for minWeightPerLDF in modelsWeightsPerLDF
        fill!(LDFArray, 0.0)
        LDFArray[:,end] .= 1.0 # tail factor 1.0
        kk += 1
        for ldfYear = 1:length(ayears) - 1
            selectedWeight = minWeightPerLDF[ldfYear]
            println("")
            println(some_stars)
            println(some_stars)
            println("Modelling LDF Year $(ldfYear). Selected minimum weight per leaf is $(selectedWeight)")
            println(some_stars)
            println(some_stars)

            resultingFiles, resM = runSingleModel(dataKnownByYE2005, dtmKnownByYE2005, selectedWeight, ldfYear, maxAY, selected_explanatory_vars, categoricalVars, folderForOutput, settOrig)
            # apply tree to the known claims
            if sett.model_type == "build_tree"
                fittedValues, leafNrs = DecisionTrees.predict(resM, dtmKnownByYE2005.features)
            else
                dfpred = DecisionTrees.predict(resM, dtmKnownByYE2005.features)
                # @show ldfYear
                fittedValues = dfpred[:RawEstimate]
                LDFArraySmoothed[:,ldfYear] .= dfpred[:SmoothedEstimate]
                LDFArrayUnSmoothed[:,ldfYear] .= dfpred[:UnsmoothedEstimate]
            end
            # save estimated ldf per observation
            LDFArray[:,ldfYear] .= copy(fittedValues)
            if (!(sett.model_type == "build_tree")) && (length(ayears) - 1 == ldfYear)
                @warn("BK: Over writing ldf factors for year=$(ldfYear)") # there seems to be an issue for the last year in the model when no split is found
                LDFArray[:,ldfYear] .= median(LDFArrayUnSmoothed[:,ldfYear])
            end
        end

    # consider estimated ultimates per claim
        estPerRow, estAgg = calculateEstimates(dataKnownByYE2005, LDFArray, clAllLOBs, paidToDatePerRow)
        treeResults[kk] = deepcopy(estPerRow)
        treeResultsAgg[kk] = deepcopy(estAgg)
    # for boosting also store the other estimates
        if !(sett.model_type == "build_tree")
            estPerRow, estAgg = calculateEstimates(dataKnownByYE2005, LDFArraySmoothed, clAllLOBs, paidToDatePerRow)
            treeResults[kk + 0.2] = deepcopy(estPerRow)
            treeResultsAgg[kk + 0.2] = deepcopy(estAgg)
            estPerRow, estAgg = calculateEstimates(dataKnownByYE2005, LDFArrayUnSmoothed, clAllLOBs, paidToDatePerRow)
            treeResults[kk + 0.4] = deepcopy(estPerRow)
            treeResultsAgg[kk + 0.4] = deepcopy(estAgg)
        end

    end

    return LDFArray
end

function runSingleModel(dataKnownByYE2005, dtmKnownByYE2005, selectedWeight, ldfYear, maxAY, selected_explanatory_vars, categoricalVars, folderForOutput, settOrig::ModelSettings)
    # a note on the data
    # dataKnownByYE2005 is a DataFrame (it contains more data such as all PayCum columns and the AY)
    # dtmKnownByYE2005 is a DTMTable object (it has a different structre than the DataFrame)
    cumulativePaymentCol = Symbol(string("PayCum", lpad(ldfYear, 2, "0")))
    cumulativePaymentColPrev = Symbol(string("PayCum", lpad(ldfYear - 1, 2, "0")))

    # subset data, conditions are:
    # 1) AY < 2005-ldfYear
    # 2) positive denominator
    # 3) cumPay is NOT zero for both ldfYear and ldfYear-1 (these rows do not add any value and can be discarded)
    subsetIdx = (dataKnownByYE2005[:AY] .<= (maxAY - ldfYear)) .& (dataKnownByYE2005[cumulativePaymentCol] .> 0) .& (.!((dataKnownByYE2005[cumulativePaymentCol] .== 0) .& (dataKnownByYE2005[cumulativePaymentColPrev] .== 0)))
    # subsetIdx is a boolean index, we need to convert it to an integer index for the next step
    subsetIdxInteger = findall(subsetIdx)
    # this step also creates a copy of the data (which is intended)
    dtmSubset = dtmKnownByYE2005[subsetIdxInteger]
    # redefine numerator and denominator (we are considering different PayCum columns of the data for each development year)
    dtmSubset.numerator = dataKnownByYE2005[cumulativePaymentColPrev][subsetIdxInteger]
    dtmSubset.denominator = dataKnownByYE2005[cumulativePaymentCol][subsetIdxInteger]

    sett = deepcopy(settOrig)
    critSTR = string(sett.crit)[15:22]
    updateSettingsMod!(sett, minWeight=selectedWeight)
    resultingFiles, resM = dtm(dtmSubset, sett, file=joinpath(folderForOutput, string(sett.model_type[1:7], "_LDF_Year_", ldfYear, "crt_", critSTR, "_minw_", sett.minWeight, ".txt")))

    return resultingFiles, resM
end

function addTotal!(d::DataFrame)
    sm = [sum(d[i]) for i in 1:size(d, 2)]
    totRow = d[1,:]
    nm = names(totRow)
    for j = 1:size(d, 2)
        totRow[j] = 0.0
        if in(nm, [:LDFS,:factorsToUltimate])
        else
            totRow[j] = sm[j]
        end
    end
    totRow[:LDFS] = 1.0
    totRow[:factorsToUltimate] = totRow[:ultimate][end] / totRow[:paidToDate][end]
    append!(d, totRow)
    return nothing
end

function CL_est_per_variable(fullData, cldf::DataFrame;writeResultToTemp=true) # (dataKnownByYE2005,CLLDFperRow,CLUltimatePerRow,paidToDatePerRow)
    # @assert size(dataKnownByYE2005,1)==length(CLUltimatePerRow)==length(CLLDFperRow)==length(paidToDatePerRow)
    comparisons = Dict{String,DataFrame}()
    # byVars=["cc","inj_part","age","LoB"]
    byVars = ["cc","inj_part","age","LoB","tenGroupsSortedByUltimate", "twentyGroupsSortedByUltimate"]
    for vv in byVars
        # attach explanatory data to the data
        cldf[Symbol(vv)] = dataKnownByYE2005[Symbol(vv)]
        # create comparison table
        comparison = DataFrames.aggregate(fullData[vcat(Symbol(vv), :PayCum11, :paidToDate)], Symbol(vv), sum)
        # rename table
        names!(comparison, vcat(Symbol(vv), :PayCum11, :paidToDate))
        # calculate reserves
        comparison[Symbol("True Outstanding Amount")] = comparison[:PayCum11] .- comparison[:paidToDate]
        sort!(comparison, Symbol(vv))
        # calculate estimated reserves
        est = DataFrames.aggregate(cldf[vcat(Symbol(vv), :ultimate, :reserves, :paidToDate)], Symbol(vv), sum)
        sort!(est, Symbol(vv))
        # rename table
        names!(est, vcat(Symbol(vv), :ultimate, :reserves, :paidToDate))
        # check
        @assert all(isapprox.(est[:paidToDate] .- comparison[:paidToDate], 0))
        # attach estimated data (this should be a join (but as long as all values occurr in the data, we are fine))
        comparison[Symbol("Estimated Reserves")] = est[:reserves]
        # reorder columns
        comparison = comparison[[1,4,5,3,2]]
        # sort the table
        if !isa(comparison[1,Symbol(vv)], Number)
            comparison[:sort] = Meta.parse.(comparison[Symbol(vv)])
        else
            comparison[:sort] = deepcopy(comparison[Symbol(vv)])
        end
        sort!(comparison, :sort)
        # add variable name column
        comparison[:variable] = vv
        comparison[:value] = string.(comparison[Symbol(vv)])
        delete!(comparison, Symbol(vv))
        # reorder columns
        comparison = comparison[[7,1,2,3,4,5,6]]

        comparisons[vv] = deepcopy(comparison)
        if writeResultToTemp
            if isdir("C:\\temp\\")
                CSV.write(string("C:\\temp\\", vv, ".csv"), comparison)
            end
        end
        println(vv, " : ", ceil(Int, sum(abs.(comparison[Symbol("True Outstanding Amount")] .- comparison[Symbol("Estimated Reserves")])) / 1000))
    end

    return comparisons
end

function write_results(treeResults, treeResultsAgg, folderForOutput, trueReservesPerLOB, clPerLOB)
    # if true
    @info("Deriving statistics...")
    # compare results (CL versus Tree versus Truth)
    modelIndices = sort(collect(keys(treeResultsAgg)))
    modelStatistics = Dict{Float64,ModelStats}()
    for thiskey in modelIndices
        @show thiskey
        resultTuple = customSummary(treeResultsAgg[thiskey], treeResults[thiskey], trueReservesPerLOB, clPerLOB)
        obj = ModelStats(resultTuple[1], resultTuple[2], resultTuple[3], resultTuple[4])
        deltaTotalPerLOB = map(x->1 / 1e6 .* (obj.comparisonByLOB[x][end,2] .- obj.comparisonByLOB[x][end,3]), string.(collect(1:4)))
        @show deltaTotalPerLOB
        modelStatistics[thiskey] = obj

        @info("Writing data to disk...")
        # write tables to csv file
        summaryByVar = vcat(collect(values(modelStatistics[thiskey].comparisons)))
        fi = joinpath(folderForOutput, string("results", thiskey, ".csv"))
        CSV.write(fi, summaryByVar)
    end

    for kk in keys(modelStatistics)
        @show kk
        obj = modelStatistics[kk]
        trueTotal = map(x->1 / 1e6 .* (obj.comparisonByLOB[x][end,2]), string.(collect(1:4)))
        treeTotal = map(x->1 / 1e6 .* (obj.comparisonByLOB[x][end,3]), string.(collect(1:4)))
        clTotal = map(x->1 / 1e6 .* (obj.comparisonByLOB[x][end,4]), string.(collect(1:4)))
        deltaTotalPerLOB = map(x->1 / 1e6 .* (obj.comparisonByLOB[x][end,2] .- obj.comparisonByLOB[x][end,3]), string.(collect(1:4)))
        @show deltaTotalPerLOB, sum(deltaTotalPerLOB)
        @show clTotal, sum(clTotal)
        @show treeTotal, sum(treeTotal)
        @show trueTotal, sum(trueTotal)
    end
    return modelStatistics
end

function test_on_subset(dataKnownByYE2005, ldfYear, folderForOutput, settOrig::ModelSettings, trnsize::Integer, selectedWeight::Integer, cvsampler)
    cumulativePaymentCol = Symbol(string("PayCum", lpad(ldfYear, 2, "0")))
    cumulativePaymentColPrev = Symbol(string("PayCum", lpad(ldfYear - 1, 2, "0")))

    # subset data, conditions are:
    # 1) AY < 2005-ldfYear
    # 2) positive denominator
    # 3) cumPay is NOT zero for both ldfYear and ldfYear-1 (these rows do not add any value and can be discarded)
    subsetIdx = (dataKnownByYE2005[:AY] .<= (maxAY - ldfYear)) .& (dataKnownByYE2005[cumulativePaymentCol] .> 0) .& (.! ((dataKnownByYE2005[cumulativePaymentCol] .== 0) .& (dataKnownByYE2005[cumulativePaymentColPrev] .== 0)))
    # subsetIdx is a boolean index, we need to convert it to an integer index for the next step
    subsetIdxInteger = findall(subsetIdx)

    @warn("Randomly sampling from data")
    # consider random subset
    rnd = rand(length(subsetIdxInteger))
    srtTmp = sortperm(rnd)
    subsetIdxInteger = subsetIdxInteger[srtTmp[1:min(trnsize, length(subsetIdxInteger))]]

    # this step also creates a copy of the data (which is intended)
    dtmSubset = dtmKnownByYE2005[subsetIdxInteger]
    # redefine numerator and denominator (we are considering different PayCum columns of the data for each development year)
    dtmSubset.numerator = dataKnownByYE2005[cumulativePaymentColPrev][subsetIdxInteger]
    dtmSubset.denominator = dataKnownByYE2005[cumulativePaymentCol][subsetIdxInteger]

    sett = deepcopy(settOrig)
    critSTR = string(sett.crit)[15:22]
    updateSettingsMod!(sett, minWeight=selectedWeight)
    if cvsampler.folds == 0
        @show 1
        resultingFiles, resM = dtm(dtmSubset, sett, file=joinpath(folderForOutput, string(sett.model_type[1:7], "_LDF_Year_", ldfYear, "crt_", critSTR, "_minw_", sett.minWeight, ".txt")))
        return resM
    else
        @show 2
        statsdf, settsdf, allmodel = dtm(dtmSubset, sett, cvsampler, file=joinpath(folderForOutput, string(sett.model_type[1:7], "_LDF_Year_", ldfYear, "crt_", critSTR, "_minw_", sett.minWeight, ".txt")))
        return  statsdf, settsdf, allmodel
    end
end

"""
for a sorted vector vsorted this function returns the percentage of the aggregate amount which is represented by the smallest "claim_count_pct" percent of elements (e.g. claims)
"""
function aggLossPctSortedInput(claim_count_pct, vsorted, vsum)
    @assert 1 >= claim_count_pct >= 0
    # @assert issorted(vsorted)
    claim_count = floor(Int, length(vsorted) * claim_count_pct)
    if claim_count > 0
        p = sum(view(vsorted, 1:claim_count)) / vsum # /sum(vsorted)
        # @show claim_count,p
    else
        p = zero(eltype(vsorted))
    end
    return p
end

function count_vs_agg_amount(vsorted, nsteps::Int)
    @assert issorted(vsorted)
    @assert nsteps > 0
    res = zeros(eltype(vsorted), nsteps, 4)
    vsize = length(vsorted)
    vsum = sum(vsorted)
    startPrev = vsize
    for i = nsteps:-1:1
        j = nsteps - i + 1
        start = floor(Int, (i / nsteps) * vsize)
        @inbounds res[j,3] = sum(view(vsorted, start:startPrev))
        startPrev = start
        @inbounds res[i,1] = i / nsteps
        if j > 1
            res[j,3] += res[j - 1,3]
        end
    end
    res[:,4] = res[:,3] ./ vsum
    res[:,2] = res[:,1] .* vsize
    # output format: row count in percentage, row count, aggregate sum, aggregate sum in percentage
    return res
end

# count_vs_agg_amount(vsorted,10)
#=
tbl=count_vs_agg_amount(vsorted,10);
@btime count_vs_agg_amount($vsorted,1000); =#

"""
for a sorted vector vsorted this function returns the percentage of the aggregate amount which is represented by the smallest "claim_count_pct" percent of elements (e.g. claims)
"""
function aggLossPctSortedInput(claim_count_pct, vsorted, vsum)
    @assert 1 >= claim_count_pct >= 0
    # @assert issorted(vsorted)
    claim_count = floor(Int, length(vsorted) * claim_count_pct)
    if claim_count > 0
        p = sum(view(vsorted, 1:claim_count)) / vsum # /sum(vsorted)
        # @show claim_count,p
    else
        p = zero(eltype(vsorted))
    end
    return p
end


text:
using AlgebraicDynamics
using AlgebraicDynamics.UWDDynam
using AlgebraicDynamics.DWDDynam

using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra

using DynamicalSystems
using OrdinaryDiffEq

using Test
const UWD = UndirectedWiringDiagram

approx_equal(u, v) = abs(maximum(u - v)) < .1

# Resource sharer tests
dx(x) = [1 - x[1]^2, 2*x[1]-x[2]]
dy(y) = [1 - y[1]^2]

r = ContinuousResourceSharer{Real}(2, (u,p,t) -> dx(u))

u0 = [-1.0, -2.0]
tspan = (0.0, 100.0)
prob = ODEProblem(r, u0, tspan)
@test prob isa ODEProblem
sol = solve(prob, Tsit5())
@test reduce(*, map(1:length(sol.t)) do i
  sol[i] == u0
end)

u0 = [.1, 10.0]
prob = ODEProblem(r, u0, tspan)
sol = solve(prob, Tsit5())
@test approx_equal(last(sol), [1.0, 2.0])

h = 0.2
dr = euler_approx(r, h)
dds = DiscreteProblem(dr, u0, tspan, nothing) # nothing because no parameters
t = solve(dds, FunctionMap())
@test approx_equal(last(t), [1.0, 2.0])

dr = DiscreteResourceSharer{Real}(1, (u,p,t) -> dy(u))
u0 = [1.0]
dds = DiscreteProblem(dr, u0, tspan, nothing)
t = solve(dds, FunctionMap())
@test reduce(*, map(1:length(sol.t)) do i
  t[i][1] == i%2
end)

# Machine tests
uf(u, x, p, t) = [x[1] - u[1]*x[2]]
rf(u,p,t) = u
u0 = [1.0]
m = ContinuousMachine{Float64}(2,1, 1, uf, rf)
xs = [t -> 1, t -> 1]
prob = ODEProblem(m, u0, xs, tspan)
@test prob isa ODEProblem

sol = solve(prob, Tsit5())
@test reduce(*, map(1:length(sol.t)) do i
  sol[i] == u0
end)

xs = t->t
prob = ODEProblem(m, u0, xs, tspan)
sol = solve(prob, Tsit5())
@test reduce(*, map(1:length(sol.t)) do i
    sol[i] == u0
  end)

xs = [2, 1]
prob = ODEProblem(m, u0, xs, tspan)
sol = solve(prob, Tsit5())
@test approx_equal(last(sol), [2.0])

dm = euler_approx(m,h)
prob = DiscreteProblem(dm, u0, xs, (0, 100.0), nothing)
s = solve(prob, FunctionMap())
@test approx_equal(last(s), [2.0])

# machines - oapply
uf(u, x, p, t) = [x[1] - u[1], 0.0]
rf(u,p,t) = u
mf = ContinuousMachine{Float64}(2,2,2, uf, rf)

d_id = singleton_diagram(Box(:f, [:A, :A], [:A, :A]))
m_id = oapply(d_id, [mf])
xs = [0.0, 0.0]
u0 = [10.0, 2.0]
prob1 = ODEProblem(mf, u0, xs, tspan)
prob2 = ODEProblem(m_id, u0, xs, tspan)

sol1 = solve(prob1, Tsit5(); dtmax = 1)
sol2 = solve(prob2, Tsit5(); dtmax = 1)
@test approx_equal(last(sol1), [0.0, 2.0])
@test approx_equal(last(sol2), [0.0, 2.0])

# Lokta Volterra models

# as resource sharers
params = [0.3, 0.015, 0.015, 0.7]
dotr(u,p,t) = p[1]*u
dotrf(u,p,t) = [-p[2]*u[1]*u[2], p[3]*u[1]*u[2]]
dotf(u,p,t) = -p[4]*u

r = ContinuousResourceSharer{Real}(1, dotr)
rf_pred = ContinuousResourceSharer{Real}(2, dotrf)
f = ContinuousResourceSharer{Real}(1, dotf)

rf_pattern = UWD(0)
add_box!(rf_pattern, 1); add_box!(rf_pattern, 2); add_box!(rf_pattern, 1)
add_junctions!(rf_pattern, 2)
set_junction!(rf_pattern, [1,1,2,2])

lv = oapply(rf_pattern, [r, rf_pred, f])

params = [0.3, 0.015, 0.015, 0.7]
u0 = [10.0, 100.0]
prob = ODEProblem(lv, u0, tspan, params)
sol = solve(prob, Tsit5())

@test reduce(*, map(1:length(sol.t)) do i
  ( 0 < sol[i][1] < 200.0 ) && (0 < sol[i][2] < 150.0)
end)

h = 0.01
lv_discrete = oapply(rf_pattern, euler_approx([r, rf_pred, f], h))
dds = DiscreteProblem(lv_discrete, u0, (0,100.0), params)
sol2 = solve(dds, FunctionMap(); dt = h)
@test reduce(*, map(1:length(sol.t)) do i
  ( 0 < sol[i][1] < 200.0 ) && (0 < sol[i][2] < 150.0)
end)

# as machines
dotr(u, x, p, t) = [p[1]*u[1] - p[2]*u[1]*x[1]]
dotf(u, x, p, t) = [p[3]*u[1]*x[1] - p[4]*u[1]]

rmachine = ContinuousMachine{Real}(1,1,1, dotr, (r,p,t) -> r)
fmachine = ContinuousMachine{Real}(1,1,1, dotf, (f,p,t) -> f)

rf_pattern = WiringDiagram([],[])
boxr = add_box!(rf_pattern, Box(nothing, [nothing], [nothing]))
boxf = add_box!(rf_pattern, Box(nothing, [nothing], [nothing]))
add_wires!(rf_pattern, Pair[
  (boxr, 1) => (boxf, 1),
  (boxf, 1) => (boxr, 1)
])

u0 = [10.0, 100.0]
params = [0.3, 0.015, 0.015, 0.7]

rf_machine = oapply(rf_pattern, [rmachine, fmachine])
prob = ODEProblem(rf_machine, u0, tspan, params)
sol = solve(prob, Tsit5(); dtmax = .1)
@test reduce(*, map(1:length(sol.t)) do i
  ( 0 < sol[i][1] < 200.0 ) && (0 < sol[i][2] < 150.0)
end)

h = 0.005
lv_discrete = oapply(rf_pattern, euler_approx([rmachine, fmachine], h))
dds = DiscreteProblem(lv_discrete, u0, (0.0, 10.0),params)
sol = solve(dds, FunctionMap(); dt = h)
@test reduce(*, map(1:length(sol.t)) do i
  ( 0 < sol[i][1] < 200.0 ) && (0 < sol[i][2] < 150.0)
end)

# ocean
α, β, γ, δ, β′, γ′, δ′ = 0.3, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35
params = [α, β, γ, δ, β′, γ′, δ′]

dotfish(f, x, p, t) = [p[1]*f[1] - p[2]*x[1]*f[1]]
dotFISH(F, x, p, t) = [p[3]*x[1]*F[1] - p[4]*F[1] - p[5]*x[2]*F[1]]
dotsharks(s, x, p, t) = [-p[7]*s[1] + p[6]*s[1]*x[1]]

fish   = ContinuousMachine{Real}(1,1,1, dotfish,   (f,p,t) ->f)
FISH   = ContinuousMachine{Real}(2,1,2, dotFISH,   (F,p,t)->[F[1], F[1]])
sharks = ContinuousMachine{Real}(1,1,1, dotsharks, (s,p,t)->s)

ocean_pat = WiringDiagram([], [])
boxf = add_box!(ocean_pat, Box(nothing, [nothing], [nothing]))
boxF = add_box!(ocean_pat, Box(nothing, [nothing, nothing], [nothing, nothing]))
boxs = add_box!(ocean_pat, Box([nothing], [nothing]))
add_wires!(ocean_pat, Pair[
  (boxf, 1) => (boxF, 1),
  (boxs, 1) => (boxF, 2),
  (boxF, 1) => (boxf, 1),
  (boxF, 2) => (boxs, 1)
])

u0 = [100.0, 10.0, 2.0]
ocean = oapply(ocean_pat, [fish, FISH, sharks])
prob = ODEProblem(ocean, u0, tspan, params)
sol = solve(prob, Tsit5(); dtmax = 0.1)
@test reduce(*, map(1:length(sol.t)) do i
  0 < sol[i][1] < 200.0 && (0 < sol[i][2] < 75) && (0 < sol[i][3] < 20)
end)

h = 0.01
dds = DiscreteProblem(euler_approx(ocean, h), u0, (0.0, 10.0), params)
sol = solve(dds, FunctionMap(), dt = h)
@test reduce(*, map(1:length(sol.t)) do i
  0 < sol[i][1] < 200.0 && (0 < sol[i][2] < 75) && (0 < sol[i][3] < 20)
end)

ext: jl


text:
export snn_controller

# A fully-connected snn baseline controller with
# + LIF neurons
# + neurogenesis
# + synaptogenesis (on neurogenesis)
# + RM-STDP

function snn_gen_cell_division(cfg::Dict)
    function snn_cell_division(cell_params::Array{Float64})
        ((cell_params[4] > cfg["divide_age"]) &&
         (cell_params[5] > cfg["divide_act"]))
     end
     snn_cell_division
end

function snn_new_cell_params(cell_params::Array{Float64})
    vcat(copy(cell_params[1:3]) + 0.01 * randn(3), 1.0, 0.0)
end

function snn_gen_cell_death(cfg::Dict)
    function snn_cell_death(cell_params::Array{Float64})
        cell_params[5] < cfg["cell_apop"]
    end
    snn_cell_death
end

function snn_gen_cell_state_update(cfg::Dict)
    function snn_cell_state_update(inputs::Array{Float64})
        # inputs: cell_params, cell_state, inputs
        cell_state = inputs[cfg["n_cell_params"] .+ (1:cfg["n_cell_state"])]
        einputs = @view inputs[(cfg["n_cell_params"]+cfg["n_cell_state"]+1):end]
        spike = cell_state[1] > (cfg["v_thresh"] + cell_state[2])
        outputs = zeros(cfg["n_channels"])
        if spike
            cell_state[1] = cfg["v_reset"]
            cell_state[2] = min(1.0, cell_state[2] + cfg["v_thresh_inc"])
            outputs[1] = cfg["spike_exc"]
            outputs[2] = -cfg["spike_inhib"]
        else
            cell_state[1] = min(1.0, max(-1.0,
                cfg["v_leak"] * cell_state[1] + einputs[1] + einputs[2]))
            cell_state[2] = cfg["v_thresh_dec"] * cell_state[2]
        end
        vcat(cell_state, outputs)
    end
    snn_cell_state_update
end

function snn_gen_cell_param_update(cfg::Dict)
    function snn_cell_param_update(inputs::Array{Float64})
        cell_params = copy(inputs[1:cfg["n_cell_params"]])
        cell_state = @view inputs[(cfg["n_cell_params"]+1):end]
        cell_params[4] *= 0.999
        cell_params[5] = cell_state[2]
        cell_params
    end
    snn_cell_param_update
end

function snn_connect(inputs::Array{Float64})
    true
end

function snn_gen_new_conn_params(cfg::Dict)
    function snn_new_conn_params(inputs::Array{Float64})
        # inputs: c1 params, c2 params
        c1 = @view inputs[1:cfg["n_cell_params"]]
        c2 = @view inputs[(cfg["n_cell_params"]+1):end]
        params = zeros(cfg["n_conn_params"])
        if c2[3] > c1[3]
            params[1] = randn()*cfg["weight_std"] + cfg["weight_mean"]
        end
        params[2] = exp(-euclidean(c1[1:3], c2[1:3]))
        max.(-1.0, min.(1.0, params))
    end
    snn_new_conn_params
end

function snn_disconnect(inputs::Array{Float64})
    false
end

function snn_gen_conn_state_update(cfg::Dict)
    function snn_conn_state_update(inputs::Array{Float64})
        # inputs: conn params, state, inputs
        params = @view inputs[1:cfg["n_conn_params"]]
        state = inputs[cfg["n_conn_params"] .+ (1:cfg["n_conn_state"])]
        einputs = @view inputs[(cfg["n_conn_params"]+cfg["n_conn_state"]+1):end]
        state .*= cfg["trace_dec"]
        nstate = min.(1.0, max.(-1.0, einputs + state))
        outputs = params .* einputs
        vcat(nstate, outputs)
    end
    snn_conn_state_update
end

function snn_gen_conn_param_update(cfg::Dict)
    function snn_conn_param_update(inputs::Array{Float64})
        # inputs: c1 params, c2 params, conn params, state
        params = inputs[2*cfg["n_cell_params"] .+ (1:cfg["n_conn_params"])]
        state = @view inputs[((2*cfg["n_cell_params"] + cfg["n_conn_params"])
                              .+ (1:cfg["n_conn_state"]))]
        dw = state[2] * cfg["stdp_lr"] * (
            (state[1] - cfg["pre_target"]) * (1.0 - params[1]))
        params[1] += dw
        max.(-1.0, min.(1.0, params))
    end
    snn_conn_param_update
end

function snn_controller(cfg::Dict)
    Controller(snn_gen_cell_division(cfg),
               snn_new_cell_params,
               snn_gen_cell_death(cfg),
               snn_gen_cell_state_update(cfg),
               snn_gen_cell_param_update(cfg),
               snn_connect,
               snn_gen_new_conn_params(cfg),
               snn_disconnect,
               snn_gen_conn_state_update(cfg),
               snn_gen_conn_param_update(cfg))
end


text:
using DataFrames
using Gadfly
using Cairo
using Fontconfig

# Input parameters
path = "/home/cns/Dropbox/collapse_code/"
abundanceFile = "structure_evolution.dat"
plotname = "profiles.pdf"
maxtime = 1.0e7
mintime = 1.0e3

plottheme = Gadfly.Theme(
                        default_point_size=1.5pt,
                        line_width=2pt,
                        minor_label_font_size=10pt,
                        minor_label_color=colorant"black",
                        major_label_font_size=12pt,
                        major_label_color=colorant"black",
                        key_label_color=colorant"black",
                        key_title_color=colorant"black",
                        key_title_font_size=12pt,
                        key_label_font_size=10pt,
                        grid_color=colorant"black",
                        highlight_width=0pt)

println("Reading input file")
structure = readtable(
                  path*abundanceFile,
                  names = [:time, :Av, :n, :tgas, :tgrain],
                  eltypes = [Float64, Float64, Float64, Float64, Float64],
                  skipstart = 1,
                  separator = ' '
                 )

structure[:Av]     = 10.^structure[:Av]
structure[:n]      = 10.^structure[:n]
structure[:tgas]   = 10.^structure[:tgas]
structure[:tgrain] = 10.^structure[:tgrain]

structure = structure[mintime.<structure[:time].<maxtime,:]

println("Now plotting data")

# Plot time vs. Av
p1 = plot(
           structure,
           x="time",
           y="Av",
           Coord.cartesian(aspect_ratio=1),
           Geom.line,
           Guide.xlabel("Time (yr)"),
           Guide.ylabel("Av"),
           Scale.x_log10(maxvalue=maxtime,minvalue=mintime),
           Scale.y_log10,
           )
push!(p1,plottheme)

# Plot time vs. n
p2 = plot(
           structure,
           x="time",
           y="n",
           Coord.cartesian(aspect_ratio=1),
           Geom.line,
           Guide.xlabel("Time (yr)"),
           Guide.ylabel("n (cm-3)"),
           Scale.x_log10(maxvalue=maxtime,minvalue=mintime),
           Scale.y_log10,
           )
push!(p2,plottheme)

# Plot time vs. Tgas
p3 = plot(
           structure,
           x="time",
           y="tgas",
           Coord.cartesian(aspect_ratio=1),
           Geom.line,
           Guide.xlabel("Time (yr)"),
           Guide.ylabel("Tgas (K)"),
           Scale.x_log10(maxvalue=maxtime,minvalue=mintime),
           Scale.y_log10,
           )
push!(p3,plottheme)

# Plot time vs. Tgrain
p4 = plot(
           structure,
           x="time",
           y="tgrain",
           Coord.cartesian(aspect_ratio=1),
           Geom.line,
           Guide.xlabel("Time (yr)"),
           Guide.ylabel("Tgrain (K)"),
           Scale.x_log10(maxvalue=maxtime,minvalue=mintime),
           Scale.y_log10,
           )
push!(p4,plottheme)

println("Now exporting plot")
draw(PDF(path*"av.pdf", 6inch, 5inch), p1)
draw(PDF(path*"density.pdf", 6inch, 5inch), p2)
draw(PDF(path*"gastemp.pdf", 6inch, 5inch), p3)
draw(PDF(path*"graintemp.pdf", 6inch, 5inch), p4)
draw(PDF(path*"profiles.pdf", 6inch, 5inch), vstack(hstack(p1,p2),hstack(p3,p4)))
println("Now ending script")


text:
import PredictMDExtra
import PredictMD

function _predictmdfull_import_all(env::AbstractDict = ENV)::Bool
    raw_env_value = strip(
        lowercase(
            strip(
                get(env, "PREDICTMDFULL_IMPORT_ALL", "")
                )
            )
        )
    return length(raw_env_value) == 0 || raw_env_value == "true"
end

_import_all_on_init() = _import_all_on_init(Main)

function _import_all_on_init(m::Module)::Nothing
    if _predictmdfull_import_all()
        import_all(m)
    end
end

_import_all() = _import_all(Main)
import_all() = import_all(Main)

function _import_all(m::Module)::Nothing
    package_list::Vector{String} = sort(unique(strip.(_package_list())))
    for p in package_list
        Base.eval(
            m,
            Base.Meta.parse(
                string(
                    "try ",
                    "import $(string(p)); ",
                    "catch e2 ",
                    "@debug(\"ignoring exception: \", e2,); ",
                    "try ",
                    "import Pkg; ",
                    "Pkg.add(\"$(string(p))\"); ",
                    "catch e3 ",
                    "@debug(\"ignoring exception: \", e3,); ",
                    "end ",
                    "end ",
                    ),
                ),
            )
        Base.eval(
            m,
            Base.Meta.parse(
                string(
                    "try ",
                    "import $(string(p)); ",
                    "@debug(\"imported $(string(p))\"); ",
                    "catch e4 ",
                    "@debug(\"ignoring exception: \", e4,); ",
                    "end ",
                    ),
                ),
            )
    end
    return nothing
end

function import_all(m::Module)::Nothing
    _import_all(m)
    PredictMDExtra.import_all(m)
    PredictMD.import_all(m)
    return nothing
end


text:
function f(x::Int64)
  x ^ 2
end

n = 10
println(typeof(n)) # -> Int64
println(f(n))

f2(n::Integer) = n ^ 3

# println(f2(10.1)) # Error
println(f2(10)) # Error

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
"""
# Lagging
```julia
lag(ts::TS, lag_value::Int = 1)
```

Lag the `ts` object by the specified `lag_value`. The rows corresponding
to lagged values will be rendered as `missing`. Negative values of lag are
also accepted (see `TSx.lead`).

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random, Statistics;

julia> random(x) = rand(MersenneTwister(123), x);

julia> dates = collect(Date(2017,1,1):Day(1):Date(2017,1,10));

julia> ts = TS(random(length(dates)), dates);
julia> show(ts)
(10 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
───────────────────────
 2017-01-01  0.768448
 2017-01-02  0.940515
 2017-01-03  0.673959
 2017-01-04  0.395453
 2017-01-05  0.313244
 2017-01-06  0.662555
 2017-01-07  0.586022
 2017-01-08  0.0521332
 2017-01-09  0.26864
 2017-01-10  0.108871

julia> lag(ts)
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64?
─────────────────────────────
 2017-01-01  missing
 2017-01-02        0.768448
 2017-01-03        0.940515
 2017-01-04        0.673959
 2017-01-05        0.395453
 2017-01-06        0.313244
 2017-01-07        0.662555
 2017-01-08        0.586022
 2017-01-09        0.0521332
 2017-01-10        0.26864

julia> lag(ts, 2) # lags by 2 values
(10 x 1) TS with Date Index

 Index       x1
 Date        Float64?
─────────────────────────────
 2017-01-01  missing
 2017-01-02  missing
 2017-01-03        0.768448
 2017-01-04        0.940515
 2017-01-05        0.673959
 2017-01-06        0.395453
 2017-01-07        0.313244
 2017-01-08        0.662555
 2017-01-09        0.586022
 2017-01-10        0.0521332

```
"""
function lag(ts::TS, lag_value::Int = 1)
    sdf = DataFrame(ShiftedArrays.lag.(eachcol(ts.coredata[!, Not(:Index)]), lag_value), TSx.names(ts))
    insertcols!(sdf, 1, :Index => ts.coredata[!, :Index])
    TS(sdf, :Index)
end


text:
using Test
using LIKWID
using LinearAlgebra

# julia should NOT have been started with `likwid-perfctr -C ... -g ... -m` for this test file
@assert !Marker.isactive()

# init
@test isnothing(Marker.init())

N = 100_000_000
a = 3.141f0
z = zeros(Float32, N)
x = rand(Float32, N)
y = rand(Float32, N)

saxpy!(z,a,x,y) = z .= a .* x .+ y

# regular workflow
saxpy!(z,a,x,y)
@test !Marker.startregion("saxpy!")
saxpy!(z,a,x,y)
@test !Marker.stopregion("saxpy!")

# nextgroup
@test isnothing(Marker.nextgroup())

# registerregion
A = rand(100,100)
B = rand(100,100)
@test !Marker.registerregion("mul")
@test !Marker.startregion("mul")
for _ in 1:10
    A * B
end
@test !Marker.stopregion("mul")

# close
@test isnothing(Marker.close())


text:
module HPack

import HttpCommon: Headers

type DecodeError <: Exception
    message::AbstractString
end

# package code goes here
include("table.jl")
include("huffman.jl")
include("encode.jl")
include("decode.jl")

export encode
export decode
export DecodeError

end # module


text:
function round(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    z = ArbFloat{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

function round(x::ArbReal{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    z = ArbReal{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(arb_set_round), Cint, (Ref{ArbReal}, Ref{ArbReal}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

function round(x::ArbComplex{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    z = ArbComplex{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(acb_set_round), Cint, (Ref{ArbComplex}, Ref{ArbComplex}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

for A in (:ArbFloat, :ArbReal, :ArbComplex)
  @eval begin

    round(::Type{T}, x::$A{P}, rounding::RoundingMode=RoundNearest) where {P,T} =
       T(round(x, digits=P, base=rounding))

    function round(x::$A{P}, rounding::RoundingMode=RoundNearest;
                   sigdigits::Integer = 0, digits::Integer = 0, base = 10) where {P}
        sigdigits = max(sigdigits, digits)
        if base == 10
            round(x, digits=digits4bits(sigdigits), base=rounding)
        elseif base ==  2
            round(x, digits=sigdigits, base=rounding)
        else
            throw(DomainError("base ($base) not supported"))
        end
     end
  end

end


text:
# TODO: add some checks for different boundary conditions
# TODO: add checks for analytic wavelets
# ConvFFT constructor tests
@testset "ConvFFT constructors" begin
    @testset "Utils" begin
        explicit = [1 0 0; 0 1 0; 0 0 1; zeros(2,3)]
        @test maximum(abs.(iden_perturbed_gaussian(5,3) - explicit) ) < 1
    end
    @testset "2D constructors" begin
        # normal size
        originalSize = (21,11,1,10)
        x = randn(Float32, originalSize)
        weightMatrix = randn(Float32, (21+10)>>1+1, 11+10, 1)
        weightMatrix = reshape([I zeros(16,5)], (16,21,1))
        padding = (5,5)
        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = Pad(padding))
        @test size(shears.fftPlan)== originalSize .+ (10, 10, 0, 0)
        @test shears.σ == abs
        @test shears.bias == nothing
        @test shears.bc.padBy == (5,5)

        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = Pad(padding), trainable=true)
        @test params(shears).order[1] == shears.weight
        @test length(params(shears).order)==1

        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         boundary = Pad(padding),trainable=false)
        @test isempty(params(shears))

        x = randn(21,11,1,10)
        ∇ = gradient((x)->shears(x)[1,1,1,1,3], x)
        @test minimum(∇[1][:,:,:, [1:2..., 4:10...]] .≈ 0)

        # check that the identity ConvFFT is, in fact, an identity
        weightMatrix = ones(Float32, (21+10)>>1+1, 11+10, 1)
        weightMatrix = cat(weightMatrix, weightMatrix, dims=3)
        padding = (5,5)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Pad(padding))
        x = randn(Float32,21,11,1,10)
        @test shears(x)[:,:,2,:,:] ≈ x

        # check that global multiplication in the Fourier domain is just multiplication
        weightMatrix = 2 .* ones(Float32, (21+10)>>1+1, 11+10, 1)
        padding = (5,5)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Pad(padding))
        x = randn(21,11,1,10)
        @test shears(x)[:,:,1,:,:] ≈ 2 .* x

        # internal methods tests
        x̂ = pad(x, shears.bc.padBy);
        x̂ = shears.fftPlan *ifftshift(x̂,(1,2))
        usedInds1 = shears.bc.padBy[1] .+ (1:size(x, 1))
        usedInds2 = shears.bc.padBy[2] .+ (1:size(x, 2))
        usedInds = (usedInds1, usedInds2)
        nextLayer = FourierFilterFlux.internalConvFFT(x̂, shears.weight, usedInds,
                                                      shears.fftPlan, shears.bias,
                                                      shears.analytic)
        ∇ = gradient((x̂)->FourierFilterFlux.internalConvFFT(x̂,
                                                            shears.weight,
                                                            usedInds,
                                                            shears.fftPlan,
                                                            shears.bias, shears.analytic)[1,1,1,1,1],
                     x̂)
        @test minimum(abs.(diag(∇[1][:, :, 1,1])).≈ 2f0/31/21)

        ax = axes(x̂)[3:end-1]
        ∇ = gradient((x̂) -> FourierFilterFlux.applyWeight(x̂, shears.weight[:,:,1], usedInds,
                                                          shears.fftPlan,
                                                          shears.bias, FourierFilterFlux.NonAnalyticMatching())[1,1,1,1,1], x̂)
        @test minimum(abs.(diag(∇[1][:, :, 1,1])).≈ 2f0/31/21)

        ∇ = gradient((x̂) -> (shears.fftPlan \ (x̂ .* shears.weight[:,:,1]))[1,1,1,1], x̂)
        @test minimum(abs.(diag(∇[1][:, :, 1,1])).≈ 1f0/31*2/21)
        sheared = shears(x)
        @test size(sheared) == (21,11,1,1,10)

        weightMatrix = 2 .* ones(Float32, (21+10)>>1+1, 11+10, 1)
        padding = (5,5)
        originalSize = (21,11,1,10)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Pad(padding))

        # convert to a gpu version
        if CUDA.functional()
            gpuVer = shears |> gpu
            @test typeof(gpuVer.weight) <: CuArray
            @test typeof(gpuVer.fftPlan) <: CUFFT.rCuFFTPlan
            if !(typeof(gpuVer.weight) <: CuArray)
                println("gpuVer.weight is of type $(typeof(gpuVer.weight))")
            end
            if !(typeof(gpuVer.fftPlan) <: CUFFT.rCuFFTPlan)
                println("gpuVer.fftPlan is of type $(typeof(gpuVer.fftPlan))")
            end
            typeof(gpuVer.fftPlan)
        end
        # extra channel dimension
        originalSize = (20,10,16, 1,10)
        shears = ConvFFT(randn(Float32, 16, 20, 3), nothing, originalSize, abs,
                         plan=true, boundary = Pad((5,5)))
        @test size(shears.fftPlan)== originalSize .+ (10, 10, 0, 0, 0)
        @test shears.σ == abs
        @test shears.bias == nothing
        @test shears.bc.padBy == (5,5)
        wSize = originalSize[1:2] .+ (10, 10)
        wSize = (wSize[1]>>1+1, wSize[2], 3)
        @test size(shears.weight) == wSize

        # random initialization
        originalSize = (20, 10, 16, 1, 10)
        shears = ConvFFT(originalSize, 3, abs,
                         plan=true, boundary = Pad((5,5)))
        @test size(shears.fftPlan)== originalSize .+ (10, 10, 0, 0, 0)
        @test shears.σ == abs
        @test size(shears.bias) == (originalSize[3:4]..., 3)
        @test shears.bc.padBy == (5,5)
        wSize = originalSize[1:2] .+ (10, 10)
        wSize = (wSize[1]>>1+1, wSize[2], 3)
        @test size(shears.weight) == wSize
end

@testset "1D constructors" begin
    # normal size
    originalSize = (21,1,10)
    x = randn(Float32, originalSize)
    @testset "basic tests across boundary conditions" begin
        weightMatrix = randn(Float32, (21+10)>>1 + 1, 1)
        padding = (5,)
        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = Pad(padding), trainable=true)
        @test size(shears.fftPlan)== originalSize .+ (10, 0, 0)
        @test shears.σ == abs
        @test shears.bias == nothing
        @test shears.bc.padBy == (5,)
        @test params(shears).order[1] == shears.weight

        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = Pad(padding), trainable=false)
        @test isempty(params(shears))

        x = randn(21,1,10)
        ∇ = gradient((x)->shears(x)[1,1,1,3], x)
        @test minimum(∇[1][:,:, [1:2..., 4:10...]] .≈ 0)

        # Sym test
        weightMatrix = randn(Float32, (21+1), 1)
        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = FourierFilterFlux.Sym())
        @test size(shears.fftPlan)== originalSize .*(2,1,1)
        @test shears.σ == abs
        @test shears.bias == nothing
        @test typeof(shears.bc)<:Sym
        @test params(shears).order[1] == shears.weight
        x = randn(21,1,10)
        ∇ = gradient((x)->shears(x)[1,1,1,3], x)
        @test minimum(∇[1][:,:, [1:2..., 4:10...]] .≈ 0)
        @test minimum(abs.(∇[1][:,1,3]))>0
        weightMatrix = randn(Float32, 21 >>1 + 1, 1)
        shears = ConvFFT(weightMatrix, nothing, originalSize, abs,
                         plan=true, boundary = FourierFilterFlux.Periodic())
        @test size(shears.fftPlan)== originalSize
        @test shears.σ == abs
        @test shears.bias == nothing
        @test typeof(shears.bc)<:FourierFilterFlux.Periodic
        @test params(shears).order[1] == shears.weight
        x = randn(21,1,10)
        ∇ = gradient((x)->shears(x)[1,1,1,3], x)
        @test minimum(∇[1][:,:, [1:2..., 4:10...]] .≈ 0)
        @test minimum(abs.(∇[1][:,1,3]))>0
    end

    # check that the identity ConvFFT is, in fact, an identity
    @testset "Identity tests" begin
        weightMatrix = ones(Float32, (21+10)>>1+1, 1)
        weightMatrix = cat(weightMatrix, weightMatrix, dims=2)
        padding = 5
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Pad(padding))
        x = randn(Float32,21,1,10)
        @test shears(x)[:,1,:,:] ≈ x

        weightMatrix = ones(Float32, 21+1, 1)
        weightMatrix = cat(weightMatrix, weightMatrix, dims=2)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = FourierFilterFlux.Sym())
        x = randn(Float32,21,1,10)
        @test shears(x)[:,1,:,:] ≈ x
        weightMatrix = ones(Float32, 21>>1+1, 1)
        weightMatrix = cat(weightMatrix, weightMatrix, dims=2)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = FourierFilterFlux.Periodic())
        x = randn(Float32,21,1,10)
        @test shears(x)[:,1,:,:] ≈ x
    end

    # check that global multiplication in the Fourier domain is just multiplication
    @testset "times 2" begin

        weightMatrix = 2 .* ones(Float32, (21+10)>>1+1, 1)
        padding = 5
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Pad(padding))
        x = randn(21,1,10)
        @test shears(x)[:,1,:,:] ≈ 2 .* x

        weightMatrix = 2 .* ones(Float32, 21+1, 1)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = Sym())
        x = randn(21,1,10)
        @test shears(x)[:,1,:,:] ≈ 2 .* x

        weightMatrix = 2 .* ones(Float32, 21>>1+1, 1)
        shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                         plan=true, boundary = FourierFilterFlux.Periodic())
        x = randn(Float32,21,1,10)
        @test shears(x)[:,1,:,:] ≈ 2 .* x
    end

    using FourierFilterFlux:applyWeight, applyBC, internalConvFFT
    weight = 2 .* ones(Complex{Float32}, (21+10)>>1+1, 1)
    bc = Pad(5)
    x = randn(Float32,21,1,10)
    xbc, usedInds = applyBC(x, bc, 1)
    x̂ = rfft(xbc, (1,))
    fftPlan = plan_rfft(xbc, (1,))
    An = map(x->FourierFilterFlux.NonAnalyticMatching(), (1:size(weight)[end]...,))
    nextLayer = internalConvFFT(x̂, weight, usedInds, fftPlan,
                                nothing, An);
    ∇ = gradient((x̂) -> internalConvFFT(x̂, weight, usedInds, fftPlan,
                                        nothing, An)[1,1,1,1,1],
                 x̂)
    y,∂ = pullback((x̂) -> internalConvFFT(x̂, weight, usedInds, fftPlan,
                                        nothing, An)[1,1,1,1,1],
                   x̂)
    ∂(y)
    ∂(y) # repeated calls to the derivative were causing errors while argWrapper
    # was in use
    @test minimum(abs.(∇[1][:, 1,1]).≈ 2f0/31)
    # no bias, not analytic and real valued output

    # no bias, analytic (so complex valued)
    fftPlan = plan_fft(xbc, (1,))
    ∇ = gradient((x̂) -> abs(applyWeight(x̂,weight[:,1], usedInds, fftPlan,
                                        nothing, FourierFilterFlux.AnalyticWavelet())[1,1,1,1]),
                 x̂)
    @test minimum(abs.(∇[1][:, 1,1]).≈ 2f0/31)

    # no bias, not analytic, complex valued, but still symmetric
    fftPlan = plan_fft(xbc, (1,))
    ∇ = gradient((x̂) -> real(applyWeight(x̂,weight[:,1], usedInds, fftPlan,
                                         nothing, FourierFilterFlux.RealWaveletRealSignal())[1,1,1,1]),
                 x̂)
    @test minimum(abs.(∇[1][2:end, 1,1]).≈ 2*2f0/31)
    @test abs(∇[1][1, 1,1]) ≈ 2f0/31

    # internal methods tests
    weightMatrix = 2 .* ones(Float32, (21+10)>>1+1, 1)
    padding = 5
    shears = ConvFFT(weightMatrix, nothing, originalSize, identity,
                    plan=true, boundary = Pad(padding))
    x = randn(21,1,10)
    x̂ = pad(x, shears.bc.padBy);
    x̂ = shears.fftPlan *ifftshift(x̂,(1,2))
    usedInds = (shears.bc.padBy[1] .+ (1:size(x, 1)),)
    nextLayer = FourierFilterFlux.internalConvFFT(x̂, shears.weight, usedInds,
                                                shears.fftPlan,
                                                shears.bias, shears.analytic);
    ∇ = gradient((x̂)->FourierFilterFlux.internalConvFFT(x̂,
                                                        shears.weight,
                                                        usedInds,
                                                        shears.fftPlan,
                                                        shears.bias, shears.analytic)[1,1,1,1,1],
                x̂)
    @test minimum(abs.(∇[1][:, 1,1]).≈ 2f0/31)
    #

    # no bias, not analytic and real valued output
    # no bias, analytic (so complex valued)
    # no bias, not analytic, complex valued, but still symmetric
    # biased (and one of the others, doesn't matter which)

    ∇ = gradient((x̂) -> (shears.fftPlan \ (x̂ .* shears.weight[:,:,1]))[1,1,1,1], x̂)
    @test minimum(abs.(∇[1][:, :, 1,1]).≈ 1f0/31*2)
    sheared = shears(x)
    @test size(sheared) == (21,1,1,10)

    # convert to a gpu version
    if CUDA.functional()
        gpuVer = shears |> gpu
        @test typeof(gpuVer.weight) <: CuArray
        @test typeof(gpuVer.fftPlan) <: CUFFT.rCuFFTPlan
    end
    # extra channel dimension
    originalSize = (20,16, 1,10)
    shears = ConvFFT(randn(Float32, 16, 3), nothing, originalSize, abs,
                    plan=true, boundary = Pad(5), trainable=false)
    @test size(shears.fftPlan)== originalSize .+ (10, 0, 0, 0)
    @test shears.σ == abs
    @test shears.bias == nothing
    @test shears.bc.padBy == (5,)
    wSize = ((originalSize[1] + (10)) >>1 +1, 3)
    @test size(shears.weight) == wSize

    # random initialization
    originalSize = (20, 16, 1, 10)
    shears = ConvFFT(originalSize, 3, abs,
                    plan=true, boundary = Pad(5), nConvDims=1)
    @test size(shears.fftPlan)== originalSize .+ (10, 0, 0, 0)
    @test shears.σ == abs
    @test size(shears.bias) == (originalSize[2:end-1]..., 3)
    @test shears.bc.padBy == (5,)
    wSize = ((originalSize[1] + 10)>>1 + 1, 3)
    @test size(shears.weight) == wSize
    end
end


text:
# Hessenberg factorization

if VERSION < v"1.3"
    struct HessenbergFactorization{T, S<:StridedMatrix, U} <: Factorization{T}
        data::S
        τ::Vector{U}
    end
    function HessenbergFactorization(A::S, τ::AbstractVector{U}
                                     ) where {S <: AbstractMatrix{T}, U} where T
        HessenbergFactorization{T,S,U}(A, τ)
    end

    Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)

    function _hessenberg!(A::StridedMatrix{T}) where T
        n = LinearAlgebra.checksquare(A)
        τ = Vector{Householder{T}}(undef, n - 1)
        for i = 1:n - 1
            xi = view(A, i + 1:n, i)
            t  = _reflector!(xi)
            H  = Householder{T,typeof(xi)}(view(xi, 2:n - i), t)
            τ[i] = H
            lmul!(H', view(A, i + 1:n, i + 1:n))
            rmul!(view(A, :, i + 1:n), H)
        end
        return HessenbergFactorization{T, typeof(A), eltype(τ)}(A, τ)
    end

    function _materializeQ(H::HessenbergFactorization{T}) where T
        n = size(H.data, 1)
        Z = Matrix{T}(I, n, n)
        for j=n-1:-1:1
            lmul!(H.τ[j], view(Z, j+1:n, j:n))
            Z[1:j-1,j] .= 0
        end
        Z
    end
else # version v1.3+

function _hessenberg!(A::StridedMatrix{T}) where T
    n = LinearAlgebra.checksquare(A)
    τ = Vector{T}(undef, n - 1)
    for i = 1:n - 1
        xi = view(A, i + 1:n, i)
        t  = _reflector!(xi)
        H  = Householder{T,typeof(xi)}(view(A, i+2:n, i), t)
        τ[i] = H.τ
        lmul!(H', view(A, i + 1:n, i + 1:n))
        rmul!(view(A, :, i + 1:n), H)
    end
    return Hessenberg(A, τ)
end

    function _hessenberg!(Ah::Union{Hermitian{T}, Symmetric{T}}) where T <: Real
        _hehessenberg!(Hermitian(Ah), Val(Symbol(Ah.uplo)))
    end
    function _hessenberg!(Ah::Hermitian{T})  where T <: Complex
        _hehessenberg!(Ah, Val(Symbol(Ah.uplo)))
    end

    function _hehessenberg!(Ah::Hermitian{T}, ::Val{:L}) where T
        # based on LAPACK zhetd2
        A = Ah.data
        n = LinearAlgebra.checksquare(A)
        τ = Vector{T}(undef, n - 1)
        d = zeros(real(T), n)
        e = zeros(real(T), n-1)
        A[1,1] = real(A[1,1])
        for i = 1:n - 1
            # reflector to annihilate A[i+2:n,i]
            ξ = view(A, i+1:n, i)
            τi  = _reflector!(ξ)
            e[i] = real(A[i+1,i])
            if !iszero(τi)
                ξ[1] = one(T)
                # store x = τᵢ A ξ in tail of τ
                w = view(τ,i:n-1)
                mul!(w, Hermitian(view(A, i+1:n, i+1:n),:L), ξ, τi, false)
                # compute w = ξ - (τᵢ/2) (xᴴ ξ) ξ
                α = - (τi / 2) * dot(w, ξ)
                w .= w .+ α .* ξ
                # apply transformation as rank-2 update
                rank2update!(view(A, i+1:n, i+1:n), -1, ξ, w, Val(:L))
            else
                A[i+1, i+1] = real(A[i+1, i+1])
            end
            A[i+1, i] = e[i]
            d[i] = real(A[i,i])
            τ[i] = τi
        end
        d[n] = real(A[n,n])
        return Hessenberg(A, τ, SymTridiagonal(d, e), Ah.uplo)
    end

    function _hehessenberg!(Ah::Hermitian{T}, ::Val{:U}) where T
        # based on LAPACK zhetd2
        A = Ah.data
        n = LinearAlgebra.checksquare(A)
        τ = Vector{T}(undef, n - 1)
        d = zeros(real(T), n)
        e = zeros(real(T), n-1)
        A[n,n] = real(A[n,n])
        xi = zeros(T,n)
        for i=n-1:-1:1
            # reflector to annihilate A[1:i-1,i+1]
            xi[1] = A[i,i+1]
            xi[2:i] .= A[1:i-1,i+1]
            ξ = view(xi,1:i)
            τi = _reflector!(ξ)
            e[i] = real(xi[1])
            A[1:i-1,i+1] .= xi[2:i]
            if !iszero(τi)
                A[i,i+1] = one(T)
                ξ = view(A, 1:i, i+1)
                # store x = τᵢ A ξ in head of τ
                w = view(τ,1:i)
                mul!(w, Hermitian(view(A, 1:i, 1:i),:U), ξ, τi, false)
                # compute w = ξ - (τᵢ/2) (xᴴ ξ) ξ
                α = - (τi / 2) * dot(w, ξ)
                w .= w .+ α .* ξ
                # apply transformation as rank-2 update
                rank2update!(view(A, 1:i, 1:i), -1, ξ, w, Val(:U))
            else
                A[i, i] = real(A[i, i])
            end
            A[i, i+1] = e[i]
            d[i+1] = real(A[i+1,i+1])
            τ[i] = τi
        end
        d[1] = real(A[1,1])
        return Hessenberg(A, τ, SymTridiagonal(d, e), Ah.uplo)
    end

    function rank2update!(A::StridedMatrix, α::Number, x::AbstractVector, y::AbstractVector, ::Val{:U})
        n = length(x)
        for j=1:n
            if iszero(x[j]) && iszero(y[j])
                A[j,j] = real(A[j,j])
            else
                t1 = α * conj(y[j])
                t2 = conj(α * x[j])
                for i=1:j-1
                    A[i,j] += x[i] * t1 + y[i] * t2
                end
                A[j,j] = real(A[j,j]) + real(x[j] * t1 + y[j] * t2)
            end
        end
    end
    function rank2update!(A::StridedMatrix, α::Number, x::AbstractVector, y::AbstractVector, ::Val{:L})
        n = length(x)
        for j=1:n
            if iszero(x[j]) && iszero(y[j])
                A[j,j] = real(A[j,j])
            else
                t1 = α * conj(y[j])
                t2 = conj(α * x[j])
                A[j,j] = real(A[j,j]) + real(x[j] * t1 + y[j] * t2)
                for i=j+1:n
                    A[i,j] += x[i] * t1 + y[i] * t2
                end
            end
        end
    end

using LinearAlgebra: QRPackedQ
function _materializeQ(H::Hessenberg{T}) where {T}
    if H.uplo == 'L'
        return _materializeQ(H, Val(:L))
    else
        return _materializeQ(H, Val(:U))
    end
end

function _materializeQ(H::Hessenberg{T}, ::Val{:L}) where {T}
    H.uplo == 'L' || throw(ArgumentError("only implemented for uplo='L'"))
    A = copy(H.Q.factors)
    n = checksquare(A)
    # shift reflectors one column rightwards
    @inbounds for j in n:-1:2
        A[1:j-1,j] .= zero(T)
        for i in j+1:n
            A[i,j] = A[i,j-1]
        end
    end
    A[2:n,1] .= zero(T)
    A[1,1] = one(T)
    Q1 = QRPackedQ(view(A,2:n,2:n), H.Q.τ)
    A[2:n,2:n] .= Matrix{T}(Q1)
    A
end

# alas, stdlib has no QLPackedQ
function _materializeQ(H::Hessenberg{T},::Val{:U}) where {T}
    H.uplo == 'U' || throw(ArgumentError("only implemented for uplo='U'"))
    A = copy(H.Q.factors)
    n = checksquare(A)
    # shift reflectors one column leftwards
    @inbounds for j in 1:n-1
        for i in 1:j-1
            A[i,j] = A[i,j+1]
        end
        A[j:n,j] .= zero(T)
    end
    A[1:n-1,n] .= zero(T)
    A[n,n] = one(T)
    w = zeros(T,n-1)
    for i in 1:n-1
        A[i,i] = one(T)
        τi = H.τ[i]
        v = view(A,1:i,i)
        # w = A[1:i,1:i-1]' * v
        for j in 1:i-1
            w[j] = dot(view(A, 1:i, j), v)
        end
        # A[1:i,1:i-1] .-= τi * v * w'
        for j in 1:i-1
            t = τi * w[j]'
            for k in 1:i
                A[k,j] -= v[k] * t
            end
        end
        A[1:i-1, i] *= -τi
        A[i,i] = one(T) - τi
        A[i+1:n,i] .= zero(T)
    end
    A
end

end # v1.3+ branch

LinearAlgebra.hessenberg!(A::StridedMatrix{<: STypes}) = _hessenberg!(A)

using LinearAlgebra: RealHermSymComplexHerm
# this should be S <: StridedMatrix, but for stdlib foolishness
function LinearAlgebra.hessenberg!(Ah::RealHermSymComplexHerm{<: STypes, <: AbstractMatrix})
    _hehessenberg!(Hermitian(Ah), Val(Symbol(Ah.uplo)))
end


text:
#=
    The file contains radicals related code
    We use the algorithm provided by basis.jl
=#

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

include("basis.jl")

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

import Nemo: base_ring, terms, monomial, derivative, gens, kernel, MatrixSpace

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

HIT = 0
Algebra = 1

function find_radical_1(vectors)
    global Algebra
    Algebra = find_basis(deepcopy(vectors))

    As = basis(Algebra)
    F = base_ring(first(As))
    n = dim(Algebra)

    # println("algebra = \n", algebra)

    traces = Dict{Tuple{Int, Int}, elem_type(F)}(
        (i, j) => tr(As[i] * As[j])
        for i in 1 : n
            for j in i : n
    )

    nnz_coords = [
        (i, j, traces[min(i, j), max(i, j)])
        for i in 1 : n
            for j in 1 : n
                if ! iszero(traces[min(i, j), max(i, j)])
    ]

    ZZ = GF(2^31 - 1)
    A = from_COO(n, n, nnz_coords, F)
    A = modular_reduction(A, ZZ)

    @info "$n×$n-dim algebra matrix of densitu $(density(A))"
    # println(AAAAAAA)

    # wiedemannchik.jl
    char_poly = minimal_polynomial(A, subspace_minpoly=__deterministic_min_poly)
    @info "minimal poly of degree $(degree(char_poly))"
    ok = iszero(evaluate(char_poly, A))
    println("m(AAAAAAA) = $(ok)")
    global HIT
    HIT += ok

    println("f = ", char_poly)

    while iszero(coeff(char_poly, 0))
        char_poly = shift_right_x(char_poly)
    end

    println("shifted f = ", char_poly)

    Image = evaluate(char_poly, A)

    Image = rational_reconstruction(Image)

    println(Image)
    radical_basis = []
    for i in [1, 3, 6, 7, 10]
        vectors = []
        for j in 1 : size(Image, 1)
            if !iszero(Image[i, j])
                push!(vectors, scale(As[j], Image[i, j]))
            end
        end
        summator(x, y) = reduce(x, y, 1)
        push!(radical_basis, to_dense(reduce(summator, vectors)))
    end

    stacked = vcat(radical_basis...)
    println(stacked)
    println(size(stacked))
    Space = MatrixSpace(QQ, size(stacked)...)
    mutual_kernel = kernel(Space(stacked))

    println("kernel =",)
    println(mutual_kernel)

end

# find_radical_1([A1, A2, A3])
#------------------------------------------------------------------------------

function uwu()
    for (i, (mfn, mdim, msz, mdata)) in enumerate(load_COO_if(from_dim=8, to_dim=15))

        @info "$i-th model : $mfn of dim : $mdim"

        if i in [2, 29, 30]
            continue
        end

        As = map(matr -> from_COO(matr..., QQ), mdata)

        @time Rad = find_radical_1(As)

        println("-------------------------")

        #break
    end
end

function OWO()
    As = [
        from_dense([
            1 0 0;
            0 2 0;
            0 0 3;
        ], QQ)
    ]

    find_radical_1(As)
end

# uwu()
# OWO()

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

#=
function construct_jacobians(system)
    nnz = [
        (i, var_index(v), derivative(f, v))
        for (i, f) in enumerate(system)
            for v in vars(f)
    ]
    m, n = length(system), length(degrees(first(system)))
    ring = parent(first(system))
    return from_COO(m, n, nnz, ring)
end
=#

function construct_jacobians(system)
    domain = base_ring(first(system))
    poly_ring = parent(first(system))

    jacobians = Dict()

    for (p_idx, poly) in enumerate(system)
        for (v_idx, var) in enumerate(vars(poly))
            # term is of form α*monomial
            for term in terms(derivative(poly, var))
                monom = monomial(term, 1)
                cf = coeff(term, 1)

                # how can we do it in a normal way..
                !haskey(jacobians, monom) && (jacobians[monom] = Dict())
                if !haskey(jacobians[monom], (p_idx, v_idx))
                    jacobians[monom][(p_idx, v_idx)] = zero(domain)
                end
                jacobians[monom][(p_idx, v_idx)] += cf
            end
        end
    end

    m, n = length(system), length(gens(poly_ring))
    factors = [
        from_COO(m, n, jac, domain)
        for jac in values(jacobians)
    ]

    @info "constructed a set of $(length(factors)) matrices $m×$n from the system Jacobian"

    return factors
end

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

#=
S, (x, y) = QQ["x", "y"]
f1 = x*y + 1
f2 = 2x
s = [f1, f2]

Js = construct_jacobians(s)
=#

#  2-th model : BIOMD0000000050.json of dim : 14 !
#  29-th model : MODEL1112150000.json of


text:
import brainflow

brainflow.enable_dev_board_logger()

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))
data_first_channel = data[eeg_channels[1], :]

# returns tuple of wavelet coeffs and lengths
wavelet_data = brainflow.perform_wavelet_transform(data_first_channel, "db4", 2)
restored_wavelet_data = brainflow.perform_inverse_wavelet_transform(wavelet_data, length(data_first_channel), "db4", 2)

fft_data = brainflow.perform_fft(data_first_channel)
restored_fft_data = brainflow.perform_ifft(fft_data)

println("Original Data")
println(data_first_channel)
println("Restored from Wavelet Data")
println(restored_wavelet_data)
println("Restored from FFT Data")
println(restored_fft_data)

max_stars_repo_name: kusumikakd/Brainflow

max_issues_repo_name: kusumikakd/Brainflow

max_forks_repo_name: kusumikakd/Brainflow


text:
module TestSecondDerivative
    using Base.Test
    import FiniteDiff

    negative_sin(x) = -sin(x)
    negative_cos(x) = -cos(x)
    square(x) = x * x
    two(x) = 2.0

    function check_error(y_true, y_approximate, err = 1 // 10_000)
        # TODO: Decide how to handle very small inputs.
        if y_true > eps(y_true)
            @test abs(y_true - y_approximate) / abs(y_true) < err
        end
    end

    function run_tests(n::Integer)
        output = Array(Float64, 1)

        funcs = (
            (sin, negative_sin),
            (cos, negative_cos),
            (square, two),
            (exp, exp),
        )

        for (f, f′′) in funcs
            for itr in 1:n
                x = 1.0 + itr

                y = FiniteDiff.second_derivative(f, x)
                check_error(f′′(x), y)

                FiniteDiff.second_derivative!(output, f, x)
                check_error(f′′(x), output[1])

                tmp = FiniteDiff.second_derivative(f, mutates=false)
                y = tmp(x)
                check_error(f′′(x), y)

                tmp! = FiniteDiff.second_derivative(f, mutates=true)
                tmp!(output, x)
                check_error(f′′(x), output[1])
            end
        end
    end

    @testset "second_derivative tests" begin
        run_tests(100)
    end
end


text:
macro something(args...)
    something_impl(collect(args))
end

function something_impl(args)
    ex = quote
        throw(ArgumentError("All evaluated as `nothing`s."))
    end
    for x in reverse(args)
        ex = quote
            x = $(esc(x))
            if x !== nothing
                x
            else
                $ex
            end
        end
    end
    return ex
end


text:
# # Instructions for Assignment02
#
# ## Introduction
#
# In this assignment,
# you will start to put the pieces together.
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line,
# or use the VS Code julia extension to run individual lines.

# The following code is used for set up;
# there is no need to change it.
# If you _do_ change it to figure out what it's doing,
# make sure you change it back or your assignment may not work correctly

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment02.html #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment02_code.html #jl

using Random
Random.seed!(42)

function generate_sequence(len)
    seq = join(rand(['A','C','G','T'], len))
    println("Your sequence is: ", seq)
    return seq
end

my_seq = generate_sequence(20)

# ## Question 1
#
# The code above generates a random 20nt DNA sequence,
# and assigns it to the variable `my_seq`.
# What kind of data type is `my_seq` (eg Float64, Int64, String, or something else)?
# Assign the variable `question1` to the correct type.
#
# Hint1: rather than guess, you can just use the `typeof()` function.

question1 = ""

# Hint2: If you're having trouble,
# note that your answer should not contain quotes.
# Eg. if you thing that `seq` has the type `MyType`,
# your code should read `question1 = MyType` and not `question1 = "MyType"`

# ## Question 2
#
# The function bellow has a "doc string",
# which explains what the function should do.
# Many julia functions have doc strings,
# which makes it easy to get help right from the REPL
#
# Open a julia REPL and type `?`.
# You will see the `julia>` prompt change to `help?>`.
# Now type `println` and hit enter.
# The docstring of `println` is printed
# and you're returned to a julia prompt
#
# Evaluate the following function, including the docstring, in the REPL
# (that is, copy from the tripple quotes down to the "end"
# and paste it into the REPL, then hit enter)

"""
    question2(sequence)

Checks if `sequence` is a String.

Example
≡≡≡≡≡≡≡≡≡

    julia> question2("hello")
    true

    julia> question2(1001)
    false
"""
function question2()
    ## put your code here
end

# At the moment, this function doesn't do what it says it should.
# Actually, it doesn't do anything.
# Fix the function so that it takes one argument
# and returns `true` if the argument is a String,
# and returns `false` otherwise.
# You should be able to run the examples in the docstring
# and get the correct answer,
# and to run this function using `my_seq` as the argument
# (it should return `true`).
#
# Hint: You haven't explicitly encountered a way to check if a type is another type.
# Learning to search for answers is a key programming skill!
# Try searching "check if type is string julia".
# Typically, when I see search results, I prioritize in this order:
# 1. stackoverflow.com
# 2. docs.julialang.org
# 3. discourse.julialang.org
# 4. github.com

# ## Question 3
#
# Let's try something a bit more complicated.
# The `question3()` function has a doc string,
# but is also incomplete.
# I've written some code to get you started,
# complete it so that the examples work as expected in the docstring.
# If you use this function on `my_seq`,
# you should get the answer `0.3`

"""
    question3(sequence)

Calculates the GC ratio of a DNA sequence.
The GC ratio is the total number of G and C bases divided by the total length of the sequence.
For more info about GC content, see here:

Example
≡≡≡≡≡≡≡≡≡≡

    julia> question3("AATG")
    0.25

    julia> question3("CCCGG")
    1.0

    julia> question3("ATTA")
    0.0
"""
function question3(sequence)
    ## throw an error if the string contains anything other than ACGT
    if any(c-> !in(c, ['A','C','G','T']), sequence)
        throw(ArgumentError("Sequence must only contain ACGT"))
    end

    ## change line to assign `seqlength` to the length of `sequence` instead of `1`
    ## If you're stuck, search for "length of string julia"
    seqlength = 1

    ## count the number of G's
    gs = count(==('G'), sequence)
    ## count the number of C's
    cs = count(==('C'), sequence)

    return gs + cs / seqlength # something is wrong with this line...
end

# ## Question 4
#
# This process should be familiar by now.
#
# Hint: you do not need to re-write code if you already have it.
# That is, don't copy any code from `question3`,
# just call it!

"""
    question4(sequence)

Calculates the GC content of a DNA sequence
and prints it to the screen.

Example
≡≡≡≡≡≡≡≡≡≡

    julia> question4("AATC")
    Sequence:
    AATC
    GC Content:
    0.25

    julia> question4("CCCGG")
    Sequence:
    CCCGG
    GC Content:
    1.0
"""
function question4(sequence)
    ## Your code here
end


text:
# Sum function: f(x)=sum(x)
Fsum(x)=([sum(x)],nothing)
Fsum_inplace(value,auxvalue,x)=fill!(value,sum(x))

function Dsum(derivativeIDX,f_c,faux_c,grad_c,grad_n,x_n...)
    axpy!(grad_c[1],ones(size(grad_n)),grad_n)
end

function Fsum(x...)
    tmp=0.0
    for i in 1:length(x)
        tmp+=sum(x[i])
    end
    return ([tmp],nothing)
end

function Fsum_inplace(value::Array,auxvalue,x...) # inplace
    tmp=0.0
    for i in 1:length(x)
        tmp+=sum(x[i])
    end
    fill!(value,tmp)
end

import Base.sum

if PROC=="GPU"

    function sum(A::CudaArray)
  #      return flatten(Float64,CUBLAS.gemv('T',1.,flatten(Float64,A),CudaArray(ones(length(A)))))
        return flatten(CUBLAS.gemv('T',convert(eltype(A),1),flatten(A),CudaArray(eltype(A),ones(length(A)))))
    end

    export sum

    import Base.sum!
    function sum!(out::CudaArray,A::CudaArray)
        CUBLAS.gemm!('T','N',1.,flatten(Float64,A),CudaArray(ones(length(A),1)),0.0,out)
    end
    export sum!

    function sum_update!(alpha::Float64,A::CudaArray,beta::Float64,out::CudaArray)
        CUBLAS.gemm!('T','N',alpha,flatten(Float64,A),CudaArray(ones(length(A),1)),beta,out)
    end
    export sum_update!

    function Dsum(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray...)
        tmp=CudaArray(eltype(grad_n),size(grad_n))
        fill!(tmp,convert(eltype(grad_n),1.0))
        axpy!(grad_c,tmp,grad_n)
        free(tmp)
    end

    function Fsum(x::CudaArray...)
        tmp=CudaArray(zeros(1,1))
        for i in 1:length(x)
            axpy!(1.0,sum(x[i]),tmp)
        end
        return (tmp,nothing)
    end

    function Fsum_inplace(value::CudaArray,auxvalue,x::CudaArray...) # inplace
        fill!(value,0.0)
        for i in 1:length(x)
            axpy!(1.0,sum(x[i]),value)
        end
    end

end

Derivative[Fsum]=Dsum # Define dictionary lookup
Inplace[Fsum]=Fsum_inplace

import Base.sum
sum(n::ADnode)=ADnode(Fsum,n)

function sum(n::ArrayADnode)
    return ADnode(Fsum,n)
end

# TODO: add similar as below for each unary function that can take a transposed argument
#sum(A::ADtrans)=ADnode(Fsum, ftranspose(node[A.parent]))
sum(A::ADtrans)=ADnode(Fsum, node[A.parent]) # sum(A')=sum(A)

sum(A::ADdiag)=sum(node[A.parent])

export sum

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# AUTO GENERATED FILE - DO NOT EDIT

export switch

"""
    switch(;kwargs...)

A Switch component.
Capture user feedback limited to small set of options. For more information, see: https://mantine.dev/core/switch/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `checked` (Bool; optional): State of check box
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Checkbox color
- `disabled` (Bool; optional): A checkbox can show it is currently unable to be interacted with
- `label` (String; optional): Checkbox label
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined border-radius value from theme.radius or number for border-radius in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined label font-size and checkbox width and height in px
- `style` (Dict; optional): Inline style override
"""
function switch(; kwargs...)
        available_props = Symbol[:id, :checked, :className, :color, :disabled, :label, :radius, :size, :style]
        wild_props = Symbol[]
        return Component("switch", "Switch", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
using Media, Lazy, Hiccup

import Media: render
import Hiccup: div

import Juno: Inline, Clipboard, Editor, Console

Media.@defpool Editor
Media.@defpool Console

setdisplay(Editor(), Any, Console())
setdisplay(Console(), Any, Console())

include("plots.jl")
include("view.jl")
include("lazy.jl")
include("errors.jl")


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

include(joinpath("$(@__DIR__)", "..", "example", "FMI2", "BouncingBall", "src", "BouncingBall.jl"))
@test isfile(fmu_save_path)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   J2 orbit propagator algorithm.
#
#       This algorithm propagates the orbit considering the perturbed two-body
#       equations as presented in [1, p. 690-692]. It uses the first-order
#       approximation of Kepler's problem, considering the effects of secular
#       gravitational and drag perturbations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Wertz, J. R (1978). Spacecraft attitude determination and control.
#       Kluwer Academic Publishers, Dordrecht, The Netherlands.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export j2_gc_egm08, j2_gc_egm96, j2_gc_jgm02, j2_gc_jgm03
export j2_init, j2!

################################################################################
#                                     TODO
################################################################################
#
# 1. Analyze the reference frame representation of the inputs for this
#    algorithm.
#
#   The SGP4 algorithm expects that the input parameters are represented in the
#   TEME (true equator, mean equinox) reference frame. This J2 orbit propagator
#   model requires that the input parameters are consistent with the
#   gravitational perturbation theory in which the `J2` coefficient was
#   computed. Looking at [1, p. 642], it appears that the perturbations are
#   considering a frame in which the Z-axis is aligned with the CIP (Celestial
#   Intermediate Pole, or the Earth rotation axis). Hence, the J2 parameter is
#   defined based on the PEF. Since no rotations or adaptations are programmed,
#   then the input parameters for this propagator should be represented in any
#   reference frame with a true Equator, because of the symmetry.
#
#   This need to be further analyzed and confirmed.
#
################################################################################

################################################################################
#                                  Overloads
################################################################################

# Copy for J2_Structure
Base.copy(m::J2_Structure) = J2_Structure([ getfield(m, k) for k = 1:length(fieldnames(m)) ]...)

# Deepcopy for J2_Structure.
Base.deepcopy(m::J2_Structure) = J2_Structure([ deepcopy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...)

################################################################################
#                                  Constants
################################################################################

# These constants were obtained from the GFC files. Remember that:
#
#   J_n = -C_n,0 * sqrt(2n+1)
#

# EGM-08 Gravitational constants.
const j2_gc_egm08 = J2_GravCte(
        6378137.0,
        sqrt(3.986004415e14/6378137.0^3),
        0.0010826261738522227
       )

# EGM-96 Gravitational constants.
const j2_gc_egm96 = J2_GravCte(
        6378136.3,
        sqrt(3.986004415e14/6378136.3^3),
        0.0010826266835531513
       )

# JGM-02 Gravitational constants.
const j2_gc_jgm02 = J2_GravCte(
        6378136.3,
        sqrt(3.986004415e14/6378136.3^3),
        0.0010826269256388149
       )

# JGM-03 Gravitational constants.
const j2_gc_jgm03 = J2_GravCte(
        6378136.3,
        sqrt(3.986004415e14/6378136.3^3),
        0.0010826360229829945
       )

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

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

Initialize the data structure of J2 orbit propagator algorithm.

# Args

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

# Returns

The structure `J2_Structure` with the initialized parameters.

# Remarks

The inputs are the mean orbital elements.

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

    # Unpack the gravitational constants to improve code readability.
    @unpack_J2_GravCte j2_gc

    # Initial values.
    al_0 = a_0/R0            # Normalized semi-major axis [er].
    n_0  = μm/al_0^(3/2)     # Unperturbed mean motion [rad/s].
    p_0  = al_0*(1-e_0^2)    # Semi-latus rectum [er].
    M_0  = f_to_M(e_0, f_0)  # Initial mean anomaly [rad].

    # Auxiliary variables.
    dn   = 2dn_o2            # Time-derivative of the mean motion [rad/s²].
    p_0² = p_0^2
    e_0² = e_0^2

    sin_i_0, cos_i_0 = sincos(i_0)
    sin_i_0² = sin_i_0^2

    # First-order time-derivative of the orbital elements.
    #
    # See [1, p 692].

    δa   = +2/3*al_0*dn/n_0
    δe   = +2/3*(1-e_0)*dn/n_0
    δΩ   = -3/2*n_0*J2/p_0²*cos_i_0
    δω   = +3/4*n_0*J2/p_0²*(4 - 5sin_i_0²)
    δM_0 = +3/4*n_0*J2/p_0²*sqrt(1-e_0²)*(2 - 3sin_i_0²)

    # Create the output structure with the data.
    J2_Structure{T}(
        epoch, al_0, n_0, e_0, i_0, Ω_0, ω_0, f_0, M_0, 0, dn_o2, ddn_o6, al_0,
        e_0, i_0, Ω_0, ω_0, f_0, M_0, δa, δe, δΩ, δω, δM_0, j2_gc
    )
end

"""
    function j2!(j2d::J2_Structure{T}, t::Number) where T

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

# Returns

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

###### Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. If the orbit parameters are obtained
from a TLE, then the inertial frame will be TEME. Notice, however, that the
perturbation theory requires an inertial frame with true equator.

"""
function j2!(j2d::J2_Structure{T}, t::Number) where T
    # Unpack the variables.
    @unpack_J2_Structure j2d
    @unpack_J2_GravCte   j2_gc

    # Time elapsed since epoch.
    Δt = t

    # Propagate the orbital elements.
    al_k = al_0 - δa*Δt
    e_k  = e_0 - δe*Δt
    i_k  = i_0
    Ω_k  = mod(Ω_0 + δΩ*Δt, 2π)
    ω_k  = mod(ω_0 + δω*Δt, 2π)

    # In [1, p. 692], it is mentioned that the mean anomaly must be updated
    # considering the equation:
    #                           .           ..
    #                          n_0          n_0
    #   M_k = M_0 + n_0 ⋅ Δt + ---- ⋅ Δt² + ---- ⋅ Δt³
    #                           2            6
    #
    # However, the time derivative of M_0 for perturbed orbits considering terms
    # up to J2 is not 0. This is clear in [1, p. 652].
    #
    # Hence, we will use the equation:
    #                                     .           ..
    #                .                    n_0          n_0
    #   M_k = M_0 + M_0 ⋅ Δt + n_0 ⋅ Δt + ---- ⋅ Δt² + ---- ⋅ Δt³
    #                                      2            6
    #

    M_k = mod(@evalpoly(Δt, M_0, (n_0 + δM_0), dn_o2, ddn_o6), 2π)
    f_k = M_to_f(e_k, M_k)

    # Make sure that eccentricity is not lower than 0.
    (e_k < 0) && (e_k = T(0))

    # Compute the position and velocity vectors given the orbital elements.
    (r_i_k, v_i_k) = kepler_to_rv(al_k*R0, e_k, i_k, Ω_k, ω_k, f_k)

    # Update the J2 orbit propagator structure.
    j2d.Δt   = Δt
    j2d.al_k = al_k
    j2d.e_k  = e_k
    j2d.i_k  = i_k
    j2d.Ω_k  = Ω_k
    j2d.ω_k  = ω_k
    j2d.M_k  = M_k
    j2d.f_k  = f_k

    # Return the position and velocity vector represented in the inertial
    # reference frame.
    (r_i_k, v_i_k)
end


text:
const IDER_TYPE = Union{AbstractString, Integer}

metindex(metnet::MetNet, ider::Int) =
    0 < ider <= length(metnet.mets) ? ider :
        error("Index ($ider) out of range, mets [1:$(length(metnet.mets))]")
metindex(metnet::MetNet, ider::AbstractString) =
    (indx = findfirst(isequal(ider), metnet.mets);
    isnothing(indx) ? error("ider ($ider) not found!!!") : indx)

rxnindex(metnet::MetNet, ider::Int) =
    0 < ider <= length(metnet.rxns) ? ider :
        error("Index ($ider) out of range, rxns [1:$(length(metnet.rxns))]")
rxnindex(metnet::MetNet, ider::AbstractString) =
    (indx = findfirst(isequal(ider), metnet.rxns);
    isnothing(indx) ? error("ider ($ider) not found!!!") : indx)


text:
module GreenFlux

using Flux: Chain, Recur, Dense, Conv, MaxPool, GlobalMaxPool, MeanPool, GlobalMeanPool,
            DepthwiseConv, ConvTranspose, CrossCor, GRUCell, LSTMCell,RNNCell, Maxout
using Flux: params, outdims, sigmoid, rrelu, elu, celu, softsign, softplus, tanh, gelu,
            hardsigmoid, logsigmoid,swish, selu, softmax, logsoftmax, hardtanh,
            leakyrelu, relu6, lisht, tanhshrink, logcosh, mish, relu, trelu,
            softshrink, identity
using Statistics: mean
using CUDAapi: has_cuda_gpu

export avgpowerdraw, modelflops

include("power/powerdraw.jl")
include("neflops/measureutils.jl")
include("neflops/layerflops.jl")
include("neflops/gradientflops.jl")
include("neflops/modelflops.jl")

function __init__()
   @info "Finished loading GreenFlux..."
end

end # module


text:
function applybranchics!(system::LegSystem,old=Dict("a"=>0),restart="no")
    if restart == "no"
        Pstart = 70*mmHgToPa;
        for i in 1:length(system.branches.ID)
            system.branches.A[i][:,1] = (Pstart/system.branches.beta[i][end] +
                sqrt(system.branches.A0[i][end]))^2;
            system.branches.Q[i][:,1] = 0;
            system.branches.P[i][:,1] = Pstart/mmHgToPa;
        end
    elseif restart == "yes"
        branches = old["branches"];
        for i in 1:length(system.branches.ID)
            if system.solverparams.JL == old["solverparams"]["JL"]
                system.branches.A[i][:,1] = branches["A"][i][:,end];
                system.branches.Q[i][:,1] = branches["Q"][i][:,end];
                system.branches.P[i][:,1] = branches["P"][i][:,end];
            else
                A = branches["A"][i][:,end];
                Q = branches["Q"][i][:,end];
                P = branches["P"][i][:,end];
                Aitp = Interpolations.interpolate(A, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
                Qitp = Interpolations.interpolate(Q, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
                Pitp = Interpolations.interpolate(P, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
                x = 0:1:(old["solverparams"]["JL"]-1);
                x = x*branches["k"][1];
                sAitp = scale(Aitp,x);
                sQitp = scale(Qitp,x);
                sPitp = scale(Pitp,x);
                xq = 0:1:(system.solverparams.JL-1);
                xq = xq*system.branches.k[i];
                Aq = [sAitp[j] for j in xq];
                Qq = [sQitp[j] for j in xq];
                Pq = [sPitp[j] for j in xq];
                system.branches.A[i][:,1] = Aq;
                system.branches.Q[i][:,1] = Qq;
                system.branches.P[i][:,1] = Pq;
            end
        end
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
"""
`prior(parameters::ParameterVector{T}) where {T<:Number}`

Calculates log joint prior density of m.parameters.
"""
function prior(parameters::ParameterVector{T}) where {T<:Number}
    free_params = Base.filter(x -> !_filter_all_fixed_para(x), parameters)
    logpdfs = map(logpdf, free_params)
    return sum(logpdfs)
end

_filter_all_fixed_para(p::Parameter) = haskey(p.regimes, :fixed) ? all(values(p.regimes[:fixed])) : p.fixed

"""
`prior(parameters::ParameterVector{T}, x::AbstractVector{T}; isfree::Bool = false) where {T<:Number}`
Calculates log joint prior density of the vector `x ∈ ℝⁿ` using
the priors and transformations specified in `parameters`,
appropriately adjusting the prior using the Jacobian of the
transformations.
If the keyword `isfree` is true, then we assume that `x` only
corresponds to free parameters already, hence its length is potentially
shorter than the length of `parameters`.
"""
function prior(parameters::ParameterVector, x::AbstractVector{T};
               isfree::Bool = false) where {T <: Number}
    if isfree
        free_params = parameters
        free_x      = x
    else
        unfixed     = map(θ -> !θ.fixed, parameters)
        free_params = Base.filter(x -> !_filter_all_fixed_para(x), parameters)
        free_x      = x[unfixed]
    end

    logpdfs = zero(T)
    @inbounds for i = 1:length(free_params)
        if hasprior(free_params[i])
            logpdfs += logpdf(get(free_params[i].prior), transform_to_model_space(free_params[i], free_x[i])) +
                log(abs(differentiate_transform_to_model_space(free_params[i], free_x[i])))
        end
    end
    logpdfs
end

"""
```
posterior(loglikelihood::Function, parameters::ParameterVector,
                   data::AbstractArray; ϕ_smc::Float64 = 1.)
```

Calculates and returns the log of the posterior distribution for `m.parameters`:

```
log posterior  = log likelihood + log prior + const
log Pr(Θ|data) = log Pr(data|Θ) + log Pr(Θ) + const
```

### Arguments

- `loglikelihood`: Function which takes vector of parameters and data, and
        outputs the log-loglikelihood of a parameters given data.
- `parameters`: Vector of parameters whose likelihood is to be evaluated
- `data`: Matrix of data for observables

### Optional Arguments

- `φ_smc`: a tempering factor to change the relative weighting of the prior and
     the likelihood when calculating the posterior. It is used primarily in SMC.
"""
function posterior(loglikelihood::Function, parameters::ParameterVector,
                   data::AbstractArray; ϕ_smc::Float64 = 1.)
    like = loglikelihood(parameters, data)
    post = ϕ_smc * like + prior(parameters)
    return post
end

"""
```
posterior!(loglikelihood::Function, parameters::ParameterVector,
           para_draw::Vector{T}, data::Matrix{T};
           sampler::Bool = false, catch_errors::Bool = false,
           φ_smc = 1) where {T<:AbstractFloat}
```

Evaluates the log posterior density at `parameters`.

### Arguments

- `loglikelihood`: a function which takes vector of parameters and data, and
                   outputs the log-loglikelihood of a parameters given data.

- `para_draw`: New values for the model parameters
- `data`: Matrix of input data for observables

### Optional Arguments
- `sampler`: Whether metropolis_hastings or smc is the caller. If `sampler=true`,
     the log likelihood and the transition matrices for the zero-lower-bound
     period are also returned.
- `catch_errors`: Whether to catch errors of type `GensysError` or `ParamBoundsError`
     If `sampler = true`, both should always be caught.
- `φ_smc`: a tempering factor to change the relative weighting of the prior and
     the likelihood when calculating the posterior. It is used primarily in SMC.
"""
function posterior!(loglikelihood::Function, parameters::ParameterVector,
                    para_draw::Vector{T}, data::Matrix{T};
                    sampler::Bool = false, catch_errors::Bool = false,
                    ϕ_smc::Float64 = 1.) where {T<:AbstractFloat}
    catch_errors = catch_errors | sampler
    if sampler
        try
            update!(parameters, para_draw)
        catch err
            if isa(err, ParamBoundsError)
                return -Inf
            else
                throw(err)
            end
        end
    else
        update!(parameters, para_draw)
    end
    return posterior(loglikelihood, parameters, data; ϕ_smc = ϕ_smc)
end

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


text:
using Documenter
using RiskMeasures

push!(LOAD_PATH,"../src/")
makedocs(sitename="RiskMeasures.jl Documentation",
         pages = [
            "Home" => "index.md",
            "Tutorials" => Any["basic.md",
				    "coherent.md",
				    "convex.md"
            ],
            "Function Library" => "library.md",
         ],
         format = Documenter.HTML(prettyurls = false)
)
# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/rubsc/RiskMeasures.jl.git",
    devbranch = "main"
)


text:
#####
##### stochastic simulation
#####

using Dyno

rootdir = Pkg.dir("Dyno")
filename = joinpath(rootdir,"examples","rbc.mod")

@time model = import_modfile(filename)
@time sol = solve(model)
@time sim = simulate(sol)
@time sim = simulate(model)

using Gadfly
plot(x=sim[:t], y=sim[:y], Geom.line)

#####
##### deterministic solve
#####

exo  = [ 0.0 0.1 0.00]'
@time sol = deterministic(model, exo, verbose=true,N=200)

using Gadfly
plot(x=sol[:t], y=sol[:y], Geom.line)