text:
function test_gray_over_gray(bm)
    @test all(c->gray(bm(c[1], c[2])) === red(bm(RGB(c[1]), RGB(c[2]))),
              ((c1, c2) for c1 in gray_f64, c2 in gray_f64))
    @test all(c->gray(bm(c[1], c[2])) === red(bm(RGB(c[1]), RGB(c[2]))),
              ((c1, c2) for c1 in gray_f32, c2 in gray_f32))
    @test all(c->gray(bm(c[1], c[2])) === red(bm(RGB(c[1]), RGB(c[2]))),
              ((c1, c2) for c1 in gray_n0f8, c2 in gray_n0f8))
    @test all(c->gray(bm(c[1], c[2])) === red(bm(RGB(c[1]), RGB(c[2]))),
              ((c1, c2) for c1 in gray_n0f16, c2 in gray_n0f16))
end

function test_gray_over_rgb(bm)
    @test all(c->bm(c[1], c[2]) === bm(c[1], RGB(c[2])),
              ((c1, c2) for c1 in rgb2_f64, c2 in gray_f64))
    @test all(c->bm(c[1], c[2]) === bm(c[1], RGB(c[2])),
              ((c1, c2) for c1 in rgb2_f32, c2 in gray_f32))
    @test all(c->bm(c[1], c[2]) === bm(c[1], RGB(c[2])),
              ((c1, c2) for c1 in rgb2_n0f8, c2 in gray_n0f8))
    @test all(c->bm(c[1], c[2]) === bm(c[1], RGB(c[2])),
              ((c1, c2) for c1 in rgb2_n0f16, c2 in gray_n0f16))
end

function test_rgb_over_gray(bm)
    @test all(c->bm(c[1], c[2]) === bm(c[1], Gray(c[2])),
              ((c1, c2) for c1 in gray_f64, c2 in rgb2_f64))
    @test all(c->bm(c[1], c[2]) === bm(c[1], Gray(c[2])),
              ((c1, c2) for c1 in gray_f32, c2 in rgb2_f32))
    @test all(c->bm(c[1], c[2]) === bm(c[1], Gray(c[2])),
              ((c1, c2) for c1 in gray_n0f8, c2 in rgb2_n0f8))
    @test all(c->bm(c[1], c[2]) === bm(c[1], Gray(c[2])),
              ((c1, c2) for c1 in gray_n0f16, c2 in rgb2_n0f16))
end

@testset "Gray" begin
    @testset "$(keyword(m))" for m in separable_modes
        @testset "Gray over Gray" begin
            test_gray_over_gray(m)
        end

        @testset "GrayA over Gray: $T" for T in Ts
            rgb = m(RGB{T}(1, 0.75, 0), RGBA{T}(0, 0.5, 1, 0.6))
            actual = (m(Gray{T}(1.0),  GrayA{T}(0.0, 0.6)),
                      m(Gray{T}(0.75), GrayA{T}(0.5, 0.6)),
                      m(Gray{T}(0.0),  GrayA{T}(1.0, 0.6)))
            expected = Gray{T}.((red(rgb), green(rgb), blue(rgb)))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "Gray over GrayA: $T" for T in Ts
            rgb = m(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}(0, 0.5, 1))
            actual = (m(GrayA{T}(1.0,  0.6), Gray{T}(0.0)),
                      m(GrayA{T}(0.75, 0.6), Gray{T}(0.5)),
                      m(GrayA{T}(0.0,  0.6), Gray{T}(1.0)))
            expected = GrayA{T}.((red(rgb), green(rgb), blue(rgb)), alpha(rgb))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "GrayA over GrayA: $T" for T in Ts
            rgb = m(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}(0, 0.5, 1, 0.6))
            actual = (m(GrayA{T}(1.0,  0.6), GrayA{T}(0.0, 0.6)),
                      m(GrayA{T}(0.75, 0.6), GrayA{T}(0.5, 0.6)),
                      m(GrayA{T}(0.0,  0.6), GrayA{T}(1.0, 0.6)))
            expected = GrayA{T}.((red(rgb), green(rgb), blue(rgb)), alpha(rgb))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end
    end
end

@testset "Gray over RGB" begin
    @testset "$(keyword(m))" for m in separable_modes
        @testset "Gray over RGB" begin
            test_gray_over_rgb(m)
        end

        @testset "GrayA over RGB: $T" for T in Ts
            grays = (GrayA{T}(0.0, 0.6), GrayA{T}(0.5, 0.6), GrayA{T}(1.0, 0.6))
            actual = m.(RGB{T}(1, 0.75, 0), grays)
            expected = m.(RGB{T}(1, 0.75, 0), RGBA{T}.(grays))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "Gray over RGBA: $T" for T in Ts
            grays = (Gray{T}(0.0), Gray{T}(0.5), Gray{T}(1.0))
            actual = m.(RGBA{T}(1, 0.75, 0, 0.6), grays)
            expected = m.(RGBA{T}(1, 0.75, 0, 0.6), RGB{T}.(grays))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "GrayA over RGBA: $T" for T in Ts
            grays = (GrayA{T}(0.0, 0.6), GrayA{T}(0.5, 0.6), GrayA{T}(1.0, 0.6))
            actual = m.(RGBA{T}(1, 0.75, 0, 0.6), grays)
            expected = m.(RGBA{T}(1, 0.75, 0, 0.6), RGBA{T}.(grays))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end
    end
end

@testset "RGB over Gray" begin
    @testset "$(keyword(m))" for m in separable_modes
        @testset "RGB over Gray" begin
            test_rgb_over_gray(m)
        end

        @testset "RGBA over Gray: $T" for T in Ts
            grays = (Gray{T}(1.0), Gray{T}(0.75), Gray{T}(0.0))
            actual = m.(grays, RGBA{T}(0.0, 0.6, 1.0, 0.6))
            expected = m.(grays, GrayA{T}(0.464, 0.6))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "RGB over GrayA: $T" for T in Ts
            grays = (GrayA{T}(1.0, 0.6), GrayA{T}(0.75, 0.6), GrayA{T}(0.0, 0.6))
            actual = m.(grays, RGB{T}(0, 0.6, 1))
            expected = m.(grays, Gray{T}(0.464))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end

        @testset "RGBA over GrayA: $T" for T in Ts
            grays = (GrayA{T}(1.0, 0.6), GrayA{T}(0.75, 0.6), GrayA{T}(0.0, 0.6))
            actual = m.(grays, RGBA{T}(0.0, 0.6, 1.0, 0.6))
            expected = m.(grays, GrayA{T}(0.464, 0.6))
            @test all(t -> t[1] ≈ t[2], zip(actual, expected))
        end
    end
end

@testset "RGB24/ARGB32/Gray24/AGray32" begin
    @test BlendHue(       RGB24(1, 0.8, 0), RGB24(0, 0.4, 1))       === RGB24(0.651, 0.792, 1.0)
    @test BlendSaturation(RGB24(1, 0.8, 0), ARGB32(0, 0.4, 1, 0.6)) === RGB24(1, 0.8, 0)
    @test BlendColor(     RGB24(1, 0.8, 0), Gray24(0.4))            === RGB24(0.773, 0.773, 0.773)
    @test BlendLuminosity(RGB24(1, 0.8, 0), AGray32(0.4, 0.6))      === RGB24(0.710, 0.569, 0.000)
    @test BlendNormal(  ARGB32(1, 0.8, 0, 0.6), RGB24(0, 0.4, 1))       === ARGB32(0, 0.4, 1, 1)
    @test BlendMultiply(ARGB32(1, 0.8, 0, 0.6), ARGB32(0, 0.4, 1, 0.6)) === ARGB32(0.286, 0.482, 0.286, 0.84)
    @test BlendScreen(  ARGB32(1, 0.8, 0, 0.6), Gray24(0.4))            === ARGB32(0.761, 0.686, 0.400, 1)
    @test BlendOverlay( ARGB32(1, 0.8, 0, 0.6), AGray32(0.4, 0.6))      === ARGB32(0.827, 0.667, 0.114, 0.84)
    @test BlendDarken(    Gray24(0.8), RGB24(0, 0.4, 1))       === Gray24(0.345)
    @test BlendLighten(   Gray24(0.8), ARGB32(0, 0.4, 1, 0.6)) === Gray24(0.8)
    @test BlendColorDodge(Gray24(0.8), Gray24(0.4))            === Gray24(1.0)
    @test BlendColorBurn( Gray24(0.8), AGray32(0.4, 0.6))      === Gray24(0.62)
    @test BlendHardLight( AGray32(0.8, 0.6), RGB24(0, 0.4, 1))       === AGray32(0.471, 1)
    @test BlendSoftLight( AGray32(0.8, 0.6), ARGB32(0, 0.4, 1, 0.6)) === AGray32(0.647, 0.84)
    @test BlendDifference(AGray32(0.8, 0.6), Gray24(0.4))            === AGray32(0.4, 1)
    @test BlendExclusion( AGray32(0.8, 0.6), AGray32(0.4, 0.6))      === AGray32(0.584, 0.84)
end


text:
"""
```
Module       FinanceLib.Statements
Description  Implement Statements modules for the FinanceLib library
Copyright    (c) 2022 Kishaloy Neogi
License      MIT
Maintainer   Kishaloy Neogi
Email        nkishaloy@yahoo.com
```
The module describes the base modules of Statements. These includes Balance Sheets,
Income Statements and Cash Flow Statements and their utiities.

You may see the github repository at <https//github.com/n-kishaloy/FinanceLib.jl>
"""
module Statements

import FinanceLib as Fl

@enum BsTyp begin
  Cash
  CurrentReceivables
  CurrentLoans
  CurrentAdvances
  OtherCurrentAssets
  CurrentInvestments
  Inventories
  RawMaterials
  WorkInProgress
  FinishedGoods
  CurrentAssets
  AccountReceivables
  LongTermLoanAssets
  LongTermAdvances
  LongTermInvestments
  OtherLongTermAssets
  PlantPropertyEquipment
  AccumulatedDepreciation
  NetPlantPropertyEquipment
  LeasingRentalAssets
  AccumulatedAmortizationLeaseRental
  NetLeaseRentalAssets
  Goodwill
  CapitalWip
  OtherTangibleAssets
  IntangibleAssets
  IntangibleAssetsDevelopment
  AccumulatedAmortization
  NetIntangibleAssets
  LongTermAssets
  Assets
  CurrentPayables
  CurrentBorrowings
  CurrentNotesPayable
  OtherCurrentLiabilities
  InterestPayable
  CurrentProvisions
  CurrentTaxPayables
  LiabilitiesSaleAssets
  CurrentLeasesLiability
  CurrentLiabilities
  AccountPayables
  LongTermBorrowings
  BondsPayable
  DeferredTaxLiabilities
  LongTermLeasesLiability
  DeferredCompensation
  DeferredRevenues
  CustomerDeposits
  OtherLongTermLiabilities
  PensionProvision
  TaxProvision
  LongTermProvision
  LongTermLiabilities
  Liabilities
  CommonStock
  PreferredStock
  PdInCapAbovePar
  PdInCapTreasuryStock
  RevaluationReserves
  Reserves
  RetainedEarnings
  AccumulatedOci
  MinorityInterests
  Equity
end

@enum PlTyp begin
  OperatingRevenue
  NonOperatingRevenue
  ExciseStaxLevy
  OtherIncome
  Revenue
  CostMaterial
  DirectExpenses
  COGS
  Salaries
  AdministrativeExpenses
  ResearchNDevelopment
  OtherOverheads
  OtherOperativeExpenses
  OtherExpenses
  ExceptionalItems
  GrossProfit
  Pbitda
  Depreciation
  AssetImpairment
  LossDivestitures
  Amortization
  Pbitx
  InterestRevenue
  InterestExpense
  CostDebt
  OtherFinancialRevenue
  Pbtx
  ExtraordinaryItems
  PriorYears
  Pbt
  TaxesCurrent
  TaxesDeferred
  Pat
  NetIncomeDiscontinuedOps
  NetIncome
  GainsLossesForex
  GainsLossesActurial
  GainsLossesSales
  FvChangeAvlSale
  OtherDeferredTaxes
  OtherComprehensiveIncome
  TotalComprehensiveIncome
end

@enum CfTyp begin
  DeferredIncomeTaxes
  ChangeInventories
  ChangeReceivables
  ChangePayables
  ChangeLiabilities
  ChangeProvisions
  OtherCfOperations
  StockCompensationExpense
  StockCompensationTaxBenefit
  AccretionDebtDiscount
  CashFlowOperations
  InvestmentsPpe
  InvestmentsCapDevp
  InvestmentsLoans
  AcqEquityAssets
  DisEquityAssets
  DisPpe
  ChangeInvestments
  CfInvestmentInterest
  CfInvestmentDividends
  OtherCfInvestments
  CashFlowInvestments
  StockSales
  StockRepurchase
  DebtIssue
  DebtRepay
  InterestFin
  Dividends
  DonorContribution
  OtherCfFinancing
  CashFlowFinancing
  NetCashFlow
end

balanceSheetCalcMap =
  [
    (Inventories,
    [RawMaterials, WorkInProgress, FinishedGoods],
    BsTyp[]
    )
  , (CurrentAssets,
    [Cash, CurrentReceivables, CurrentLoans, CurrentAdvances, OtherCurrentAssets, CurrentInvestments, Inventories],
    BsTyp[]
    )
  , (NetPlantPropertyEquipment,
    [PlantPropertyEquipment],
    [AccumulatedDepreciation]
    )
  , (NetLeaseRentalAssets,
    [LeasingRentalAssets],
    [AccumulatedAmortizationLeaseRental]
    )
  , (NetIntangibleAssets,
    [IntangibleAssets, IntangibleAssetsDevelopment],
    [AccumulatedAmortization]
    )
  , (LongTermAssets,
    [AccountReceivables, LongTermLoanAssets, LongTermAdvances, LongTermInvestments, OtherLongTermAssets, NetPlantPropertyEquipment, NetLeaseRentalAssets, Goodwill, CapitalWip, OtherTangibleAssets, NetIntangibleAssets],
    BsTyp[]
    )
  , (Assets,
    [CurrentAssets, LongTermAssets],
    BsTyp[]
    )
  , (CurrentLiabilities,
    [CurrentPayables, CurrentBorrowings, CurrentNotesPayable, OtherCurrentLiabilities, InterestPayable, CurrentProvisions, CurrentTaxPayables, LiabilitiesSaleAssets, CurrentLeasesLiability],
    BsTyp[]
    )
  , (LongTermLiabilities,
    [AccountPayables, LongTermBorrowings, BondsPayable, DeferredTaxLiabilities, LongTermLeasesLiability, DeferredCompensation, DeferredRevenues, CustomerDeposits, OtherLongTermLiabilities, PensionProvision, TaxProvision, LongTermProvision],
    BsTyp[]
    )
  , (Liabilities,
    [CurrentLiabilities, LongTermLiabilities],
    BsTyp[]
    )
  , (Equity,
    [CommonStock, PreferredStock, PdInCapAbovePar, PdInCapTreasuryStock, RevaluationReserves, Reserves, RetainedEarnings, AccumulatedOci, MinorityInterests],
    BsTyp[]
    )

  ]  :: Vector{Tuple{BsTyp, Vector{BsTyp}, Vector{BsTyp}}}

profitLossCalcMap =
  [ (Revenue,
    [OperatingRevenue, NonOperatingRevenue],
    [ExciseStaxLevy]
  )
  , (COGS,
    [CostMaterial, DirectExpenses],
    PlTyp[]
  )
  , (GrossProfit,
    [Revenue],
    [COGS]
  )
  , (Pbitda,
    [GrossProfit, OtherIncome],
    [Salaries, AdministrativeExpenses, ResearchNDevelopment, OtherOverheads, OtherOperativeExpenses, OtherExpenses, ExceptionalItems]
  )
  , (Pbitx,
    [Pbitda],
    [Depreciation, AssetImpairment, LossDivestitures, Amortization]
  )
  , (Pbtx,
    [Pbitx, InterestRevenue, OtherFinancialRevenue],
    [InterestExpense, CostDebt]
  )
  , (Pbt,
    [Pbtx],
    [ExtraordinaryItems, PriorYears]
  )
  , (Pat,
    [Pbt],
    [TaxesCurrent, TaxesDeferred]
  )
  , (NetIncome,
    [Pat, NetIncomeDiscontinuedOps],
    PlTyp[]
  )
  , (OtherComprehensiveIncome,
    [GainsLossesForex, GainsLossesActurial, GainsLossesSales, FvChangeAvlSale],
    [OtherDeferredTaxes]
  )
  , (TotalComprehensiveIncome,
    [NetIncome, OtherComprehensiveIncome],
    PlTyp[]
  )
  ]  :: Vector{Tuple{PlTyp, Vector{PlTyp}, Vector{PlTyp}}}

# cashFlowCalcMap =

const BsDict = Dict{BsTyp, Float64}
const PlDict = Dict{PlTyp, Float64}
const CfDict = Dict{CfTyp, Float64}

using Dates: Date

struct Account
  dateBegin         :: Date
  dateEnd           :: Date

  balanceSheetBegin :: Union{BsDict, Nothing}
  balanceSheetEnd   :: Union{BsDict, Nothing}
  profitLoss        :: Union{PlDict, Nothing}
  cashFlow          :: Union{CfDict, Nothing}

  others            :: Union{Dict{Symbol, Float64}, Nothing}
end

struct Params
  U   :: Float64  # Unlevered
  S   :: Float64  # Tax shield
  E   :: Float64  # Equity
  D   :: Float64  # Debt
  V   :: Float64  # Valuation

end

struct Company
  code          :: Symbol
  affiliation   :: Dict{Symbol, Float64}
  consolidated  :: Bool
  period        :: Vector{Date}
  balanceSheet  :: Dict{Date, BsDict}
  profitLoss    :: Dict{Tuple{Date, Date}, PlDict}
  cashFlow      :: Dict{Tuple{Date, Date}, CfDict}
  others        :: Dict{Tuple{Date, Date}, Dict{Symbol, Float64}}
  sharePrice    :: Union{Fl.DateSeries, Nothing}
  rate          :: Union{Vector{Params}, Nothing}
  beta          :: Union{Vector{Params}, Nothing}
end

function createCalcSets(bM :: Vector{Tuple{T, Vector{T}, Vector{T}}}) where T
  bCl = [x for (x, _, _) ∈ bM] :: Vector{T}
  bE = mapreduce( ((x,y,z),) -> (y...,z...),
      (z,y) -> [z...,filter(x -> x ∉ bCl , y)...] :: Vector{T}, bM, init = T[])
  Set(bCl), Set(bE)
end

balanceSheetCalcItems, balanceSheetEntries = createCalcSets(balanceSheetCalcMap)
profitLossCalcItems, profitLossEntries = createCalcSets(profitLossCalcMap)

getAccount(x::Company, d1::Date, d2::Date) = Account(
  d1, d2, get(x.balanceSheet,d1,nothing), get(x.balanceSheet,d2,nothing),
  get(x.profitLoss,(d1,d2),nothing), get(x.cashFlow,(d1,d2),nothing),
  get(x.others,(d1,d2),nothing)
)

putAccount!(x::Company, z::Account) = error("TODO: Add Implementation")

putDict!(x::Company, z::BsDict, d1) = error("TODO: Add Implementation")

putDict!(x::Company, z::PlDict, d1, d2) = error("TODO: Add Implementation")

putDict!(x::Company, z::CfDict, d1, d2) = error("TODO: Add Implementation")

cleanDict!(x) = error("TODO: Add Implementation")

cleanDict!(x::Company) = error("TODO: Add Implementation")

calcDict!(x::BsDict) = error("TODO: Add Implementation")

calcDict!(x::PlDict) = error("TODO: Add Implementation")

calcDict!(x::CfDict, b1, b2) = error("TODO: Add Implementation")

calcDict!(x::Company) = error("TODO: Add Implementation")

saveParquet(x) = error("TODO: Add Implementation")

readParquet(s) = error("TODO: Add Implementation")

eps(x) = error("TODO: Add Implementation")
dilutedEPS(;earn, prefDiv = 0.0, shares, sharePrice, options) = error("TODO: Add Implementation")

end

ext: jl


text:
#=
The Computer Language Benchmarks Game
 https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
 direct transliteration of the swift#3 program by Ralph Ganszky and Daniel Muellenborn:
 https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/mandelbrot-swift-3.html
 modified for Julia 1.0 by Simon Danisch
=#
const zerov8 = ntuple(x-> 0f0, 8)

# Calculate mandelbrot set for one Vec8 into one byte
Base.@propagate_inbounds function mand8(cr, ci)
    Zr = zerov8
    Zi = zerov8
    Tr = zerov8
    Ti = zerov8
    t = zerov8
    for i in 0:49
        Zi = 2f0 .* Zr .* Zi .+ ci
        Zr = Tr .- Ti .+ cr
        Tr = Zr .* Zr
        Ti = Zi .* Zi
        t = Tr .+ Ti
        all(x-> x > 4f0, t) && break
    end
    byte = UInt8(0)
    t[1] <= 4f0 && (byte |= 0x80)
    t[2] <= 4f0 && (byte |= 0x40)
    t[3] <= 4f0 && (byte |= 0x20)
    t[4] <= 4f0 && (byte |= 0x10)
    t[5] <= 4f0 && (byte |= 0x08)
    t[6] <= 4f0 && (byte |= 0x04)
    t[7] <= 4f0 && (byte |= 0x02)
    t[8] <= 4f0 && (byte |= 0x01)
    return byte
end

function mandel_inner(rows, ci, y, N, xvals)
    @simd for x in 1:8:N
        @inbounds begin
            cr = ntuple(i-> xvals[x + (i - 1)], 8)
            rows[((y-1)*N÷8+(x-1)÷8) + 1] = mand8(cr, ci)
        end
    end
end

function mandelbrot(n = 200)
    inv_ = 2.0 / n
    N = n
    xvals = zeros(Float32, n)
    yvals = zeros(Float32, n)
    Threads.@threads for i in 0:(N-1)
        @inbounds xvals[i + 1] = i * inv_ - 1.5
        @inbounds yvals[i + 1] = i * inv_ - 1.0
    end
    rows = zeros(UInt8, n*N÷8)
    Threads.@threads for y in 1:N
        @inbounds ci = yvals[y]
        mandel_inner(rows, ci, y, N, xvals)
    end
    write(stdout, "P4\n$n $n\n")
    write(stdout, rows)
end

mandelbrot(parse(Int, ARGS[1]))


text:
module PSLQ
using LinearAlgebra

export pslq

include("pslq.jl")
include("iterations.jl")

end # PSLQ


text:
#=
prompt.jl
prompt user to hit a key to continue
=#

export prompt

import MIRTjim # prompt

@deprecate prompt(args... ; kwargs...) MIRTjim.prompt(args... ; kwargs...)


text:
#=
Use these methods to convert .ts files to .csv
"""
    `TSdataset(dataset::Array)`
TSdataset takes subset of the `univariate_datasets` avilabe on the `timeseriesclassification`
website and adds `csv` files to the `data` folder after conversion (`ts` -> `csv`).
eg.
```julia
   TSdataset(["ACSF1", "Adiac"])
```
"""
function TSdataset(dataset::Array) #check on win
    for Dir in dataset
        exdir = string("data/", Dir)
        _link = string("http://timeseriesclassification.com/Downloads/", Dir ,".zip")
        Base.download(_link, Dir)
        fileFullPath = isabspath(Dir) ?  Dir : joinpath(pwd(), Dir)
        basePath = dirname(fileFullPath)
        outPath = (exdir == "" ? basePath : (isabspath(exdir) ? exdir : joinpath(pwd(),exdir)))
        isdir(outPath) ? "" : mkdir(outPath)
        zarchive = ZipFile.Reader(fileFullPath)
        for f in zarchive.files
            if f.name[end-1:end] == "ts"
                f.name = string(f.name[1:end-2], "csv")
                fullFilePath = joinpath(outPath,f.name)
                write(fullFilePath, read(f))
                open(fullFilePath) do input
                    readuntil(input, "@data")
                    write(fullFilePath, read(input))
                end
                run(`sed -i'.original' 's/:/,/g' $fullFilePath`)
                fullFilePath1 = string(fullFilePath, ".original")
                run(`rm $fullFilePath1`)
            end
        end
        close(zarchive)
        run(`rm $fileFullPath`)
    end
end

function TSdataset(filepath::String)
    exdir = string("data/", basename(filepath))
    fileFullPath = isabspath(filepath) ?  filepath : joinpath(pwd(), filepath)
    basePath = dirname(fileFullPath)
    outPath = (exdir == "" ? basePath : (isabspath(exdir) ? exdir : joinpath(pwd(),exdir)))
    isdir(outPath) ? "" : mkdir(outPath)
    files = readdir(fileFullPath, join=true)
    for f in files
        if f[end-1:end] == "ts"
            fullFilePath = joinpath(outPath, string(f[1:end-2], "csv"))
            write(fullFilePath, read(f))
            open(fullFilePath) do input
                readuntil(input, "@data")
                write(fullFilePath, read(input))
            end
            run(`sed -i'.original' 's/:/,/g' $fullFilePath`)
            fullFilePath1 = string(fullFilePath, ".original")
            run(`rm $fullFilePath1`)
        end
    end
end
=#


text:
# test module wrapping

# addprocs(1) # TODO reinsert once Julia 0.7 tests pass

using Test
using IncrementalInference

module First

export solve, AnotherInferenceType, Container

abstract type AnotherInferenceType end

mutable struct Container
  col::Dict{Symbol, Function}
  Container(::Nothing) = new()
  Container(;col=Dict{Symbol, Function}()) = new(col)
end

function registerCallback!(col::Container, fnc::Function)
  # get module
  # Symbol(string(m))
  m = Symbol(typeof(fnc).name.module)
  col.col[m] = fnc
end

function solve(ctl::Container, val::AnotherInferenceType) # m::Symbol,
  m = Symbol(typeof(val).name.module)
  if false
    # doesnt work
    evalPotential(val)
  else
    evalPotential = ctl.col[m]
    evalPotential(val)
  end
end

end

module SecondModule

using Main.First
export SecondType, SecondAgain, evalPotential, solve, registerCallback!, Container

struct SecondType <: Main.First.AnotherInferenceType
  x::Int
end
struct SecondAgain <: Main.First.AnotherInferenceType
  x::Int
end

function evalPotential(x::SecondType)
  println("evalPotential sees $(x)")
  return x.x
end

function evalPotential(x::SecondAgain)
  println("evalPotential also sees $(x)")
  return x.x
end

end

# get module
# typeof(f).name.mt.name

# Okay lets see
using Main.SecondModule

@testset "out of module evalPotential..." begin

Col = First.Container()

# running on first process only
First.registerCallback!(Col, SecondModule.evalPotential)

CCol = deepcopy(Col)

stA = 1.0
saA = 3.0

st = SecondModule.SecondType(stA)
sa = SecondModule.SecondAgain(saA)

# Symbol(typeof(st).name.module)
@test First.solve(CCol, st) == stA
@test First.solve(CCol, sa) == saA

@elapsed SecondModule.evalPotential(st)
t1 = @elapsed SecondModule.evalPotential(st)
t2 = @elapsed solve(Col, sa)

println("Check the speed is reasonable")
@test t2 < 15.0*t1 # should actually be about equal, slack for threading uncertainty

end

# expand tests to include multiprocessor
# println("Tesing call of function on separate process...")
# fr = remotecall(evalPotential,procs()[2], st)
# @test fetch(fr) == stA
#
# fr = remotecall(solve,procs()[2], CCol,sa)
# @test fetch(fr) == saA
#
#
#
# println("Stopping all but first process...")
# rmprocs(procs()[2:end])
# @show procs()

#


text:
using ArgParse
using Dates
using Base.Filesystem

function parse_commandline_extract_pyrohelio()
     s = ArgParseSettings( description = "Make pyrohelio.csv from L0 FITS files for one day.")
     @add_arg_table! s begin
         "date"
            help = "Date YYYYMMDD to extract pyroheliometer data for."
            arg_type = String
            required = true
         "--user"
            help = "NExScI Archive username"
            arg_type = String
         "--password"
            help = "NExScI Archive password"
            arg_type = String
         "--output"
            help = "Path for outputs: pyrohelio.csv"
            arg_type = String
            default = joinpath(pwd(),"pyrohelio.csv")
         "--work"
            help = "Working directory for NExScI query and L0 downloads."
            arg_type = String
         "--cookie"
            help = "Filename for NExScI cookie."
            arg_type = String
            default = "neidadmincookie.txt"
         "--query"
            help = "Filename for L0 query results."
            arg_type = String
            default = "meta_l0.csv"
         "--pyrohelio"
            help = "Path for pyroheliometer inputs: *.tel"
            arg_type = String
            #default = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/pyrohelio/202110_fromChad/"
         "--verbose"
            help = "Verbose outputs"
            action = :store_true
      end

  return parse_args(s)
end
args = parse_commandline_extract_pyrohelio()
m = match(r"(\d{4})(\d{2})(\d{2})", args["date"])
year  = parse(Int64,m[1])
month = parse(Int64,m[2])
day   = parse(Int64,m[3])
@assert 2020 <= year <= 2040
@assert 1 <= month <= 12
@assert 1 <= day <= 31

if haskey(args,"output")     output_fn          = args["output"]      end

tmp_path = !isnothing(args["work"]) ? args["work"] : Filesystem.mktempdir()
cookie_fn          = args["cookie"]
query_fn           = args["query"]
pyrohelio_dir      = args["pyrohelio"]
verbose = args["verbose"]

using FITSIO
using CSV, DataFrames
using Statistics
using NeidArchive

cookie_nexsci = joinpath(pwd(),cookie_fn) # "./neidadmincookie.txt"
query_result_file = joinpath(pwd(),query_fn) # "./criteria1.csv"
println("tmp_path = >",tmp_path,"<.")
#cookie_nexsci = joinpath(tmp_path,cookie_fn) # "./neidadmincookie.txt"
#query_result_file = joinpath(tmp_path,query_fn) # "./criteria1.csv"

if !(filesize(cookie_nexsci)>0)
   if !isnothing(args["user"]) && !isnothing(args["password"])
      NeidArchive.login(userid=args["user"], password=args["password"], cookiepath=cookie_nexsci)
   else
      NeidArchive.login(cookiepath=cookie_nexsci)
   end
end

param = Dict{String,String}()
 param["datalevel"] = "solarl0"
 param["piname"] = "Mahadevan"
 param["object"] = "Sun"
 param["datetime"] = NeidArchive.datetime_one_day_solar(Date(year,month,day))

NeidArchive.query(param, cookiepath=cookie_nexsci, outpath=query_result_file)# , outdir=tmp_path)
num_lines = countlines(query_result_file) - 1
println("# Query resulted in file with ", num_lines, " entries.")

start_dir = pwd()
#cd(tmp_path)
#NeidArchive.download(query_result_file, param["datalevel"], cookiepath=cookie_nexsci, start_row=1, end_row=3)
NeidArchive.download(query_result_file, "l0", cookiepath=cookie_nexsci, outdir=tmp_path) # , start_row=1, end_row=3)

touch(output_fn)
exampleFileIOStream =  open(output_fn,"w")
#cd(start_dir)
close(exampleFileIOStream)
exit(0)

fnin = ""
for (i,line) in enumerate(lines)
   m = match(r"(neidL0_\d{4}\d{2}\d{2}T\d+\.fits)", line)
   if (m==nothing) continue end
   fnin = string(m[1])
   fninfull = joinpath(tmp_path,fnin)
   if !(filesize(fninfull) > 0)
      m = match(r"neidL0_(\d{4})(\d{2})(\d{2})T(\d+)\.fits", fnin)
      if (m==nothing) continue end
      println("# Downloading $fnin.\n")
      try
        yyyy = m[1]
        mm = m[2]
        dd = m[3]
        tttttt = m[4]
        run(`wget --load-cookies neid_cookie.txt -O $fninfull "https://neid.ipac.caltech.edu/get_file.php?filehand=raw/$yyyy/$mm/$dd/$fnin&solar"`)
      catch ex
         println("# Failed to download $fninfull.\n")
         continue
      end
    end
    fnout = replace(fnin,"neidL0_"=>"neid_pyro_")
    fnout = replace(fnout,"fits"=>"csv")
    df = DataFrame
    if  filesize(fnout) >0
        df = CSV.read(fnout,DataFrame)
    else
    try
       f = FITS(fninfull)
       t = read(f[21],"Time")
       v = read(f[21],"Voltage")
       fd = read(f[21],"FluxDensity")
       df = DataFrame(:Time=>t, :Voltage=>v, :FluxDensity=>fd)
       CSV.write(fnout, df)
     catch ex
        println("# Failed to write $fnout.\n")
         continue
     end
     end
     try
       mean_flux = mean(df.FluxDensity)
       rms_flux = sqrt(var(df.FluxDensity, corrected=false))
       outputline = fnin * "," * string(mean_flux) * "," * string(rms_flux) * "\n"
       write(exampleFileIOStream, outputline)
       if mod(i,10) == 0 flush(exampleFileIOStream)  end
     catch ex
        println("# Failed to write summary data for $fn_out\n.")
        continue
     end
     rm(fninfull)
end
close(exampleFileIOStream)


text:
using Photon

λ = 1.55;

Lx =  λ;
Ly =  λ;
Lz =  8λ;
sim = Simulation(λ₀ = λ ;  grid = Grid(extent = (Lx,Ly,Lz) , spacing =  (λ/15,λ/15,λ/15) ))
sim(setpml!,ẑ,λ)
sim(setsrc!, PlaneWave(k̂ = (0,0,2pi/λ), ê = (0,1,0)))
sim(setTFSF! , ẑ, 1.1λ)
material = Material(ϵᵣ = 13.491 + 0.036730im )
inSphere(x) = sqrt((x[1]-Lx/2)^2 + (x[2]-Ly/2)^2+(x[3]-Lz/2)^2)  <=  1.2
sim(setmaterial!,material,inSphere)
sim(solve!)


text:
type VehicleSystem
    H::Matrix{Float64} # observation Jacobian
    R::MvNormal # process noise
    Q::MvNormal # observation noise
    Δt::Float64
    n_integration_steps::Int
    control_noise_accel::Float64
    control_noise_turnrate::Float64

    function VehicleSystem(;
        process_noise::Float64 = 0.077,
        observation_noise::Float64 = 16.7,
        control_noise_accel::Float64 = 16.7,
        control_noise_turnrate::Float64 = 0.46,
        )

        Δt = 0.1
        H = [1.0 0.0 0.0 0.0;
             0.0 1.0 0.0 0.0]
        r = process_noise
        R = MvNormal(diagm([r*0.01, r*0.01, r*0.00001, r*0.1])) # process, TODO: tune this
        q = observation_noise
        Q = MvNormal(diagm([q, q])) # obs, TODO: tune this

        n_integration_steps = 10

        new(H, R, Q, Δt, n_integration_steps, control_noise_accel, control_noise_turnrate)
    end
end

draw_proc_noise(ν::VehicleSystem) = rand(ν.R)
draw_obs_noise(ν::VehicleSystem) = rand(ν.Q)
get_process_noise_covariance(ν::VehicleSystem) = ν.R.Σ.mat
get_observation_noise_covariance(ν::VehicleSystem) = ν.Q.Σ.mat

"""
    To derive the covariance of the additional motion noise,
    we first determine the covariance matrix of the noise in control space
    Inputs:
        - ν is the vehicle concrete type
        - u is the control [a,ω]ᵀ
"""
function get_control_noise_in_control_space(ν::VehicleSystem, u::Vector{Float64})
    [ν.control_noise_accel 0.0;
     0.0 ν.control_noise_turnrate]
end

"""
    To derive the covariance of the additional motion noise,
    we transform the covariance of noise in the control space
    by a linear approximation to the derivative of the motion function
    with respect to the motion parameters
    Inputs:
        - ν is the vehicle concrete type
        - u is the control [a,ω]ᵀ
        - x is the state estimate [x,y,θ,v]ᵀ
"""
function get_transform_control_noise_to_state_space(ν::VehicleSystem, u::Vector{Float64}, x::Vector{Float64})

    x, y, θ, v = x[1], x[2], x[3], x[4]
    a, ω = u[1], u[2]
    ω² = ω*ω
    Δt = ν.Δt

    ϕ = θ + ω*Δt

    if abs(ω) < 1e-6
        [0.0                  0.0;
         0.0                  0.0;
         0.0                   Δt;
          Δt                  0.0]
    else
        [0.0   v/ω²*(sin(θ) - sin(ϕ))+v/ω*cos(ϕ)*Δt;
         0.0  -v/ω²*(cos(θ) - cos(ϕ))+v/ω*sin(ϕ)*Δt;
         0.0                   Δt;
          Δt                  0.0]
    end
end

"""
    Vehicle dynamics, return the new state
    Inputs:
        - ν is the vehicle concrete type
        - x is the state estimate [x,y,θ,v]ᵀ
        - u is the control [a,ω]ᵀ
"""
function step(ν::VehicleSystem, x::Vector{Float64}, u::Vector{Float64})

    a, ω = u[1], u[2]
    x, y, θ, v = x[1], x[2], x[3], x[4]

    δt = ν.Δt/ν.n_integration_steps

    for i in 1 : ν.n_integration_steps

        if abs(ω) < 1e-6 # simulate straight
            x += v*cos(θ)*δt
            y += v*sin(θ)*δt
        else # simulate with an arc
            x += (v/ω)*(sin(θ + ω*δt) - sin(θ))
            y += (v/ω)*(cos(θ) - cos(θ + ω*δt))
        end

        θ += ω*δt
        v += a*δt
    end

    [x,y,θ,v]
end

"""
    Vehicle observation, returns a saturated observation
    Inputs:
        - ν is the vehicle concrete type
        - x is the state estimate [x,y,θ,v]ᵀ
"""
observe(ν::VehicleSystem, x::Vector{Float64}) = [x[1], x[2]]

"""
    Computes the observation Jacobian (H matrix)
    Inputs:
        - ν is the vehicle concrete type
        - x is the state estimate [x,y,θ,v]ᵀ
"""
compute_observation_jacobian(ν::VehicleSystem, x::Vector{Float64}) = ν.H

"""
    Computes the dynamics Jacobian
    Inputs:
        - ν is the vehicle
        - x is the vehicle state [x,y,θ,v]ᵀ
        - u is the control [a,ω]
"""
function compute_dynamics_jacobian(ν::VehicleSystem, x::Vector{Float64}, u::Vector{Float64})
    Δt = ν.Δt
    θ, v = x[3], x[4]
    ω = u[2]

    if abs(ω) < 1e-6
        # drive straight
        [1.0 0.0  -v*sin(θ)*Δt  cos(θ)*Δt;
         0.0 1.0   v*cos(θ)*Δt  sin(θ)*Δt;
         0.0 0.0        1.0         0.0;
         0.0 0.0        0.0         1.0]
    else
        # drive in an arc
        ϕ = θ+ω*Δt
        [1.0 0.0  v/ω*(-cos(θ) + cos(ϕ))   -1/ω*sin(θ) + 1/ω*sin(ϕ);
         0.0 1.0  v/ω*(-sin(θ) + sin(ϕ))    1/ω*cos(θ) - 1/ω*cos(ϕ);
         0.0 0.0              1.0                     0.0;
         0.0 0.0              0.0                     1.0]
    end
end

function EKF(
    ν::VehicleSystem,
    μ::Vector{Float64}, # mean of belief at time t-1
    Σ::Matrix{Float64}, # cov of belief at time t-1
    u::Vector{Float64}, # next applied control
    z::Vector{Float64}, # observation for time t
    )

    G = compute_dynamics_jacobian(ν, μ, u)
    μbar = step(ν, μ, u)
    R = get_process_noise_covariance(ν)
    M = get_control_noise_in_control_space(ν, u)
    V = get_transform_control_noise_to_state_space(ν, u, μbar)
    # R = V*M*V'
    Σbar = G*Σ*G' + R + V*M*V'
    H = compute_observation_jacobian(ν, μbar)
    K = Σbar * H' / (H*Σbar*H' + get_observation_noise_covariance(ν))
    μ_next = μbar + K*(z - observe(ν, μbar))
    Σ_next = Σbar - K*H*Σbar
    (μ_next, Σ_next)
end

type SimulationResults
    x_arr::Matrix{Float64}
    z_arr::Matrix{Float64}
    u_arr::Matrix{Float64}
    μ_arr::Matrix{Float64}
    Σ_arr::Array{Float64, 3}
end
function simulate(ν::VehicleSystem, nsteps::Int64, x₀::Vector{Float64})
    x_arr = fill(NaN, 4, nsteps+1)
    z_arr = fill(NaN, 2, nsteps)
    u_arr = fill(NaN, 2, nsteps)
    μ_arr = fill(NaN, 4, nsteps+1)
    Σ_arr = fill(NaN, 4, 4, nsteps+1)

    # initial belief
    μ = copy(x₀)
    Σ = copy(ν.R.Σ.mat)

    x_arr[:, 1] = x₀
    μ_arr[:, 1] = μ
    Σ_arr[:, :, 1] = Σ

    x = x₀
    for i in 1 : nsteps

        # move system forward and make observation
        u = [sin(i*0.01), cos(i*0.01)+0.01]
        xₚ = step(ν, x, u) + draw_proc_noise(ν)
        z = observe(ν, xₚ) + draw_obs_noise(ν)

        # record trajectories
        x_arr[:,i+1] = xₚ
        z_arr[:,i] = z
        u_arr[:,i] = u

        # apply Kalman filter
        μ_next, Σ_next = EKF(ν, μ, Σ, u, z)
        μ_arr[:,i+1] = μ_next
        Σ_arr[:,:,i+1] = Σ_next

        copy!(x, xₚ)
        copy!(μ, μ_next)
        copy!(Σ, Σ_next)
    end

    SimulationResults(x_arr, z_arr, u_arr, μ_arr, Σ_arr)
end


text:
using Plots
using Test
using LinearAlgebra

@testset "Test everything" begin
    include("test_maxwell_sovers.jl")
    include("test_vlasov_solvers.jl")
    include("test_boltzmann_solvers.jl")
end


text:
# Based on England data (CHESS and NHS England)

# I want a way to keep this as the "average" disease progression, but modify it such that old people have less favorable outcomes (as observed)
# But correspondingly I want people at lower risk to have more favorable outcome on average

# For calculations see data_cleaning_py.ipynb, calculations from NHS England dataset as per 05 Apr
relativeDeathRisk_given_COVID_by_age = [-0.99742186, -0.99728639, -0.98158438, -0.9830432 , -0.82983414,
       -0.84039294,  0.10768979,  0.38432409,  5.13754904]

#ageRelativeDiseaseSeverity = np.array([-0.8, -0.6, -0.3, -0.3, -0.1, 0.1, 0.35, 0.4, 0.5]) # FIXED (above) - this is a guess, find data and fix
#ageRelativeRecoverySpeed = np.array([0.2]*5+[-0.1, -0.2, -0.3, -0.5]) # TODO - this is a guess, find data and fix
ageRelativeRecoverySpeed = [0.0 for i=1:9]# For now we make it same for everyone, makes calculations easier

# For calculations see data_cleaning_py.ipynb, calculations from NHS England dataset as per 05 Apr
caseFatalityRatioHospital_given_COVID_by_age = ([0.00856164, 0.03768844, 0.02321319, 0.04282494, 0.07512237,
       0.12550367, 0.167096  , 0.37953452, 0.45757006])

function trFunc_diseaseProgression(
    # Basic parameters to adhere to
    nonsymptomatic_ratio = 0.86,

    # number of days between measurable events
    infect_to_symptoms = 5.,
    #symptom_to_death = 16.,
    symptom_to_recovery = 10., # 20.5, #unrealiticly long for old people
    symptom_to_hospitalisation = 5.76,
    hospitalisation_to_recovery = 14.51,
    IgG_formation = 15.,

    # Age related parameters
    # for now we'll assume that all hospitalised cases are known (overall 23% of hospitalised COVID patients die. 9% overall case fatality ratio)
    caseFatalityRatioHospital_given_COVID_by_age = caseFatalityRatioHospital_given_COVID_by_age,
    ageRelativeRecoverySpeed = ageRelativeRecoverySpeed,

    # Unknown rates to estimate
    nonsymp_to_recovery = 15.,
    inverse_IS1_IS2 = 4.,

    )
# HACK

    nonsymptomatic_ratio = 0.86;

    # number of days between measurable events
    infect_to_symptoms = 5.;
    #symptom_to_death = 16.;
    symptom_to_recovery = 10.; # 20.5; #unrealiticly long for old people
    symptom_to_hospitalisation = 5.76;
    hospitalisation_to_recovery = 14.51;
    IgG_formation = 15.;

    # Age related parameters
    # for now we'll assume that all hospitalised cases are known (overall 23% of hospitalised COVID patients die. 9% overall case fatality ratio)
    caseFatalityRatioHospital_given_COVID_by_age = caseFatalityRatioHospital_given_COVID_by_age;
    ageRelativeRecoverySpeed = ageRelativeRecoverySpeed;

    # Unknown rates to estimate
    nonsymp_to_recovery = 15.;
    inverse_IS1_IS2 = 4.;
    # Now we have all the information to build the age-aware multistage SIR model transition matrix
    # The full transition tensor is a sparse map from the Age x HealthState x isolation state to HealthState,
        # and thus is a 4th order tensor itself, representing a linear mapping
        # from "number of people aged A in health state B and isolation state C to health state D.
    trTensor_diseaseProgression = zeros((nAge, nHS, nIso, nHS))

    # Use basic parameters to regularise inputs
    E_IS1 = 1.0/infect_to_symptoms
    # Numbers nonsymptomatic is assumed to be 86% -> E->IN / E-IS1 = 0.86/0.14
    E_IN = 0.86/0.14 * E_IS1

    # Nonsymptomatic recovery
    IN_R1 = 1.0/nonsymp_to_recovery

    IS1_IS2  = 1.0/inverse_IS1_IS2

    IS2_R1 = 1.0/(symptom_to_recovery-inverse_IS1_IS2)

    R1_R2 = 1.0/IgG_formation

    # Disease progression matrix # TODO - calibrate (together with transmissionInfectionStage)
    # rows: from-state, cols: to-state (non-symmetric!)
    # - this represent excess deaths only, doesn't contain baseline deaths!

    # Calculate all non-serious cases that do not end up in hospitals.
    # Note that we only have reliable death data from hospitals (NHS England), so we do not model people dieing outside hospitals
    diseaseProgBaseline = [
    # to: E,   IN,   IS1,   IS2,    R1,   R2,   D
          0.0  E_IN E_IS1    0   0     0   0    # from E
          0   0     0   0    IN_R1   0   0    # from IN
          0   0     0 IS1_IS2  0     0    0  # from IS1
          0   0     0    0  IS2_R1   0   0   # from IS2
          0   0     0    0    0    R1_R2  0    # from R1
          0   0     0    0    0     0   0    # from R2
          0   0     0    0    0     0   0    # from D
    ]
# TODO can be improved
#   ageAdjusted_diseaseProgBaseline = copy.deepcopy(np.repeat(diseaseProgBaseline[np.newaxis],nAge,axis=0))

ageAdjusted_diseaseProgBaseline=cat([ reshape(diseaseProgBaseline,(1,nHS-1,nHS-1)) for i=1:nAge ]...,dims=1)
ageAdjusted_diseaseProgBaseline=reshape(ageAdjusted_diseaseProgBaseline,(nAge,nHS-1,1,nHS-1))
for i=1:nIso
    trTensor_diseaseProgression[:,2:8,i,2:8]=ageAdjusted_diseaseProgBaseline
end

# # Modify all death and R1 rates:
    # for ii in range(ageAdjusted_diseaseProgBaseline.shape[1]):
    #     # Adjust death rate by age dependent disease severity
    #     ageAdjusted_diseaseProgBaseline[:,ii,end-1] = adjustRatesByAge_KeepAverageRate(
    #         ageAdjusted_diseaseProgBaseline[0,ii,end-1],
    #         ageRelativeAdjustment=relativeDeathRisk_given_COVID_by_age
    #     )
    #
    #     # Adjust recovery rate by age dependent recovery speed
    #     ageAdjusted_diseaseProgBaseline[:,ii,-3] = adjustRatesByAge_KeepAverageRate(
    #         ageAdjusted_diseaseProgBaseline[0,ii,-3],
    #         ageRelativeAdjustment=ageRelativeRecoverySpeed,
    #         agePopulationRatio=agePopulationRatio
    #     )
    #
    # ageAdjusted_diseaseProgBaseline_Hospital = copy.deepcopy(ageAdjusted_diseaseProgBaseline)
    # # Calculate hospitalisation based rates, for which we do have data. Hospitalisation can end up with deaths

    # Make sure that the ratio of recoveries in hospital honour the case fatality ratio appropriately
    # IS2 -> death
    # ageAdjusted_diseaseProgBaseline_Hospital[:,3,-1] = (
    #     # IS2 -> recovery
    #     ageAdjusted_diseaseProgBaseline_Hospital[:,3,-3] * (
    #         # multiply by cfr / (1-cfr) to get correct rate towards death
    #         caseFatalityRatioHospital_given_COVID_by_age/(
    #              1 -  caseFatalityRatioHospital_given_COVID_by_age)
    #     )
    # )

    # # TODO - time to death might be incorrect overall without an extra delay state, especially for young people
    #
    # # Non-hospitalised disease progression
    # for i1 in [0,1,3]:
    #     trTensor_diseaseProgression[:,1:,i1,1:] = ageAdjusted_diseaseProgBaseline
    #
    # # hospitalised disease progression
    # trTensor_diseaseProgression[:,1:,2,1:] = ageAdjusted_diseaseProgBaseline_Hospital

    return trTensor_diseaseProgression
end

max_stars_repo_name: vollmersj/coexist-julia

max_forks_repo_name: vollmersj/coexist-julia


text:
@testset "Algebraic datatypes" begin

@test Piecewise.supertype_all(Real[1, 2.0]) == Float64
@test Int * Float64 == Tuple{Int, Float64}
@test Int | Float64 == Union{Int, Float64}

test_union_typevar_1(::Type{T}, ::Type{Float64}) where T = T | Float64
test_union_typevar_2(::Type{Float64}, ::Type{T}) where T = Float64 | T
test_union_typevar_3(::Type{T}, ::Type{S}) where {T, S} = T | S

@test test_union_typevar_1(Int, Float64) == Int | Float64
@test test_union_typevar_2(Float64, Int) == Float64 | Int
@test test_union_typevar_3(Float64, Int) == Float64 | Int

@test (:Num1 => Int) * (:Num2 => Float64) == NamedTuple{(:Num1, :Num2), Int * Float64}
p = (:Num1 => Int) * (:Num2 => Float64)
testp = p((10, 1.0))
@test testp.Num1 === 10
@test testp.Num2 === 1.0

end


text:
struct SDivHandler end

gethandler(::Bool, ::Type{typeof(/)}, ::Type{<:SBitstreamLike}, ::Type{<:SBitstreamLike}) = SDivHandler()
init_state(::SDivHandler) = (id = 0,)

function (handler::SDivHandler)(buffer, netlist, state, inputs, outputs)
    # update netlist with inputs
    setsigned!(netlist, inputs[1], true)
    setsigned!(netlist, inputs[2], true)

    # compute output size
    lname, rname = handle_broadcast_name(name(inputs[1]), name(inputs[2]),
                                         netsize(inputs[1]), netsize(inputs[2]))
    outsize = netsize(outputs[1])

    # update netlist with output
    setsigned!(netlist, outputs[1], true)

    # add internal nets to netlist
    push!(netlist, Net(name = "div$(state.id)_out_pp", size = outsize))
    push!(netlist, Net(name = "div$(state.id)_out_mp", size = outsize))

    write(buffer, """
        $stdcomment
        // BEGIN div$(state.id)
        stoch_div_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(state.id)_pp (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_p"))),
                .B($(rname("_p"))),
                .Y(div$(state.id)_out_pp)
            );
        stoch_div_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(state.id)_mp (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_m"))),
                .B($(rname("_p"))),
                .Y(div$(state.id)_out_mp)
            );
        stoch_sat_sub_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(state.id)_p (
                .CLK(CLK),
                .nRST(nRST),
                .A(div$(state.id)_out_pp),
                .B(div$(state.id)_out_mp),
                .Y($(name(outputs[1]))_p)
            );
        stoch_sat_sub_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(state.id)_m (
                .CLK(CLK),
                .nRST(nRST),
                .A(div$(state.id)_out_mp),
                .B(div$(state.id)_out_pp),
                .Y($(name(outputs[1]))_m)
            );
        // END div$(state.id)
        \n""")

    return buffer, (id = state.id + 1,)
end


text:
export
    ZmqServer

mutable struct ZmqServer
    fg::GraphsDFG
    config::Dict{String, Any}
    isServerActive::Bool
    binding::String
end


text:
module TestEnsembleSampler

using Ensemble

using LinearAlgebra: dot, factorize, eigvals, eigvecs
using Statistics
using Test: @test, @testset

function testgaussian()
    ndim = 10

    mu = randn(ndim)

    sigma = randn(ndim, ndim)
    sigma = sigma*transpose(sigma)

    sigmafact = factorize(sigma)

    function lnprob(x)
        y = x - mu
        -0.5*dot(y, sigmafact\y)
    end

    ps = randn(ndim, 100)
    lnps = EnsembleSampler.lnprobs(ps, lnprob)
    # Burnin 1000 steps
    for i in 1:1000
        ps, lnps = EnsembleSampler.update(ps, lnps, lnprob)
    end

    # Run 1000 steps
    chain, chainlnp = EnsembleSampler.run_mcmc(ps, lnps, lnprob, 1000, thin=10)

    chain_mean = zeros(ndim)
    chain_sigma = zeros(ndim,ndim)

    for i in 1:100
        for j in 1:100
            for k in 1:ndim
                chain_mean[k] += chain[k,j,i]
            end
        end
    end
    chain_mean /= 100*100

    for i in 1:100
        for j in 1:100
            for k in 1:ndim
                x = chain[k,j,i] - chain_mean[k]
                for l in 1:ndim
                    y = chain[l,j,i] - chain_mean[l]
                    chain_sigma[k,l] += x*y
                end
            end
        end
    end
    chain_sigma /= 100*100 - 1

    @testset "Gaussian distribution tests" begin
        for i in 1:ndim
            @test isapprox(chain_mean[i], mu[i], atol=0.1*sigma[i,i])
        end

        chain_evals = eigvals(chain_sigma)
        evals = eigvals(sigma)

        for i in 1:ndim
            @test isapprox(chain_evals[i], evals[i], atol=0.1*maximum(evals))
        end

        chain_evecs = eigvecs(chain_sigma)
        evecs = eigvecs(sigma)

        for i in 1:ndim
            @test isapprox(abs(dot(chain_evecs[:,i], evecs[:,i])), 1.0, atol=0.1)
        end
    end
end

sqr(x) = x*x
sqr(x::AbstractArray{T,1}) where {T} = x.*x

# SNR = A/2*sqrt(n)/sigma

logit(x, a=0.0, b=1.0) = log(x-a) - log(b-x)

function invlogit(y::Float64, a=0.0, b=1.0)
    if y > 0
        ey = exp(-y)
        (a*ey + b)/(ey + 1.0)
    else
        ey = exp(y)
        (a + b*ey)/(1.0 + ey)
    end
end

function logitlj(y::Float64, a=0.0, b=1.0)
    if y > 0
        log(b-a) - y - 2.0*log1p(exp(-y))
    else
        log(b-a) + y - 2.0*log1p(exp(y))
    end
end

function make_sinusoidlnprob(ts::Array{Float64, 1}, data::Array{Float64, 1})
    function lnprob(x)
        lna::Float64 = x[1]
        lnP::Float64 = x[2]
        logitphi::Float64 = x[3]

        n = size(ts,1)

        phi = invlogit(logitphi, 0.0, 2*pi)

        a = exp(lna)
        P = exp(lnP)

        sig = zeros(n)
        for i in 1:n
            sig[i] = a*cos.(2.0*pi*ts[i]/P + phi)
        end

        resid = data - sig

        ll = -0.5*dot(resid,resid)
        # flat prior in a, P, phi
        lp = lna + lnP + logitlj(logitphi, 0.0, 2*pi)

        ll+lp
    end
end

function testsinusoid()
    snr = 5.0

    n = 100
    sigma = 1.0

    a = 2.0*sigma*snr/sqrt(n)
    P = 1.0

    phi0 = 2*pi*rand()

    ts = sort(100.0*rand(n))

    noise = sigma*randn(n)
    signal = a.*cos.(2.0*pi.*ts./P .+ phi0)

    data::Array{Float64, 1} = noise+signal

    lnprob = make_sinusoidlnprob(ts, data)

    ptrue = [log(a), log(P), logit(phi0, 0, 2*pi)]
    p0 = zeros(3, 100)
    for i in 1:100
        for j in 1:3
            p0[j,i] = ptrue[j] + 1e-4*randn()
        end
    end
    lnp0 = EnsembleSampler.lnprobs(p0, lnprob)

    for i in 1:1000
        p0, lnp0 = EnsembleSampler.update(p0, lnp0, lnprob)
    end

    chain, chainlnp = EnsembleSampler.run_mcmc(p0, lnp0, lnprob, 1000, thin=10)

    as = exp.(chain[1,:,:])
    Ps = exp.(chain[2,:,:])
    phis = [invlogit(c, 0, 2*pi) for c in chain[3,:,:]]

    @testset "sinusoidal probability tests" begin
        @test isapprox(mean(as), a, atol=3*std(as))
        @test isapprox(mean(Ps), P, atol=3*std(Ps))
        @test isapprox(mean(phis), phi0, atol=3*std(phis))
    end
end

function testall()
    @testset "EnsembleSampler tests" begin
        testgaussian()
        testsinusoid()
    end
end

end


text:
### definitely need to update this away from something so quick and dirty, ideally use plot recipes

### A letter is a Vector{Vector{Tuple{Float64, Float64}}}.

### Plots.Shape requires a Vector{Tuple} functions to convert this to and from a matrix
lettermatrix(L) =  [mapreduce(t -> [t[1] t[2]], vcat, t) for t in L]
lettertuple(M)  =  [map(r -> (r[1], r[2]), eachrow(T)) for T in M]

### Scale letter to unit square
function normletter(l)
    TM = lettermatrix(l)
    tmax = maximum(mapreduce(t -> maximum(t, dims=1), vcat, TM), dims=1)
    TN = [T./tmax for T in TM]
    lettertuple(TN)
end

### offset a letter in x and y
function offsetletter(L, x=0, y=0)
    M = lettermatrix(L)
    MO = [s .+ [x y] for s in M]
    lettertuple(MO)
end

### scale a letter in x and y
function scaleletter(L, x=1, y=1)
    M = lettermatrix(L)
    MO = [s .* [x y] for s in M]
    lettertuple(MO)
end

### plot a letters scaling in x and y
function plotletter!(off, scale, letter; kwargs...)
    plotletter!(offsetletter(scaleletter(letter, scale[1], scale[2]), off[1], off[2]) ; kwargs...)
end

### plot the letters
function plotletter!(letter ; bgc=:white, kwargs...)
    p = plot!(Shape(letter[1]), lab="", line=stroke(0); kwargs...)
    for l in letter[2:end]
        plot!(Shape(l), lab="", c=bgc, line=stroke(0))
    end
    p
end

### plot seq logo
function seqlogo(mot; rc=false, xo=0.0, yo=0.0, kwargs...)
    plot()
    seqlogo!(mot, rc=rc, xo=xo, yo=yo; kwargs...)
end

### plot seq logo at specified offset
function seqlogo!(mot; weight="normal", s=10, rc=false, xo=0.0, yo=0.0, label=:none, pseudo=1e-6, kwargs...)
    acgt = motif_letter_data(weight=weight, s=s)
    pwm = ifelse(rc, rcm(mot.pwm), mot.pwm) .+ pseudo
    H = log2(4) .+ sum(pwm.*log2.(pwm), dims=1)
    bases = (DNA_A, DNA_C, DNA_G, DNA_T)
    S = pwm.*H
    p = plot!(; kwargs...)
    off = xo .+ 0.5
    cc = palette(:default)[1:4]
    cc = [:green, :steelblue, :orange, :red]
    for c in eachcol(S)
        si = sortperm(c)
        yoff = yo
        for s in si
            v = c[s]
            b = bases[s]
            cl = cc[s]
           plotletter!([off, yoff], [1.0, v], acgt[b], c=cl)
           yoff += v
        end

        off += 1
    end

    if label == :left
        annotate!((xo, yo .+ 1, text(mot.name, font(8, :right, :middle))))
    elseif label == :right
        annotate!((off, yo .+ 1, text(mot.name, font(8, :left, :middle))))

    elseif label == :top
        annotate!(((off .+ xo)/2, yo + 2, text(mot.name, font(8, :top, :center))))
    end
    p
end

### for plotting sequences
function plotseq(seq; weight="normal", s=10, xo = 0.5, yo = 0, kwargs...)
    acgt=motif_letter_data(weight=weight, s=s)
    p = plot(; kwargs...)

    offx = xo
    offy = yo
    for (i, s) in enumerate(seq)
        plotletter!([offx, offy], [1.0, 0.25], acgt[s], c=:black)
        offx += 1
    end
    p
end


text:
# Julia wrapper for header: pigpio.h
# Automatically generated using Clang.jl

function gpioInitialise()
    ccall((:gpioInitialise, pigpio), Cint, ())
end

function gpioTerminate()
    ccall((:gpioTerminate, pigpio), Cvoid, ())
end

function gpioSetMode(gpio, mode)
    ccall((:gpioSetMode, pigpio), Cint, (UInt32, UInt32), gpio, mode)
end

function gpioGetMode(gpio)
    ccall((:gpioGetMode, pigpio), Cint, (UInt32,), gpio)
end

function gpioSetPullUpDown(gpio, pud)
    ccall((:gpioSetPullUpDown, pigpio), Cint, (UInt32, UInt32), gpio, pud)
end

function gpioRead(gpio)
    ccall((:gpioRead, pigpio), Cint, (UInt32,), gpio)
end

function gpioWrite(gpio, level)
    ccall((:gpioWrite, pigpio), Cint, (UInt32, UInt32), gpio, level)
end

function gpioPWM(user_gpio, dutycycle)
    ccall((:gpioPWM, pigpio), Cint, (UInt32, UInt32), user_gpio, dutycycle)
end

function gpioGetPWMdutycycle(user_gpio)
    ccall((:gpioGetPWMdutycycle, pigpio), Cint, (UInt32,), user_gpio)
end

function gpioSetPWMrange(user_gpio, range)
    ccall((:gpioSetPWMrange, pigpio), Cint, (UInt32, UInt32), user_gpio, range)
end

function gpioGetPWMrange(user_gpio)
    ccall((:gpioGetPWMrange, pigpio), Cint, (UInt32,), user_gpio)
end

function gpioGetPWMrealRange(user_gpio)
    ccall((:gpioGetPWMrealRange, pigpio), Cint, (UInt32,), user_gpio)
end

function gpioSetPWMfrequency(user_gpio, frequency)
    ccall((:gpioSetPWMfrequency, pigpio), Cint, (UInt32, UInt32), user_gpio, frequency)
end

function gpioGetPWMfrequency(user_gpio)
    ccall((:gpioGetPWMfrequency, pigpio), Cint, (UInt32,), user_gpio)
end

function gpioServo(user_gpio, pulsewidth)
    ccall((:gpioServo, pigpio), Cint, (UInt32, UInt32), user_gpio, pulsewidth)
end

function gpioGetServoPulsewidth(user_gpio)
    ccall((:gpioGetServoPulsewidth, pigpio), Cint, (UInt32,), user_gpio)
end

function gpioSetAlertFunc(user_gpio, f)
    ccall((:gpioSetAlertFunc, pigpio), Cint, (UInt32, gpioAlertFunc_t), user_gpio, f)
end

function gpioSetAlertFuncEx(user_gpio, f, userdata)
    ccall((:gpioSetAlertFuncEx, pigpio), Cint, (UInt32, gpioAlertFuncEx_t, Ptr{Cvoid}), user_gpio, f, userdata)
end

function gpioSetISRFunc(gpio, edge, timeout, f)
    ccall((:gpioSetISRFunc, pigpio), Cint, (UInt32, UInt32, Cint, gpioISRFunc_t), gpio, edge, timeout, f)
end

function gpioSetISRFuncEx(gpio, edge, timeout, f, userdata)
    ccall((:gpioSetISRFuncEx, pigpio), Cint, (UInt32, UInt32, Cint, gpioISRFuncEx_t, Ptr{Cvoid}), gpio, edge, timeout, f, userdata)
end

function gpioNotifyOpen()
    ccall((:gpioNotifyOpen, pigpio), Cint, ())
end

function gpioNotifyOpenWithSize(bufSize)
    ccall((:gpioNotifyOpenWithSize, pigpio), Cint, (Cint,), bufSize)
end

function gpioNotifyBegin(handle, bits)
    ccall((:gpioNotifyBegin, pigpio), Cint, (UInt32, UInt32), handle, bits)
end

function gpioNotifyPause(handle)
    ccall((:gpioNotifyPause, pigpio), Cint, (UInt32,), handle)
end

function gpioNotifyClose(handle)
    ccall((:gpioNotifyClose, pigpio), Cint, (UInt32,), handle)
end

function gpioWaveClear()
    ccall((:gpioWaveClear, pigpio), Cint, ())
end

function gpioWaveAddNew()
    ccall((:gpioWaveAddNew, pigpio), Cint, ())
end

function gpioWaveAddGeneric(numPulses, pulses)
    ccall((:gpioWaveAddGeneric, pigpio), Cint, (UInt32, Ptr{gpioPulse_t}), numPulses, pulses)
end

function gpioWaveAddSerial(user_gpio, baud, data_bits, stop_bits, offset, numBytes, str)
    ccall((:gpioWaveAddSerial, pigpio), Cint, (UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, Cstring), user_gpio, baud, data_bits, stop_bits, offset, numBytes, str)
end

function gpioWaveCreate()
    ccall((:gpioWaveCreate, pigpio), Cint, ())
end

function gpioWaveCreatePad(pctCB, pctBOOL, pctTOOL)
    ccall((:gpioWaveCreatePad, pigpio), Cint, (Cint, Cint, Cint), pctCB, pctBOOL, pctTOOL)
end

function gpioWaveDelete(wave_id)
    ccall((:gpioWaveDelete, pigpio), Cint, (UInt32,), wave_id)
end

function gpioWaveTxSend(wave_id, wave_mode)
    ccall((:gpioWaveTxSend, pigpio), Cint, (UInt32, UInt32), wave_id, wave_mode)
end

function gpioWaveChain(buf, bufSize)
    ccall((:gpioWaveChain, pigpio), Cint, (Cstring, UInt32), buf, bufSize)
end

function gpioWaveTxAt()
    ccall((:gpioWaveTxAt, pigpio), Cint, ())
end

function gpioWaveTxBusy()
    ccall((:gpioWaveTxBusy, pigpio), Cint, ())
end

function gpioWaveTxStop()
    ccall((:gpioWaveTxStop, pigpio), Cint, ())
end

function gpioWaveGetMicros()
    ccall((:gpioWaveGetMicros, pigpio), Cint, ())
end

function gpioWaveGetHighMicros()
    ccall((:gpioWaveGetHighMicros, pigpio), Cint, ())
end

function gpioWaveGetMaxMicros()
    ccall((:gpioWaveGetMaxMicros, pigpio), Cint, ())
end

function gpioWaveGetPulses()
    ccall((:gpioWaveGetPulses, pigpio), Cint, ())
end

function gpioWaveGetHighPulses()
    ccall((:gpioWaveGetHighPulses, pigpio), Cint, ())
end

function gpioWaveGetMaxPulses()
    ccall((:gpioWaveGetMaxPulses, pigpio), Cint, ())
end

function gpioWaveGetCbs()
    ccall((:gpioWaveGetCbs, pigpio), Cint, ())
end

function gpioWaveGetHighCbs()
    ccall((:gpioWaveGetHighCbs, pigpio), Cint, ())
end

function gpioWaveGetMaxCbs()
    ccall((:gpioWaveGetMaxCbs, pigpio), Cint, ())
end

function gpioSerialReadOpen(user_gpio, baud, data_bits)
    ccall((:gpioSerialReadOpen, pigpio), Cint, (UInt32, UInt32, UInt32), user_gpio, baud, data_bits)
end

function gpioSerialReadInvert(user_gpio, invert)
    ccall((:gpioSerialReadInvert, pigpio), Cint, (UInt32, UInt32), user_gpio, invert)
end

function gpioSerialRead(user_gpio, buf, bufSize)
    ccall((:gpioSerialRead, pigpio), Cint, (UInt32, Ptr{Cvoid}, Cint), user_gpio, buf, bufSize)
end

function gpioSerialReadClose(user_gpio)
    ccall((:gpioSerialReadClose, pigpio), Cint, (UInt32,), user_gpio)
end

function i2cOpen(i2cBus, i2cAddr, i2cFlags)
    ccall((:i2cOpen, pigpio), Cint, (UInt32, UInt32, UInt32), i2cBus, i2cAddr, i2cFlags)
end

function i2cClose(handle)
    ccall((:i2cClose, pigpio), Cint, (UInt32,), handle)
end

function i2cWriteQuick(handle, bit)
    ccall((:i2cWriteQuick, pigpio), Cint, (UInt32, UInt32), handle, bit)
end

function i2cWriteByte(handle, bVal)
    ccall((:i2cWriteByte, pigpio), Cint, (UInt32, UInt32), handle, bVal)
end

function i2cReadByte(handle)
    ccall((:i2cReadByte, pigpio), Cint, (UInt32,), handle)
end

function i2cWriteByteData(handle, i2cReg, bVal)
    ccall((:i2cWriteByteData, pigpio), Cint, (UInt32, UInt32, UInt32), handle, i2cReg, bVal)
end

function i2cWriteWordData(handle, i2cReg, wVal)
    ccall((:i2cWriteWordData, pigpio), Cint, (UInt32, UInt32, UInt32), handle, i2cReg, wVal)
end

function i2cReadByteData(handle, i2cReg)
    ccall((:i2cReadByteData, pigpio), Cint, (UInt32, UInt32), handle, i2cReg)
end

function i2cReadWordData(handle, i2cReg)
    ccall((:i2cReadWordData, pigpio), Cint, (UInt32, UInt32), handle, i2cReg)
end

function i2cProcessCall(handle, i2cReg, wVal)
    ccall((:i2cProcessCall, pigpio), Cint, (UInt32, UInt32, UInt32), handle, i2cReg, wVal)
end

function i2cWriteBlockData(handle, i2cReg, buf, count)
    ccall((:i2cWriteBlockData, pigpio), Cint, (UInt32, UInt32, Cstring, UInt32), handle, i2cReg, buf, count)
end

function i2cReadBlockData(handle, i2cReg, buf)
    ccall((:i2cReadBlockData, pigpio), Cint, (UInt32, UInt32, Cstring), handle, i2cReg, buf)
end

function i2cBlockProcessCall(handle, i2cReg, buf, count)
    ccall((:i2cBlockProcessCall, pigpio), Cint, (UInt32, UInt32, Cstring, UInt32), handle, i2cReg, buf, count)
end

function i2cReadI2CBlockData(handle, i2cReg, buf, count)
    ccall((:i2cReadI2CBlockData, pigpio), Cint, (UInt32, UInt32, Cstring, UInt32), handle, i2cReg, buf, count)
end

function i2cWriteI2CBlockData(handle, i2cReg, buf, count)
    ccall((:i2cWriteI2CBlockData, pigpio), Cint, (UInt32, UInt32, Cstring, UInt32), handle, i2cReg, buf, count)
end

function i2cReadDevice(handle, buf, count)
    ccall((:i2cReadDevice, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function i2cWriteDevice(handle, buf, count)
    ccall((:i2cWriteDevice, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function i2cSwitchCombined(setting)
    ccall((:i2cSwitchCombined, pigpio), Cvoid, (Cint,), setting)
end

function i2cSegments(handle, segs, numSegs)
    ccall((:i2cSegments, pigpio), Cint, (UInt32, Ptr{pi_i2c_msg_t}, UInt32), handle, segs, numSegs)
end

function i2cZip(handle, inBuf, inLen, outBuf, outLen)
    ccall((:i2cZip, pigpio), Cint, (UInt32, Cstring, UInt32, Cstring, UInt32), handle, inBuf, inLen, outBuf, outLen)
end

function bbI2COpen(SDA, SCL, baud)
    ccall((:bbI2COpen, pigpio), Cint, (UInt32, UInt32, UInt32), SDA, SCL, baud)
end

function bbI2CClose(SDA)
    ccall((:bbI2CClose, pigpio), Cint, (UInt32,), SDA)
end

function bbI2CZip(SDA, inBuf, inLen, outBuf, outLen)
    ccall((:bbI2CZip, pigpio), Cint, (UInt32, Cstring, UInt32, Cstring, UInt32), SDA, inBuf, inLen, outBuf, outLen)
end

function bscXfer(bsc_xfer)
    ccall((:bscXfer, pigpio), Cint, (Ptr{bsc_xfer_t},), bsc_xfer)
end

function bbSPIOpen(CS, MISO, MOSI, SCLK, baud, spiFlags)
    ccall((:bbSPIOpen, pigpio), Cint, (UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), CS, MISO, MOSI, SCLK, baud, spiFlags)
end

function bbSPIClose(CS)
    ccall((:bbSPIClose, pigpio), Cint, (UInt32,), CS)
end

function bbSPIXfer(CS, inBuf, outBuf, count)
    ccall((:bbSPIXfer, pigpio), Cint, (UInt32, Ptr{Cvoid}, Ptr{Cvoid}, UInt32), CS, inBuf, outBuf, count)
end

function spiOpen(spiChan, baud, spiFlags)
    ccall((:spiOpen, pigpio), Cint, (UInt32, UInt32, UInt32), spiChan, baud, spiFlags)
end

function spiClose(handle)
    ccall((:spiClose, pigpio), Cint, (UInt32,), handle)
end

function spiRead(handle, buf, count)
    ccall((:spiRead, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function spiWrite(handle, buf, count)
    ccall((:spiWrite, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function spiXfer(handle, txBuf, rxBuf, count)
    ccall((:spiXfer, pigpio), Cint, (UInt32, Cstring, Cstring, UInt32), handle, txBuf, rxBuf, count)
end

function serOpen(sertty, baud, serFlags)
    ccall((:serOpen, pigpio), Cint, (Cstring, UInt32, UInt32), sertty, baud, serFlags)
end

function serClose(handle)
    ccall((:serClose, pigpio), Cint, (UInt32,), handle)
end

function serWriteByte(handle, bVal)
    ccall((:serWriteByte, pigpio), Cint, (UInt32, UInt32), handle, bVal)
end

function serReadByte(handle)
    ccall((:serReadByte, pigpio), Cint, (UInt32,), handle)
end

function serWrite(handle, buf, count)
    ccall((:serWrite, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function serRead(handle, buf, count)
    ccall((:serRead, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function serDataAvailable(handle)
    ccall((:serDataAvailable, pigpio), Cint, (UInt32,), handle)
end

function gpioTrigger(user_gpio, pulseLen, level)
    ccall((:gpioTrigger, pigpio), Cint, (UInt32, UInt32, UInt32), user_gpio, pulseLen, level)
end

function gpioSetWatchdog(user_gpio, timeout)
    ccall((:gpioSetWatchdog, pigpio), Cint, (UInt32, UInt32), user_gpio, timeout)
end

function gpioNoiseFilter(user_gpio, steady, active)
    ccall((:gpioNoiseFilter, pigpio), Cint, (UInt32, UInt32, UInt32), user_gpio, steady, active)
end

function gpioGlitchFilter(user_gpio, steady)
    ccall((:gpioGlitchFilter, pigpio), Cint, (UInt32, UInt32), user_gpio, steady)
end

function gpioSetGetSamplesFunc(f, bits)
    ccall((:gpioSetGetSamplesFunc, pigpio), Cint, (gpioGetSamplesFunc_t, UInt32), f, bits)
end

function gpioSetGetSamplesFuncEx(f, bits, userdata)
    ccall((:gpioSetGetSamplesFuncEx, pigpio), Cint, (gpioGetSamplesFuncEx_t, UInt32, Ptr{Cvoid}), f, bits, userdata)
end

function gpioSetTimerFunc(timer, millis, f)
    ccall((:gpioSetTimerFunc, pigpio), Cint, (UInt32, UInt32, gpioTimerFunc_t), timer, millis, f)
end

function gpioSetTimerFuncEx(timer, millis, f, userdata)
    ccall((:gpioSetTimerFuncEx, pigpio), Cint, (UInt32, UInt32, gpioTimerFuncEx_t, Ptr{Cvoid}), timer, millis, f, userdata)
end

function gpioStoreScript(script)
    ccall((:gpioStoreScript, pigpio), Cint, (Cstring,), script)
end

function gpioRunScript(script_id, numPar, param)
    ccall((:gpioRunScript, pigpio), Cint, (UInt32, UInt32, Ptr{UInt32}), script_id, numPar, param)
end

function gpioUpdateScript(script_id, numPar, param)
    ccall((:gpioUpdateScript, pigpio), Cint, (UInt32, UInt32, Ptr{UInt32}), script_id, numPar, param)
end

function gpioScriptStatus(script_id, param)
    ccall((:gpioScriptStatus, pigpio), Cint, (UInt32, Ptr{UInt32}), script_id, param)
end

function gpioStopScript(script_id)
    ccall((:gpioStopScript, pigpio), Cint, (UInt32,), script_id)
end

function gpioDeleteScript(script_id)
    ccall((:gpioDeleteScript, pigpio), Cint, (UInt32,), script_id)
end

function gpioSetSignalFunc(signum, f)
    ccall((:gpioSetSignalFunc, pigpio), Cint, (UInt32, gpioSignalFunc_t), signum, f)
end

function gpioSetSignalFuncEx(signum, f, userdata)
    ccall((:gpioSetSignalFuncEx, pigpio), Cint, (UInt32, gpioSignalFuncEx_t, Ptr{Cvoid}), signum, f, userdata)
end

function gpioRead_Bits_0_31()
    ccall((:gpioRead_Bits_0_31, pigpio), UInt32, ())
end

function gpioRead_Bits_32_53()
    ccall((:gpioRead_Bits_32_53, pigpio), UInt32, ())
end

function gpioWrite_Bits_0_31_Clear(bits)
    ccall((:gpioWrite_Bits_0_31_Clear, pigpio), Cint, (UInt32,), bits)
end

function gpioWrite_Bits_32_53_Clear(bits)
    ccall((:gpioWrite_Bits_32_53_Clear, pigpio), Cint, (UInt32,), bits)
end

function gpioWrite_Bits_0_31_Set(bits)
    ccall((:gpioWrite_Bits_0_31_Set, pigpio), Cint, (UInt32,), bits)
end

function gpioWrite_Bits_32_53_Set(bits)
    ccall((:gpioWrite_Bits_32_53_Set, pigpio), Cint, (UInt32,), bits)
end

function gpioHardwareClock(gpio, clkfreq)
    ccall((:gpioHardwareClock, pigpio), Cint, (UInt32, UInt32), gpio, clkfreq)
end

function gpioHardwarePWM(gpio, PWMfreq, PWMduty)
    ccall((:gpioHardwarePWM, pigpio), Cint, (UInt32, UInt32, UInt32), gpio, PWMfreq, PWMduty)
end

function gpioTime(timetype, seconds, micros)
    ccall((:gpioTime, pigpio), Cint, (UInt32, Ptr{Cint}, Ptr{Cint}), timetype, seconds, micros)
end

function gpioSleep(timetype, seconds, micros)
    ccall((:gpioSleep, pigpio), Cint, (UInt32, Cint, Cint), timetype, seconds, micros)
end

function gpioDelay(micros)
    ccall((:gpioDelay, pigpio), UInt32, (UInt32,), micros)
end

function gpioTick()
    ccall((:gpioTick, pigpio), UInt32, ())
end

function gpioHardwareRevision()
    ccall((:gpioHardwareRevision, pigpio), UInt32, ())
end

function gpioVersion()
    ccall((:gpioVersion, pigpio), UInt32, ())
end

function gpioGetPad(pad)
    ccall((:gpioGetPad, pigpio), Cint, (UInt32,), pad)
end

function gpioSetPad(pad, padStrength)
    ccall((:gpioSetPad, pigpio), Cint, (UInt32, UInt32), pad, padStrength)
end

function eventMonitor(handle, bits)
    ccall((:eventMonitor, pigpio), Cint, (UInt32, UInt32), handle, bits)
end

function eventSetFunc(event, f)
    ccall((:eventSetFunc, pigpio), Cint, (UInt32, eventFunc_t), event, f)
end

function eventSetFuncEx(event, f, userdata)
    ccall((:eventSetFuncEx, pigpio), Cint, (UInt32, eventFuncEx_t, Ptr{Cvoid}), event, f, userdata)
end

function eventTrigger(event)
    ccall((:eventTrigger, pigpio), Cint, (UInt32,), event)
end

function shell(scriptName, scriptString)
    ccall((:shell, pigpio), Cint, (Cstring, Cstring), scriptName, scriptString)
end

function fileOpen(file, mode)
    ccall((:fileOpen, pigpio), Cint, (Cstring, UInt32), file, mode)
end

function fileClose(handle)
    ccall((:fileClose, pigpio), Cint, (UInt32,), handle)
end

function fileWrite(handle, buf, count)
    ccall((:fileWrite, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function fileRead(handle, buf, count)
    ccall((:fileRead, pigpio), Cint, (UInt32, Cstring, UInt32), handle, buf, count)
end

function fileSeek(handle, seekOffset, seekFrom)
    ccall((:fileSeek, pigpio), Cint, (UInt32, Int32, Cint), handle, seekOffset, seekFrom)
end

function fileList(fpat, buf, count)
    ccall((:fileList, pigpio), Cint, (Cstring, Cstring, UInt32), fpat, buf, count)
end

function gpioCfgBufferSize(cfgMillis)
    ccall((:gpioCfgBufferSize, pigpio), Cint, (UInt32,), cfgMillis)
end

function gpioCfgClock(cfgMicros, cfgPeripheral, cfgSource)
    ccall((:gpioCfgClock, pigpio), Cint, (UInt32, UInt32, UInt32), cfgMicros, cfgPeripheral, cfgSource)
end

function gpioCfgDMAchannel(DMAchannel)
    ccall((:gpioCfgDMAchannel, pigpio), Cint, (UInt32,), DMAchannel)
end

function gpioCfgDMAchannels(primaryChannel, secondaryChannel)
    ccall((:gpioCfgDMAchannels, pigpio), Cint, (UInt32, UInt32), primaryChannel, secondaryChannel)
end

function gpioCfgPermissions(updateMask)
    ccall((:gpioCfgPermissions, pigpio), Cint, (UInt64,), updateMask)
end

function gpioCfgSocketPort(port)
    ccall((:gpioCfgSocketPort, pigpio), Cint, (UInt32,), port)
end

function gpioCfgInterfaces(ifFlags)
    ccall((:gpioCfgInterfaces, pigpio), Cint, (UInt32,), ifFlags)
end

function gpioCfgMemAlloc(memAllocMode)
    ccall((:gpioCfgMemAlloc, pigpio), Cint, (UInt32,), memAllocMode)
end

function gpioCfgNetAddr(numSockAddr, sockAddr)
    ccall((:gpioCfgNetAddr, pigpio), Cint, (Cint, Ptr{UInt32}), numSockAddr, sockAddr)
end

function gpioCfgGetInternals()
    ccall((:gpioCfgGetInternals, pigpio), UInt32, ())
end

function gpioCfgSetInternals(cfgVal)
    ccall((:gpioCfgSetInternals, pigpio), Cint, (UInt32,), cfgVal)
end

function gpioCustom1(arg1, arg2, argx, argc)
    ccall((:gpioCustom1, pigpio), Cint, (UInt32, UInt32, Cstring, UInt32), arg1, arg2, argx, argc)
end

function gpioCustom2(arg1, argx, argc, retBuf, retMax)
    ccall((:gpioCustom2, pigpio), Cint, (UInt32, Cstring, UInt32, Cstring, UInt32), arg1, argx, argc, retBuf, retMax)
end

function rawWaveAddSPI(spi, offset, spiSS, buf, spiTxBits, spiBitFirst, spiBitLast, spiBits)
    ccall((:rawWaveAddSPI, pigpio), Cint, (Ptr{rawSPI_t}, UInt32, UInt32, Cstring, UInt32, UInt32, UInt32, UInt32), spi, offset, spiSS, buf, spiTxBits, spiBitFirst, spiBitLast, spiBits)
end

function rawWaveAddGeneric(numPulses, pulses)
    ccall((:rawWaveAddGeneric, pigpio), Cint, (UInt32, Ptr{rawWave_t}), numPulses, pulses)
end

function rawWaveCB()
    ccall((:rawWaveCB, pigpio), UInt32, ())
end

function rawWaveCBAdr(cbNum)
    ccall((:rawWaveCBAdr, pigpio), Ptr{rawCbs_t}, (Cint,), cbNum)
end

function rawWaveGetOOL(pos)
    ccall((:rawWaveGetOOL, pigpio), UInt32, (Cint,), pos)
end

function rawWaveSetOOL(pos, lVal)
    ccall((:rawWaveSetOOL, pigpio), Cvoid, (Cint, UInt32), pos, lVal)
end

function rawWaveGetOut(pos)
    ccall((:rawWaveGetOut, pigpio), UInt32, (Cint,), pos)
end

function rawWaveSetOut(pos, lVal)
    ccall((:rawWaveSetOut, pigpio), Cvoid, (Cint, UInt32), pos, lVal)
end

function rawWaveGetIn(pos)
    ccall((:rawWaveGetIn, pigpio), UInt32, (Cint,), pos)
end

function rawWaveSetIn(pos, lVal)
    ccall((:rawWaveSetIn, pigpio), Cvoid, (Cint, UInt32), pos, lVal)
end

function rawWaveInfo(wave_id)
    ccall((:rawWaveInfo, pigpio), rawWaveInfo_t, (Cint,), wave_id)
end

function getBitInBytes(bitPos, buf, numBits)
    ccall((:getBitInBytes, pigpio), Cint, (Cint, Cstring, Cint), bitPos, buf, numBits)
end

function putBitInBytes(bitPos, buf, bit)
    ccall((:putBitInBytes, pigpio), Cvoid, (Cint, Cstring, Cint), bitPos, buf, bit)
end

function time_time()
    ccall((:time_time, pigpio), Cdouble, ())
end

function time_sleep(seconds)
    ccall((:time_sleep, pigpio), Cvoid, (Cdouble,), seconds)
end

function rawDumpWave()
    ccall((:rawDumpWave, pigpio), Cvoid, ())
end

function rawDumpScript(script_id)
    ccall((:rawDumpScript, pigpio), Cvoid, (UInt32,), script_id)
end


text:
using Test, Random, Statistics, LinearAlgebra
import transD_GP
include("00_test_TransD_GP_s_ns.jl")
include("00_test_TransD_GP_ns.jl")
include("00_test_TransD_GP_s.jl")


text:
using Documenter, MLJSerialization

makedocs(
    modules = [MLJSerialization],
    sitename = "MLJSerialization.jl",
)

deploydocs(
    repo = "github.com/alan-turing-institute/MLJSerialization.jl.git",
)


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# This example contains a listing of the source code given in the first example
# of the README.md package documentation

# Load the TOAST module
using TOAST

# Load a mesh from disk.
meshdir = joinpath(TOAST._jl_toast_dir, "test", "2D", "meshes")
mesh = Mesh(joinpath(meshdir, "circle25_32.msh"));

# Construct nodal coefficient vectors representing the parameters
μ = fill(NodalCoeff, mesh, 0.01);
κ = fill(NodalCoeff, mesh, 0.3);
ζ = fill(NodalCoeff, mesh, 0.5);

# Build the system matrix in-place
S = SystemMatrix(mesh)
assemble!(S, PDD, κ);
assemble!(S, PFF, μ);
assemble!(S, BNDPFF, ζ);

# Construct a special-purpose source term
source = PointSource(mesh, [0.,0.], Isotropic)
Q = assemble(source)

# Convert the system matrix from its CSR form to a Julia CSC sparse matrix
sysmat = sparse(S)

# Solve as you like
ϕ = sparse(S)\Q

# If you have PyPlot installed...
# using PyPlot
# vtx, = data(mesh);
# plot_trisurf(vtx[:,1], vtx[:,2], ϕ)


text:
# Variable functions

export Variable
"""
    struct Variable

Holds a C pointer `adios2_variable *`.
"""
struct Variable
    ptr::Ptr{Cvoid}
    adios::Adios
    Variable(ptr::Ptr{Cvoid}, adios::Adios) = new(ptr, adios)
end

function Base.show(io::IO, variable::Variable)
    nm = name(variable)
    T = type(variable)
    sid = shapeid(variable)
    sh = shape(variable)
    print(io, "Variable(name=$nm,type=$T,shapeid=$sid,shape=$sh)")
    return
end

function Base.show(io::IO, ::MIME"text/plain", variable::Variable)
    nm = name(variable)
    print(io, "Variable($nm)")
    return
end

export set_block_selection
"""
    set_block_selection(variable::Variable, block_id::Int)
"""
function set_block_selection(variable::Variable, block_id::Int)
    err = ccall((:adios2_set_block_selection, libadios2_c), Cint,
                (Ptr{Cvoid}, Csize_t), variable.ptr, block_id)
    Error(err) ≠ error_none && return nothing
    return ()
end

export name
"""
    var_name = name(variable::Variable)
    var_name::Union{Nothing,String}

Retrieve variable name.
"""
function name(variable::Variable)
    size = Ref{Csize_t}()
    err = ccall((:adios2_variable_name, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), C_NULL, size,
                variable.ptr)
    Error(err) ≠ error_none && return nothing
    name = Array{Cchar}(undef, size[])
    err = ccall((:adios2_variable_name, libadios2_c), Cint,
                (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), name, size,
                variable.ptr)
    Error(err) ≠ error_none && return nothing
    return unsafe_string(pointer(name), size[])
end

export type
"""
    var_type = type(variable::Variable)
    var_type::Union{Nothing,Type}

Retrieve variable type.
"""
function type(variable::Variable)
    type = Ref{Cint}()
    err = ccall((:adios2_variable_type, libadios2_c), Cint,
                (Ref{Cint}, Ptr{Cvoid}), type, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return julia_type(AType(type[]))
end

# export type_string
# """
#     type_string = variable_type_string(variable::Variable)
#     type_string::Union{Nothing,String}
#
# Retrieve variable type in string form "char", "unsigned long", etc.
# This reports C type names.
# """
# function type_string(variable::Variable)
#     size = Ref{Csize_t}()
#     err = ccall((:adios2_variable_type_string, libadios2_c), Cint,
#                 (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), C_NULL, size,
#                 variable.ptr)
#     Error(err) ≠ error_none && return nothing
#     type = '\0'^size[]
#     err = ccall((:adios2_variable_type_string, libadios2_c), Cint,
#                 (Ptr{Cchar}, Ref{Csize_t}, Ptr{Cvoid}), name, size,
#                 variable.ptr)
#     Error(err) ≠ error_none && return nothing
#     return type
# end

export shapeid
"""
    var_shapeid = shapeid(variable::Variable)
    var_shapeid::Union{Nothing,ShapeId}

Retrieve variable shapeid.
"""
function shapeid(variable::Variable)
    shapeid = Ref{Cint}()
    err = ccall((:adios2_variable_shapeid, libadios2_c), Cint,
                (Ref{Cint}, Ptr{Cvoid}), shapeid, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return ShapeId(shapeid[])
end

"""
    var_ndims = ndims(variable::Variable)
    var_ndims::Union{Nothing,Int}

Retrieve current variable number of dimensions.
"""
function Base.ndims(variable::Variable)
    var_shapeid = shapeid(variable)
    var_shapeid ≡ nothing && return nothing
    var_shapeid == shapeid_local_value && return 0
    ndims = Ref{Csize_t}()
    err = ccall((:adios2_variable_ndims, libadios2_c), Cint,
                (Ref{Csize_t}, Ptr{Cvoid}), ndims, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Int(ndims[])
end

export shape
"""
    var_shape = shape(variable::Variable)
    var_shape::Union{Nothing,NTuple{N,Int} where N}

Retrieve current variable shape.
"""
function shape(variable::Variable)
    var_shapeid = shapeid(variable)
    var_shapeid ≡ nothing && return nothing
    var_shapeid ∈ (shapeid_local_value, shapeid_local_array) && return nothing
    D = ndims(variable)
    shape = Array{Csize_t}(undef, D)
    err = ccall((:adios2_variable_shape, libadios2_c), Cint,
                (Ptr{Csize_t}, Ptr{Cvoid}), shape, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Tuple(reverse!(shape))
end

export start
"""
    var_start = start(variable::Variable)
    var_start::Union{Nothing,NTuple{N,Int} where N,}

Retrieve current variable start.
"""
function start(variable::Variable)
    var_shapeid = shapeid(variable)
    var_shapeid ≡ nothing && return nothing
    var_shapeid ∈ (shapeid_local_value, shapeid_local_array) && return nothing
    D = ndims(variable)
    start = Array{Csize_t}(undef, D)
    err = ccall((:adios2_variable_start, libadios2_c), Cint,
                (Ptr{Csize_t}, Ptr{Cvoid}), start, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Tuple(reverse!(start))
end

"""
    var_count = count(variable::Variable)
    var_count::Union{Nothing,NTuple{N,Int} where N}

Retrieve current variable count.
"""
function Base.count(variable::Variable)
    var_shapeid = shapeid(variable)
    var_shapeid ≡ nothing && return nothing
    var_shapeid == shapeid_local_value && return ()
    D = ndims(variable)
    count = Array{Csize_t}(undef, D)
    err = ccall((:adios2_variable_count, libadios2_c), Cint,
                (Ptr{Csize_t}, Ptr{Cvoid}), count, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Tuple(reverse!(count))
end

export steps_start
"""
    var_steps_start = steps_start(variable::Variable)
    var_steps_start::Union{Nothing,Int}

Read API, get available steps start from available steps count (e.g.
in a file for a variable).

This returns the absolute first available step, don't use with
`adios2_set_step_selection` as inputs are relative, use `0` instead.
"""
function steps_start(variable::Variable)
    steps_start = Ref{Csize_t}()
    err = ccall((:adios2_variable_steps_start, libadios2_c), Cint,
                (Ref{Csize_t}, Ptr{Cvoid}), steps_start, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Int(steps_start[])
end

export steps
"""
    var_steps = steps(variable::Variable)
    var_steps::Union{Nothing,Int}

Read API, get available steps count from available steps count (e.g.
in a file for a variable). Not necessarily contiguous.
"""
function steps(variable::Variable)
    steps = Ref{Csize_t}()
    err = ccall((:adios2_variable_steps, libadios2_c), Cint,
                (Ref{Csize_t}, Ptr{Cvoid}), steps, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Int(steps[])
end

export selection_size
"""
    var_selection_size = selection_size(variable::Variable)
    var_selection_size::Union{Nothing,Int}

Return the minimum required allocation (in number of elements of a
certain type, not bytes) for the current selection.
"""
function selection_size(variable::Variable)
    selection_size = Ref{Csize_t}()
    err = ccall((:adios2_variable_selection_size, libadios2_c), Cint,
                (Ref{Csize_t}, Ptr{Cvoid}), selection_size, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return Int(selection_size[])
end

"""
    var_min = minimum(variable::Variable)
    var_min::Union{Nothing,T}

Read mode only: return the absolute minimum for variable.
"""
function Base.minimum(variable::Variable)
    T = type(variable)
    T ≡ nothing && return nothing
    varmin = Ref{T}()
    err = ccall((:adios2_variable_min, libadios2_c), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), varmin, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return varmin[]
end

"""
    var_max = maximum(variable::Variable)
    var_max::Union{Nothing,T}

Read mode only: return the absolute maximum for variable.
"""
function Base.maximum(variable::Variable)
    T = type(variable)
    T ≡ nothing && return nothing
    varmax = Ref{T}()
    err = ccall((:adios2_variable_max, libadios2_c), Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}), varmax, variable.ptr)
    Error(err) ≠ error_none && return nothing
    return varmax[]
end


text:
const JOB_REGISTRY = (in_progress = DFC.config_path("jobs", "active.txt"), archived = DFC.config_path("jobs", "archived.txt"))

all_known_jobs() = vcat(readlines(JOB_REGISTRY[1]), readlines(JOB_REGISTRY[2]))

function cleanup_job_registry!()
    in_progress, archived = readlines(JOB_REGISTRY[1]), readlines(JOB_REGISTRY[2])
    for (ir, REG) in enumerate((in_progress, archived))
        stale_ids = findall(x -> !ispath(x) || !ispath(joinpath(x, "job.tt")), REG)
        if stale_ids !== nothing
            jobs_to_remove = REG[stale_ids]
            deleteat!(REG, stale_ids)
            unique!(REG)
        end
        write(JOB_REGISTRY[ir], join(REG, "\n"))
    end
    return (in_progress=in_progress, archived=archived)
end

function maybe_register_job(abspath::AbstractString)
    if abspath[end] == '/'
        abspath = abspath[1:end-1]
    end
    if ispath(abspath)
        n = occursin(".archived", abspath) ? JOB_REGISTRY.archived :
              JOB_REGISTRY.in_progress
        REG = readlines(n)
        jid = findfirst(isequal(abspath), REG)
        if jid === nothing
            push!(REG, abspath)
            write(n, join(REG, "\n"))
        end
    end
end
maybe_register_job(job::Job) = maybe_register_job(main_job_dir(abspath(job)))

"""
    registered_jobs(fuzzy::AbstractString = "")

Lists all the known [`Jobs`](@ref Job) directories that contain `fuzzy`.
Intended to be used as:
```julia
job_dirs = registered_jobs("NiO")
job = Job(job_dirs[1])
```
"""
function registered_jobs(fuzzy::AbstractString = "")
    choices = filter(x -> occursin(fuzzy, x), vcat(cleanup_job_registry!()...))
    timestamps = timestamp.(choices)
    sort_ids = sortperm(timestamps)
    return [(choices[i], timestamps[i]) for i in sort_ids]
end

function timestamp(jobdir::AbstractString)
    scriptpath = joinpath(jobdir, "job.tt")
    if ispath(scriptpath)
        return unix2datetime(mtime(scriptpath))
    else
        return DateTime(0)
    end
end

"""
    archived_jobs(fuzzy::AbstractString = "")

Returns a `Vector` of pairs with all archived [`Jobs`](@ref Job) whose directory contains the fuzzy as the first, and their description as the second item of the pair.
"""
function archived_jobs(fuzzy::AbstractString = "")
    jobs = filter(x -> occursin(fuzzy, x), cleanup_job_registry!().archived)
    return [j => ispath(joinpath(j, "description.txt")) ?
                 read(joinpath(j, "description.txt"), String) : "" for j in jobs]
end

"""
    load_jobs(fuzzy::AbstractString)

Loads all the known [`Jobs`](@ref Job) whose `dir` contains `fuzzy`.
"""
load_jobs(fuzzy::AbstractString) = map(x -> Job(x[1]), registered_jobs(fuzzy))

"""
    load_running_jobs(fuzzy::AbstractString = "")

Loads all [`Jobs`](@ref Job) that are currently running.
"""
function load_running_jobs(fuzzy::AbstractString = "")
    jobs = load_jobs(fuzzy)
    return filter(isrunning, jobs)
end

"""
    running_jobs(fuzzy::AbstractString = "")

Returns all the paths to running ['Jobs'](@ref Job).
"""
running_jobs(fuzzy::AbstractString = "") = getfield.(load_running_jobs(fuzzy), :dir)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import Base: bind

"""
    Executor

An executor is a realization of a symbolic architecture defined by a `SymbolicNode`.
The actual forward and backward computation specified by the network architecture can
be carried out with an executor.
"""
mutable struct Executor
  handle :: MX_ExecutorHandle
  symbol :: SymbolicNode
  arg_arrays  :: VecOfNDArray
  grad_arrays :: Vector{Union{Cvoid,<:NDArray}}
  aux_arrays  :: VecOfNDArray
  outputs     :: VecOfNDArray
  arg_dict    :: Dict{Symbol}
  aux_dict    :: Dict{Symbol}
end

function Executor(hdl::MX_ExecutorHandle, sym::SymbolicNode,
                  arg_arrays::VecOfNDArray, grad_arrays::AbstractVector,
                  aux_arrays::VecOfNDArray)
  # get output arrays
  ref_size = Ref{MX_uint}(0)
  ref_hdls = Ref{Ptr{MX_handle}}(C_NULL)
  @mxcall(:MXExecutorOutputs, (MX_handle, Ref{MX_uint}, Ref{Ptr{MX_handle}}),
          hdl, ref_size, ref_hdls)
  out_hdrs = unsafe_wrap(Array, ref_hdls[], ref_size[])
  out_arrays = [NDArray(MX_NDArrayHandle(x)) for x in out_hdrs]

  arg_names = list_arguments(sym)
  @assert(length(arg_names) == length(unique(arg_names)), "Duplicated names in arguments: $arg_names")
  arg_dict = Dict(zip(arg_names, arg_arrays))

  aux_names = list_auxiliary_states(sym)
  @assert(length(aux_names) == length(unique(aux_names)), "Duplicated names in auxiliary states: $aux_names")
  aux_dict = Dict(zip(aux_names, aux_arrays))

  Executor(hdl, sym, arg_arrays, grad_arrays, aux_arrays, out_arrays, arg_dict, aux_dict)
end

Base.unsafe_convert(::Type{MX_handle}, obj::Executor) =
  Base.unsafe_convert(MX_handle, obj.handle)
Base.convert(t::Type{MX_handle}, obj::Executor) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::Executor) = Base.unsafe_convert(t, obj)

function _get_ndarray_inputs(arg_key::AbstractString, args::VecOfNDArray,
                             arg_names::Vector{Symbol}, allow_missing::Bool)
  @assert(length(args) == length(arg_names), "Length of $arg_key does not match number of arguments")
  return (MX_handle[args...], args)
end

function _get_ndarray_inputs(arg_key::AbstractString, args::Dict{Symbol},
                             arg_names::Vector{Symbol}, allow_missing::Bool)
  args_vec = map(arg_names) do name
    arr = get(args, name, nothing)
    if !allow_missing
      @assert(!isa(arr, Cvoid), "Must specify all arguments in $arg_key ($name is missing)")
    end
    arr
  end
  # help the type inference
  if allow_missing
    args_vec = Union{NDArray,Cvoid}[args_vec...]
  else
    args_vec = NDArray[args_vec...]
  end
  args_hdr = MX_handle[(isa(x,Cvoid) ? MX_handle(0) : x) for x in args_vec]
  return (args_hdr, args_vec)
end

"""
    bind(sym, ctx, args; args_grad=Dict(), aux_states=Dict(), grad_req=GRAD_WRITE)

Create an `Executor` by binding a `SymbolicNode` to concrete `NDArray`.

# Arguments
* `sym::SymbolicNode`: the network architecture describing the computation graph.
* `ctx::Context`: the context on which the computation should run.
* `args`: either a list of `NDArray` or a dictionary of name-array pairs. Concrete
          arrays for all the inputs in the network architecture. The inputs typically include
          network parameters (weights, bias, filters, etc.), data and labels.
          See [`list_arguments`](@ref) and [`infer_shape`](@ref).
* `args_grad`: a `Vector` of `NDArray` or a `Dict` contains `NDArray`
* `aux_states`: a `Vector` of `NDArray` or a `Dict` contains `NDArray`
* `grad_req`: single value, a `Vector` of `GRAD_REQ` or a `Dict{Symbol,GRAD_REQ}`
"""
function bind(self::SymbolicNode, ctx::Context, args;
              args_grad = Dict{Symbol,NDArray}(),
              aux_states = Dict{Symbol,NDArray}(),
              grad_req = GRAD_WRITE)

  arg_names = list_arguments(self)

  args_hdr, args           = _get_ndarray_inputs("args", args, arg_names, false)
  args_grad_hdr, args_grad = _get_ndarray_inputs("args_grad", args_grad, arg_names, true)
  aux_args_hdr, aux_states = _get_ndarray_inputs("aux_states", aux_states, list_auxiliary_states(self), false)

  if isa(grad_req, GRAD_REQ)
    reqs = MX_uint[MX_uint(grad_req) for i=1:length(args)]
  elseif isa(grad_req, Vector{GRAD_REQ})
    @assert(length(grad_req) == length(args))
    reqs = MX_uint[MX_uint.(grad_req)...]
  elseif isa(grad_req, Dict{Symbol, GRAD_REQ})
    reqs = MX_uint[MX_uint(get(grad_req, name, GRAD_NOP)) for name in arg_names]
  end

  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXExecutorBind,
          (MX_handle, Cint, Cint, MX_uint, Ptr{MX_handle}, Ptr{MX_handle}, Ptr{MX_uint},
           MX_uint, Ptr{MX_handle}, Ref{MX_handle}),
          self, ctx.device_type, ctx.device_id, length(args), args_hdr,
          args_grad_hdr, reqs, length(aux_states), aux_args_hdr, ref_hdr)
  args_grad = convert(Vector{Union{Cvoid,NDArray}}, args_grad)
  executor = Executor(MX_ExecutorHandle(ref_hdr[]), self,
                      args, args_grad, aux_states)
end

function bind(x::SymbolicNode; context::Context = cpu(), kwargs...)
  kwargs = Dict(kwargs)
  @assert(haskey(kwargs, :args), "Must specify args")
  args = pop!(kwargs, :args)
  bind(x, context, args; kwargs...)
end

function simple_bind(self::SymbolicNode, ctx::Context;
                     grad_req::Union{GRAD_REQ,Dict{Symbol,GRAD_REQ}} = GRAD_WRITE,
                     kwargs...)
  arg_shapes, out_shapes, aux_shapes = infer_shape(self; kwargs...)
  @assert(!isa(arg_shapes, Cvoid), "Information not enough to perform complete shape inference")

  arg_arrays = NDArray[zeros(shape, ctx) for shape in arg_shapes]
  arg_names  = list_arguments(self)

  grad_arrays = Dict{Symbol,NDArray}()

  if grad_req != GRAD_NOP
    shapes = zip(arg_names, arg_shapes)

    # if not in provided data, should be parameters
    provided_data_names = [x[1] for x in kwargs]
    shapes = filter(x -> !in(x[1], provided_data_names), shapes)

    # Remove all gradients for nop params
    # if isa(grad_req, Dict{Symbol, GRAD_REQ})
    #  shapes = filter(x -> grad_req[x[1]] != GRAD_NOP,shapes)
    # end

    for (name, shape) in shapes
      grad_arrays[name] = zeros(shape, ctx)
    end
  end

  aux_arrays = [zeros(shape, ctx) for shape in aux_shapes]
  return bind(self, ctx, arg_arrays, args_grad=grad_arrays, grad_req=grad_req, aux_states=aux_arrays)
end

function forward(self::Executor; is_train::Bool = false, kwargs...)
  for (k,v) in kwargs
    @assert(k ∈ self.arg_dict, "Unknown argument $k")
    @assert(isa(v, NDArray), "Keyword argument $k must be an NDArray")
    copy!(self.arg_dict[k], v)
  end

  @mxcall(:MXExecutorForward, (MX_handle, Cint), self, is_train)

  self.outputs
end

backward(x::Executor) = backward(x, NDArray[])
backward(x::Executor, out_grad::NDArray) = backward(x, [out_grad])
backward(x::Executor, out_grads::VecOfNDArray) =
  @mxcall(:MXExecutorBackward, (MX_handle, MX_uint, Ptr{MX_handle}),
          x, length(out_grads), MX_handle[out_grads...])

function copy_params_from(self::Executor, arg_params::Dict{Symbol},
                          aux_params::Dict{Symbol} = Dict{Symbol,Any}();
                          allow_extra_params::Bool = false)
  for (name, array) in arg_params
    if haskey(self.arg_dict, name)
      copy!(self.arg_dict[name], array)
    else
      @assert(allow_extra_params, "Extra params $name not in the arguments")
    end
  end

  for (name, array) in aux_params
    if haskey(self.aux_dict, name)
      copy!(self.aux_dict[name], array)
    else
      @assert(allow_extra_params, "Extra auxiliary state $name not recognized")
    end
  end
end

Base.show(io::IO, x::Executor) =
  print(io, "mx.", split(string(typeof(x)), '.')[end], " ", x.handle.value)

"""
    print([io::IO], x::Executor)

Get a debug string about internal execution plan.

Can be used to get an estimated about the memory cost.

```julia
julia> x = mx.Variable(:x)
MXNet.mx.SymbolicNode x

julia> exec = mx.bind(x + 1, mx.cpu(), Dict(:x => mx.ones(2,3)))
mx.Executor Ptr{Nothing} @0x000055c3dee9eb30

julia> print(exec)
Symbol Outputs:
        output[0]=_plus_scalar0(0)
Variable:x
--------------------
Op:_plus_scalar, Name=_plus_scalar0
Inputs:
        arg[0]=x(0) version=0
Attrs:
        scalar=1.00000000e+00
Total 0 MB allocated
Total 11 TempSpace resource requested
```
"""
Base.print(io::IO, x::Executor) = print(io, debug_str(x))
Base.print(x::Executor)         = print(STDOUT, x)

function debug_str(x::Executor)
  s_ref = Ref{Cstring}(C_NULL)
  @mxcall(:MXExecutorPrint, (MX_handle, Ptr{Cstring}), x.handle, s_ref)
  unsafe_string(s_ref[])
end

max_stars_repo_name: allisonpaily/incubator-mxnet

max_issues_repo_name: ptrendx/mxnet

max_forks_repo_name: ptrendx/mxnet


text:
import Base: show

function show(io::IO, layer::T) where T <: SlothLayer
    head = get_layer_type(layer)
    args = get_layer_args(layer)
    body = join(["$k=$v" for (k,v) in args], ", ")
    str = "$head($body)"
    print(str)
end

function show(io::IO, sequential::Sequential; level=1, num_spaces=2)
    println("")
    println(repeat(" ", level*num_spaces) * "Sequential(")
    num_layers = length(sequential.layers)
    for (i,l) in enumerate(sequential.layers)
        print(repeat(" ", (1+level)*num_spaces))
        if isa(l, Sequential)
            show(io, l; level=level+1, num_spaces=num_spaces)
        else
            show(io, l)
        end
        i == num_layers && continue
        println(",")
    end
    println("")
    print(repeat(" ", level*num_spaces) * ")")
end

get_layer_type(l::SlothLayer) = titlecase(split(string(typeof(l)), ".")[end])
get_layer_type(l::Pool) = ("MaxPool", "PaddedAvgPool", "AvgPool")[Int(l.mode)+1]
get_layer_type(l::Activation) = titlecase(split(string(l.f), ".")[end])

function get_layer_args(layer::Linear)
    args = Any[
        ("in_features", layer.in_features),
        ("out_features", layer.out_features)]
    isa(layer.bias, AbstractFloat) && push!(args, ("bias", false))
    return args
end

function get_layer_args(layer::Union{Conv, ConvTranspose})
    args = Any[
        ("in_channels", layer.in_channels),
        ("out_channels", layer.out_channels),
        ("kernel_size",layer.kernel_size)]
    layer.mode != 0 && push!(args, ("mode", layer.mode))
    layer.padding != 0 && push!(args, ("padding", layer.padding))
    layer.stride != 1 && push!(args, ("stride", layer.stride))
    layer.bias == F(0.0) && push!(args, ("bias", false))
    return args
end

function get_layer_args(layer::Pool)
    args = Any[]
    Int(layer.mode) != 0 && push!("mode", layer.mode)
    push!(args, ("window", layer.window))
    layer.padding != 0 && push!(args, ("padding", layer.padding))
    layer.stride != layer.window && push!(args, ("stride", layer.stride))
    return args
end

get_layer_args(layer::BatchNorm) = Any[("input", layer.num_features)]

get_layer_args(layer::Embedding) = Any[
    ("num_embeddings", layer.num_embeddings),
    ("embedding_dim", layer.embedding_dim)]

get_layer_args(layer::Dropout) = Any[("p", layer.p)]

get_layer_args(l::Activation) = Any[]


text:
# TODO:
# - reorganize this
# - test for `include_weave`
# - fire horrible tests
# - test for ipynb integration
# - test for integrations with other libraries, especially for Plots.jl and Gadfly.jl

# %%
using Weave, Test
using Weave: WeaveDoc, run_doc, get_format

function mock_doc(str, informat = "markdown")
    f = tempname()
    write(f, str)
    return WeaveDoc(f, informat)
end
mock_run(str, informat = "markdown"; kwargs...) = run_doc(mock_doc(str, informat); kwargs...)

function test_mock_weave(test_function, str; kwargs...)
    f = tempname()
    write(f, str)
    f = weave(f; kwargs...)
    try
        weave_body = read(f, String)
        test_function(weave_body)
    catch
        rethrow()
    finally
        rm(f)
    end
end

# %%
@testset "Weave" begin
    @testset "reader" begin
        include("reader/test_chunk_options.jl")
        include("reader/test_inline.jl")
    end

    @testset "header processing" begin
        include("test_header.jl")
    end

    @testset "run" begin
        include("run/test_module.jl")
        include("run/test_meta.jl")
        include("run/test_error.jl")
    end

    @testset "render" begin
        include("render/texformats.jl")
    end

    @testset "conversions" begin
        include("test_converter.jl")
    end

    @testset "display" begin
        include("test_display.jl")
    end

    @testset "end2end" begin
        include("end2end/test_end2end.jl")
    end

    @testset "legacy" begin
        include("markdown_test.jl")
        include("render_figures_test.jl")
        include("cache_test.jl")
    end
end