text:
module TestFMDatasetML100K

using FactorizationMachines
using Base.Test

info("Running tests on ml-100k dataset...")
(ml100kTrain, ml100kTrainRatings) = @time read_libsvm("data/ml-100k.train.libfm")
info("Train dim: $(size(ml100kTrain))")
(ml100kTest, ml100kTestRatings) = @time read_libsvm("data/ml-100k.test.libfm")
info("Test dim: $(size(ml100kTest))")

info("Test ML-100K - Training model...")
fmMl100k = @time train(ml100kTrain, ml100kTrainRatings,
        method = Methods.sgd(num_epochs = 10, alpha = 0.1),
        model_params = Models.gauss(num_factors = 4))
#info("FM Model: $fmMl100k")
info("Test ML-100K - Prediction over test dataset...")
p = @time predict(fmMl100k, ml100kTest)

@test length(p) == length(ml100kTestRatings)

rmse = sqrt(sum((p-ml100kTestRatings).^2)/length(p))
info("RMSE on test data: $rmse")

if isfile("data/ml-100k.output.libfm")
    info("File with output from libfm exists. Bencharking to its result...")
    libfmOutput = readcsv("data/ml-100k.output.libfm")
    rmseWithLibfm = sqrt(sum((p-libfmOutput).^2)/length(p))
    info("RMSE with Libfm output: $rmseWithLibfm")
end

end


text:
"""
plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1,j=length(x))

plot the result from the [`check_gradient`](@ref)
on data `x,y` with two comparison lines

1) `line_base`+ t`slope`  as the goal the plot should have
2) `a`+ `b`t` on the interval [`x[i]`, `x[j]`] for some (best fitting) comparison slope
"""
function plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1, j=length(x))
    fig = plot(
        x,
        y;
        xaxis=:log,
        yaxis=:log,
        label="\$E(t)\$",
        linewidth=3,
        legend=:topleft,
        color=:lightblue,
    )
    s_line = [exp10(line_base + t * slope) for t in log10.(x)]
    plot!(fig, x, s_line; label="slope s=2", linestyle=:dash, color=:black, linewidth=2)
    if (i != 0) && (j != 0)
        best_line = [exp10(a + t * b) for t in log10.(x[i:j])]
        plot!(
            fig,
            x[i:j],
            best_line;
            label="best slope $(@sprintf("%.4f", b))",
            color=:blue,
            linestyle=:dot,
            linewidth=2,
        )
    end
    gui(fig)
    return nothing
end

ext: jl


text:
using Decimals
using Compat.Test

@testset "Rounding" begin

@test round(Decimal(7.123456), digits=0) == Decimal(7)
@test round(Decimal(7.123456), digits=2) == Decimal(7.12)
@test round(Decimal(7.123456), digits=3) == Decimal(7.123)
@test round(Decimal(7.123456), digits=5) == Decimal(7.12346)
@test round(Decimal(7.123456), digits=6) == Decimal(7.123456)

@test round.(Decimal.([0.1111, 0.2222, 0.8888]), digits=2) == Decimal.([0.11, 0.22, 0.89])

function tet()
    a = parse(Decimal, "1.0000001")
    for i = 1:27
        a *= a
    end
    return a
end

# set DIGITS = 20 (aaljuffali's example)
@test tet() == Decimal(0, 67453047074102193157641340, -20)

end


text:
using   TexTables
import  TexTables:  TableIndex, TableDict
using   DataStructures
using   DataFrames
using   Test
using   Random
using   LinearAlgebra
using   Statistics
using   RDatasets
using   CSV

include("helper.jl")

tests = [   "tablecol", "composite_tables", "indexing", "printing",
            "examples", "quicktools", "formatted_numbers"]

@testset "TexTables" begin
    for testsuite in tests
        @testset "$testsuite" begin
            include("$testsuite.jl")
        end
    end
end


text:
using CSV
using DataStructures: SortedDict, SortedSet
using Random
using DataAPI
using StringEncodings
using StaticArrays: SVector
using CoordinateTransformations: AbstractAffineMap

Base.convert(::Type{Symbol}, elm::Element) = Symbol(uppercase(elm.symbol))

"""
Handles RJ Lee-style Zeppelin particle data sets.
"""
struct Zeppelin
    headerfile::String
    header::SortedDict{String,String}
    elms::SortedSet{Element}
    classes::SortedDict{Int, String}
    data::DataFrame

    function Zeppelin( #
        headerfile::String,
        header::AbstractDict{String,String},
        elms::Union{AbstractSet{Element}, AbstractVector{Element}},
        classes::AbstractDict{String,String},
        data::DataFrame
    )
        function _massagehdz(header)
            res = SortedDict(header)
            foreach(hi -> if hi in keys(header) delete!(res, hi) end, ( "PARTICLE_PARAMETERS", "PARAMETERS", "HEADER_FMT" ))
            foreach(f->if startswith(f,"ELEM") delete!(res,f) end, keys(header))
            foreach(f->if startswith(f,"CLASS") delete!(res,f) end, keys(header))
            return res
        end
        # Mostly converts columns to categorical as desired...
        function _massagepxz(header, pxz)::DataFrame
            res = copy(pxz)
            sortclasses(c1, c2) = isless(parse(Int, c1[6:end]), parse(Int, c2[6:end]))
            sortedkeys = sort(collect(filter(c -> !isnothing(match(r"^CLASS\d+", c)), collect(keys(header)))), lt = sortclasses)
            clsnames = append!(["--None--"], map(c -> header[c], sortedkeys))
            for col in ("CLASS", "VERIFIEDCLASS")
                ic = findfirst(nm->nm==col, names(res))
                if !isnothing(ic)
                    cls = categorical(Union{String,Missing}[ clsnames... ], ordered=true)[1:0] # empty but with correct levels
                    foreach(name->push!(cls, name), map(cl -> get(clsnames, convert(Int, cl) + 2, "####"), pxz[:, col]))
                    select!(res, Not(ic))
                    insertcols!(res, ic, col=>cls)
                end
            end
            for col in ( "FIRSTELM", "SECONDELM", "THIRDELM", "FOURTHELM" )
                ic = findfirst(nm->nm==col, names(res))
                if !isnothing(ic)
                    select!(res, Not(ic))
                    insertcols!(res, ic, col=>map(z-> z in 1:length(elements) ? elements[z] : missing, pxz[:,col]))
                end
            end
            for col in ( "XDAC", "YDAC", "TYPE_4ET_", "TYPE4ET" )
                ic = findfirst(nm->nm==col, names(res))
                if !isnothing(ic)
                    select!(res, Not(ic))
                end
            end
            return res
        end
        function _massageclasses(classes::AbstractDict{String, String})
            return SortedDict{Int, String}( parse(Int, key[6:end]) => val for (key, val) in classes )
        end
        new(
            headerfile,
            _massagehdz(header),
            SortedSet(elms),
            _massageclasses(classes),
            _massagepxz(header, data)
        )
    end

    Zeppelin( hdzfilename::String) = loadZep( hdzfilename)

    function Zeppelin( #
        headerfile::String,
        header::AbstractDict{String,String},
        data::DataFrame,
    )
        Zeppelin(
            headerfile,
            header,
            filter(elm -> uppercase(elm.symbol) in names(data), PeriodicTable.elements[1:94]),
            Dict{String,String}(),
            data,
        )
    end
end

Base.copy(z::Zeppelin) = Zeppelin(z.headerfile, copy(z.header), copy(z.elms), copy(z.classes), copy(z.data))

function loadZep( hdzfilename::String)
    remapcolumnnames = Dict(
        "PART#" => "NUMBER",
        "PARTNUM" => "NUMBER",
        "FIELD#" => "FIELD",
        "FIELDNUM" => "FIELD",
        "MAGFIELD#" => "MAGFIELD",
        "MAGFIELDNUM" => "MAGFIELD",
        "X_ABS" => "XABS",
        "Y_ABS" => "YABS",
        "X_DAC" => "XDAC",
        "Y_DAC" => "YDAC",
        "XCENT" => "XDAC",
        "YCENT" => "YDAC",
        "X_FERET" => "XFERET",
        "Y_FERET" => "YFERET",
        "DVG" => "DAVG",
        "DAVE" => "DAVG",
        "PERIM" => "PERIMETER",
        "ORIENT" => "ORIENTATION",
        "LIVE_TIME" => "LIVETIME",
        "FIT_QUAL" => "FITQUAL",
        "MAG_INDEX" => "MAGINDEX",
        "FIRST_ELEM" => "FIRSTELM",
        "SECOND_ELEM" => "SECONDELM",
        "THIRD_ELEM" => "THIRDELM",
        "FOURTH_ELEM" => "FOURTHELM",
        "ATOMICNUMBER1" => "FIRSTELM",
        "ATOMICNUMBER2" => "SECONDELM",
        "ATOMICNUMBER3" => "THIRDELM",
        "ATOMICNUMBER4" => "FOURTHELM",
        "FIRST_CONC" => "COUNTS1",
        "SECOND_CONC" => "COUNTS2",
        "THIRD_CONC" => "COUNTS3",
        "FOURTH_CONC" => "COUNTS4",
        "FIRST_PCT" => "FIRSTPCT",
        "SECOND_PCT" => "SECONDPCT",
        "THIRD_PCT" => "THIRDPCT",
        "FOURTH_PCT" => "FOURTHPCT",
        "PCT1" => "FIRSTPCT",
        "PCT2" => "SECONDPCT",
        "PCT3" => "THIRDPCT",
        "PCT4" => "FOURTHPCT",
        "TYPE(4ET)#" => "TYPE4ET",
        "TYPE(4ET)" => "TYPE4ET",
        "VOID_AREA" => "VOIDAREA",
        "RMS_VIDEO" => "RMSVIDEO",
        "FIT_QUAL" => "FITQUAL",
        "VERIFIED_CLASS" => "VERIFIEDCLASS",
        "EDGE_ROUGHNESS" => "EDGEROUGHNESS",
        "COMP_HASH" => "COMPHASH",
        "PSEM_CLASS" => "CLASS",
    )
    columnnames(cols) = uppercase.(map(cn -> get(remapcolumnnames, cn, cn), map(c -> c[1], cols)))
    header, columns, hdr = Dict{String,String}(), [], true
    open(hdzfilename, enc"WINDOWS-1252","r") do f
        for line in readlines(f)
            if hdr
                p = findfirst(c -> c == '=', line)
                if !isnothing(p)
                    (k, v) = line[1:p-1], line[p+1:end]
                    header[k] = v
                    hdr = !isequal(uppercase(k), "PARTICLE_PARAMETERS")
                end
            else
                push!(columns, string.(strip.(split(line, "\t"))))
            end
        end
    end
    pxz = CSV.File(
        replace(hdzfilename, r".[h|H][d|D][z|Z]$" => ".pxz"),
        header = columnnames(columns),
        delim='\t',
        normalizenames = true,
        missingstring="-",
        decimal='.'
    ) |> DataFrame
    #elms = filter(z -> uppercase(z.symbol) in names(pxz), PeriodicTable.elements[1:94])
    #classes = Dict{String, String}(key=>header[key] for key in filter(f->startswith(f, "CLASS") && (f ≠ "CLASSES"),keys(header)))
    return Zeppelin( hdzfilename, header, pxz)
end

function Base.show(io::IO, zep::Zeppelin)
    print(io, "Zeppelin[$(zep.headerfile),$(size(zep.data))]")
end

Base.size(z::Zeppelin) = size(z.data)

classes(zep::Zeppelin) = zep.classes

NeXLCore.elms(zep::Zeppelin) = zep.elms

header(zep::Zeppelin) = SortedDict(zep.header)

"""
    asa(::Type{DataFrame}, zep::Zeppelin; rows=missing, columns=missing, sortcol=:None, rev=false)

Create a DataFrame containing the particle data.  Optionally specify which rows and columns to include and
whether to sort the resulting table.
"""
function NeXLUncertainties.asa(::Type{DataFrame}, zep::Zeppelin; rows=missing, sortcol=:None, columns=missing, rev=false)
    rs = ismissing(rows) ? eachparticle(zep) : intersect(eachparticle(zep), rows)
    cols = ismissing(columns) ? names(zep.data) : intersect(columns,names(zep.data))
    res = zep.data[rs, cols]
    return sortcol≠:None ? sort(res, sortcol, rev=rev) : res
end

eachparticle(zep::Zeppelin) = 1:size(zep.data, 1)

function DataAPI.describe(zeps::AbstractVector{Zeppelin}; dcol=:DAVG, nelms=2)
    function build(zep)
        df=describe(zep, dcol=dcol, nelms=nelms)
        insertcols!(df, 1, :Dataset=>[zep.header["DESCRIPTION"] for _ in 1:size(df,1)])
        return df
    end
    return vcat(map(z->build(z),zeps))
end

function DataAPI.describe(zep::Zeppelin; dcol=:DAVG, nelms=3)
    function sortedstats(rows)
        tmp = [ ( elm, summarystats(zep[rows,convert(Symbol, elm)])) for elm in zep.elms ]
        return sort!(tmp, lt=(v1,v2)->isless(v1[2].mean,v2[2].mean),rev=true)
    end
    cnelms = min(length(zep.elms), nelms)
    ds = summarystats(zep[:,dcol])
    clss = String["All"]
    szs = Int[size(zep.data,1)]
    minds = Float64[ ds.min ]
    medds = Float64[ ds.median ]
    maxds = Float64[ ds.max ]
    sbcm = StatsBase.countmap(zep[:,:CLASS])
    for (cn, _) in sbcm
        rows = rowsclass(zep,String(cn))
        #if length(rows)>0
        push!(clss, String(cn))
        push!(szs, length(rows))
        ds = summarystats(zep[rows,dcol])
        push!(minds, ds.min)
        push!(medds, ds.median)
        push!(maxds, ds.max)
        #end
    end
    dfres = DataFrame(Symbol("Class")=>clss, Symbol("Count")=>szs,Symbol("Min[$dcol]")=>minds, Symbol("Median[$dcol]")=>medds, Symbol("Max[$dcol]")=>maxds)
    # Add All elemental row
    ss = sortedstats(eachparticle(zep))
    elmdfs=DataFrame[ ]
    uem, uef = Union{String,Missing}, Union{Float64,Missing}
    for ne in 1:cnelms
        cols = Symbol.( [ "Elm[$ne]", "Min[$ne]", "Median[$ne]", "Max[$ne]" ])
        sst = ss[ne][2]
        push!(elmdfs, DataFrame(cols[1]=>uem[ symbol(ss[ne][1]) ], cols[2]=>uef[sst.min], cols[3]=>uef[sst.median], cols[4]=>uef[sst.max]))
    end
    for (cn, _) in sbcm
        try
            ss = sortedstats(rowsclass(zep, String(cn)))
            for ne in 1:cnelms
                sst = ss[ne][2]
                push!(elmdfs[ne], [ symbol(ss[ne][1]), sst.min, sst.median, sst.max ])
            end
        catch
            for ne in 1:cnelms
                push!(elmdfs[ne], [ missing, missing, missing, missing ])
            end
        end
    end
    for elmdf in elmdfs
        for col in names(elmdf)
            insertcols!(dfres, size(dfres,2)+1, col=>elmdf[:,col])
        end
    end
    return sort!(dfres, [:Count, :Class], rev=true)
end

"""
    zep[123] # where zep is a Zeppelin

Returns the Spectrum (with images) associated with the particle at row
"""
Base.getindex(zep::Zeppelin, row::Int) = spectrum(zep, row, true)

Base.getindex(zep::Zeppelin, rows, cols) = Base.getindex(zep.data, rows, cols)

Base.lastindex(zep::Zeppelin, axis::Integer) = Base.lastindex(zep.data, axis)

# Replace the default in PeriodicTable because it is too verbose...
# Base.show(io::IO, elm::Element) = print(io, elm.symbol)

function filenumber(zep::Zeppelin, row::Int)::String
    tmp="$(zep.data[row,:NUMBER])"
    return repeat('0',max(0,5-length(tmp)))*tmp
end

"""
    spectrumfilename(zep::Zeppelin, row::Int, dir::AbstractString="MAG", ext::AbstractString=".tif")

Returns the name of the spectrum/image file for the particle in the specified row.
"""
function spectrumfilename(zep::Zeppelin, row::Int, dir::AbstractString = "MAG", ext::AbstractString = ".tif", relocated=false)
    mag = hasproperty(zep.data, :MAG) ? convert(Int,trunc(zep.data[row, :MAG])) : 0
    # First check if a spectrum file exists
    tmp=repr(zep.data[row,:NUMBER])
    if relocated && isdir(joinpath(dirname(zep.headerfile),"RELOCATED"))
        for n in 6:-1:4
            fn = joinpath(dirname(zep.headerfile), "RELOCATED", repeat('0',max(0,n-length(tmp)))*tmp*ext)
            if isfile(fn)
                return fn
            end
        end
    end
    for n in 6:-1:4
        fn = joinpath(dirname(zep.headerfile), "$(dir)$(mag)", repeat('0',max(0,n-length(tmp)))*tmp*ext)
        if isfile(fn)
            return fn
        end
    end
    # Default case
    return joinpath(dirname(zep.headerfile), "$(dir)$(mag)", repeat('0',max(0,5-length(tmp)))*tmp*ext)
end

"""
    ParticleClassifier

A type that implements:

    classify(zep::Zeppelin, sr::ParticleClassifier)::CategoricalArray{String}
"""
abstract type ParticleClassifier end

"""
    spectrum(zep::Zeppelin, row::Int, withImgs = true)::Union{Spectrum, missing}

Returns the Spectrum (with images) associated with the particle at row.  If withImgs
is true, the associated image or images are read.
"""
function spectrum(zep::Zeppelin, row::Int, withImgs = true, relocated=true)::Union{Spectrum,Missing}
    file, at = spectrumfilename(zep, row, "MAG", ".tif", relocated), missing
    if isfile(file)
        try
            at = loadspectrum(ASPEXTIFF, file; withImgs = withImgs)
        catch err
            showerror(stderr, err)
            @info "$(file) does not appear to be a valid ASPEX spectrum TIFF."
        end
        try
            at[:BeamEnergy] = beamenergy(zep, get(at, :BeamEnergy, 20.0e3))
            at[:ProbeCurrent] = get(at, :ProbeCurrent, probecurrent(zep, 1.0))
            at[:Signature] = filter(kv->kv[2]>0.0, Dict(elm => zep.data[row, convert(Symbol,elm)] for elm in zep.elms))
            at[:Name] = "P[$(zep.data[row, :NUMBER]), $(zep.data[row, :CLASS])]"
        catch
            @info "Error adding properties to ASPEX TIFF file."
        end
    end
    return at
end

function iszeppelin(filename::String)
    res=false
    open(filename) do ios
        seekstart(ios)
        if isequal(uppercase(String(read(ios,11))),"PARAMETERS=")
            readline(ios) # read the rest of the line
            res = isequal(uppercase(readline(ios)),"HEADER_FMT=ZEPP_1")
        end
    end
    return res
end

function writeZep(zep::Zeppelin,  hdzfilename::String)
    remapcolumnnames = Dict(
        "NUMBER"=> "PART#\t1\tINT16",
        "FIELD" => "FIELD#\t1\tINT16",
        "MAGFIELD" => "MAGFIELD#\t1\tINT16",
        "XABS" => "X_ABS\tmm\tFLOAT",
        "YABS" => "Y_ABS\tmm\tFLOAT",
        "XDAC" => "X_DAC\t1\tINT16",
        "YDAC" => "Y_DAC\t1\tINT16",
        "XFERET" => "X_FERET\tµm\tFLOAT",
        "YFERET" => "Y_FERET\tµm\tFLOAT",
        "DAVG" => "DAVE\tµm\tFLOAT",
        "DMAX" => "DMAX\tµm\tFLOAT",
        "DMIN" => "DMIN\tµm\tFLOAT",
        "DPERP" => "DPERP\tµm\tFLOAT",
        "ASPECT" =>"ASPECT\t1\tFLOAT",
        "AREA" => "AREA\tµm²\tFLOAT",
        "PERIMETER" => "PERIMETER\tµm\tFLOAT",
        "ORIENTATION" => "ORIENTATION\tdeg\tFLOAT",
        "LIVETIME" => "LIVE_TIME\ts\tFLOAT",
        "FITQUAL" => "FIT_QUAL\t1\tFLOAT",
        "MAG" => "MAG\t1\tINT16",
        "VIDEO" => "VIDEO\t1\tINT16",
        "IMPORTANCE" => "IMPORTANCE\t1\tINT16",
        "COUNTS" => "COUNTS\t1\tFLOAT",
        "MAGINDEX" => "MAG_INDEX\t1\tINT16",
        "FIRSTELM" => "FIRST_ELEM\t1\tINT16",
        "SECONDELM" => "SECOND_ELEM\t1\tINT16",
        "THIRDELM" => "THIRD_ELEM\t1\tINT16",
        "FOURTHELM" => "FOURTH_ELEM\t1\tINT16",
        "COUNTS1" => "FIRST_CONC\tcounts\tFLOAT",
        "COUNTS2" => "SECOND_CONC\tcounts\tFLOAT",
        "COUNTS3" => "THIRD_CONC\tcounts\tFLOAT",
        "COUNTS4" => "FOURTH_CONC\tcounts\tFLOAT",
        "FIRSTPCT" => "FIRST_PCT\t%\tFLOAT",
        "SECONDPCT" => "SECOND_PCT\t%\tFLOAT",
        "THIRDPCT" => "THIRD_PCT\t%\tFLOAT",
        "FOURTHPCT" => "FOURTH_PCT\t%\tFLOAT",
        "TYPE4ET" => "TYPE(4ET#)\t1\tLONG",
        "VOIDAREA" => "VOID_AREA\tµm²\tFLOAT",
        "RMSVIDEO" => "RMS_VIDEO\t1\tINT16",
        "FITQUAL" => "FIT_QUAL\t1\tFLOAT",
        "VERIFIEDCLASS" => "VERIFIED_CLASS\t1\tINT16",
        "EDGEROUGHNESS" => "EDGE_ROUGHNESS\t1\tFLOAT",
        "COMPHASH" => "COMP_HASH\t1\tLONG",
        "CLASS" => "PSEM_CLASS\t1\tINT16",
        "TYPE_4ET_" => "Type[4ET]\t1\tLONG",
    )
    merge!(remapcolumnnames, Dict( uppercase(elm.symbol) => "$(uppercase(elm.symbol))\t%(k)\tFLOAT" for elm in zep.elms))
    merge!(remapcolumnnames, Dict( "U_$(uppercase(elm.symbol))_"  => "U[$(uppercase(elm.symbol))]\t%(k)\tFLOAT" for elm in zep.elms))
    headeritems = copy(zep.header)
    # add back the element tags
    for (i,elm) in enumerate(zep.elms)
        headeritems["ELEM$(i-1)"] = "$(elm.symbol) $(z(elm)) 1"
    end
    headeritems["ELEMENTS"] = "$(length(zep.elms))"
    if "CLASS" in names(zep.data)
        clsdata=zep.data[:,:CLASS]
        for (i,cls) in enumerate(levels(clsdata))
            # println("CLASS$(i-1) => $cls")
            headeritems["CLASS$(i-1)"]=cls
            headeritems["CLASSES"]="$(i+1)"
        end
    end
    headeritems["TOTAL_PARTICLES"] = "$(size(zep.data,1))"
    # write out the header
    colnames = names(zep.data)
    open(hdzfilename,"w") do ios
        println(ios, "PARAMETERS=$(length(headeritems)+size(zep.data,2)+2)")
        println(ios, "HEADER_FMT=ZEPP_1")
        foreach(hk->println(ios,"$(hk)=$(headeritems[hk])"), sort(collect(keys(headeritems))))
        println(ios, "PARTICLE_PARAMETERS=$(size(zep.data,2))")
        foreach(colnames) do col
            println(ios, get(remapcolumnnames, col, "$col\t1\t"*uppercase(repr(eltype(zep.data[:,col])))))
        end
    end
    pxzfilename = replace( hdzfilename, r".[h|H][d|D][z|Z]"=>".pxz")
    zd = DataFrame(zep.data)
    # Replace categorical data with the index (either into CLASS# in header or z(elm))
    for (ic, col) in enumerate(names(zd))
        if col in ("CLASS", "VERIFIEDCLASS", "FIRSTELM", "SECONDELM", "THIRDELM", "FOURTHELM" )
            coldata=zd[:,col]
            select!(zd,Not(ic))
            if col in ( "CLASS", "VERIFIEDCLASS" )
                @assert coldata isa CategoricalArray "Whoops! $col is a $(typeof(zd[:,col]))"
                insertcols!(zd, ic, col=>[ coldata.refs[i]-1  for i in eachindex(coldata) ]) # CLASS
            elseif col in ("FIRSTELM", "SECONDELM", "THIRDELM", "FOURTHELM" )
                insertcols!(zd, ic, col=>[ ismissing(elm) ? 0 : z(elm)  for elm in coldata ]) # Element
            else
                @info "Unanticipated column $col"
            end
        end
    end
    CSV.write(pxzfilename, zd, delim="\t", missingstring="-", header=false)
end

function beamenergy(zep::Zeppelin, def=missing)
    number(v) = parse(Float64, match(r"([+-]?[0-9]+[.]?[0-9]*)",v)[1])
    val = def
    try
        v = get(zep.header, "ACCELERATING_VOLTAGE", missing)
        val = !ismissing(v) ? number(v)*1.0e3 : def
    catch
        # Ignore it...
    end
    return val
end

function probecurrent(zep::Zeppelin, def=missing)
    number(v) = parse(Float64, match(r"([+-]?[0-9]+[.]?[0-9]*)",v)[1])
    val = def
    try
        v = get(zep.header, "PROBE_CURRENT", missing)
        val = !ismissing(v) ? number(v) : def
    catch
        # Ignore it...
    end
    return val
end

function magdata(zep::Zeppelin, index::Int)
    h = split(zep.header["MAG_FMT"],isspace)
    v = split(zep.header["MAG$index"],isspace)
    return Dict( string(h[i])=>parse(Float64, v[i]) for i in 1:min(length(h),length(v)) )
end

"""
    randomsubset(zep::Zeppelin, maxRows::Int)

Creates an ordered random subselection of `maxRows` without replacement of the rows in the `zep` dataset.  If maxRows
is larger than the number of rows in `zep` then a UnitRange with all rows is returned.
"""
randomsubset(zep::Zeppelin, maxRows::Int) =
    return maxRows<size(zep.data,1) ? #
        sort(Random.shuffle(collect(eachparticle(zep)))[1:maxRows]) :
        eachparticle(zep)

"""
    residual(zep::Zeppelin, row::Int, withImgs=false)::Union{Spectrum,Missing}

Returns the residual spectrum for the particle at `row` or missing if it does not exist. Never returns images.
"""
function residual(zep::Zeppelin, row::Int, withImgs=false)::Union{Spectrum,Missing}
    res = missing
    try
        filename = spectrumfilename(zep, row, "Residual", ".msa")
        res = isfile(filename) ? loadspectrum(filename) : missing
    catch
        # Ignore errors, just return missing
    end
    return res
end

"""
    maxparticle(zep::Zeppelin, rows::Union{AbstractVector{Int},UnitRange{Int}})
    maxresidual(zep::Zeppelin, rows::Union{AbstractVector{Int},UnitRange{Int}})

Computes the maxparticle spectrum from the specified particles by `rows`.
"""
function maxparticle(zep::Zeppelin, rows::Union{AbstractVector{Int},UnitRange{Int}})
    mp, firstspec = missing, missing
    for spec in filter(s->!ismissing(s), map(r->spectrum(zep,r,false),rows))
        @assert spec isa Spectrum
        firstspec = ismissing(firstspec) ? spec : firstspec
        cx = NeXLSpectrum.counts(spec)
        mp = ismissing(mp) ? cx : max.(mp, cx)
    end
    # Now copy it out as a spectrum
    props = copy(firstspec.properties)
    props[:Name] = "MaxParticle"
    return Spectrum(firstspec.energy, mp, props)
end

"""
    maxresidual(zep::Zeppelin, rows::Union{AbstractVector{Int},UnitRange{Int}}=1:1000000)

Computes the max-residual spectrum from the specified particles by `rows`.  The residual must have been
previously computed in quantify(zep, ....)
"""
function maxresidual(zep::Zeppelin, rows::Union{AbstractVector{Int},UnitRange{Int}})
    mp, firstspec = missing, missing
    for spec in filter(s->!ismissing(s), map(r->residual(zep,r,false),rows))
        @assert spec isa Spectrum
        firstspec = ismissing(firstspec) ? spec : firstspec
        cx = NeXLSpectrum.counts(spec)
        mp = ismissing(mp) ? cx : max.(mp, cx)
    end
    # Now copy it out as a spectrum
    props = copy(firstspec.properties)
    props[:Name] = "MaxResidual"
    return Spectrum(firstspec.energy, mp, props)
end

"""
    rowsmax(zep::Zeppelin, col::Symbol; n::Int=10000000, classname::Union{Missing,AbstractString}=missing)

Returns the row indices associated with the `n` maximum values in the `col` column.

Examples:

    # Plot ten largest by :DAVG
    plot(zep, rowsmax(zep, :DAVG, 10))
    # data from 100 most Iron-rich sorted by DAVG
    rowsmax(DataFrame, zep, rowsmax(zep, :FE, 100), sortCol=:DAVG)
"""
function rowsmax(zep::Zeppelin, col::Symbol; n::Int=10000000, classname::Union{Missing,AbstractString}=missing)
    if ismissing(classname)
        d = collect(zip(1:250,zep.data[:,col]))
        s = sort(d, lt=(a1,a2)->isless(a1[2],a2[2]),rev=true)
        return getindex.(s[intersect(eachparticle(zep),1:n)],1)
    else
        cr = filter(r->zep.data[r,:CLASS]==classname, eachparticle(zep))
        d = collect(zip(cr, zep.data[cr,col]))
        s = sort(d, lt=(a1,a2)->isless(a1[2],a2[2]),rev=true)
        return getindex.(s[intersect(eachindex(cr),1:n)],1)
    end
end
"""
    rowsmin(zep::Zeppelin, col::Symbol; n::Int=10000000, classname::Union{Missing,AbstractString}=missing)

Returns the row indices associated with the `n` minimum values in the `col` column.

Examples:

    # Plot ten smallest by :DAVG
    plot(zep, rowsmin(zep, :DAVG, 10))
    # data from 100 most Iron-rich sorted by DAVG
    rowsmin(DataFrame, zep, rowsmax(zep, :FE, 100), sortCol=:DAVG)
"""
function rowsmin(zep::Zeppelin, col::Symbol; n::Int=20, classname::Union{Missing,AbstractString}=missing)
    if ismissing(classname)
        d = collect(zip(eachparticle(zep),zep.data[:,col]))
        s = sort(d, lt=(a1,a2)->isless(a1[2],a2[2]),rev=false)
        return getindex.(s[intersect(eachparticle(zep),1:n)],1)
    else
        cr = filter(r->zep.data[r,:CLASS]==classname, eachparticle(zep))
        d = collect(zip(cr, zep.data[cr,col]))
        s = sort(d, lt=(a1,a2)->isless(a1[2],a2[2]),rev=false)
        return getindex.(s[intersect(eachindex(cr),1:n)],1)
    end
end

"""
    rowsclass(zep::Zeppelin, classname::AbstractString; shuffle=false, n=1000000)

Returns the row indices associated with the specified `classname`.  If `shuffle` is true,
the row indices are shuffled so that a randomized n can be plucked off to plot or other.

Example:

    # Plot a randomized selection of 10 "Calcite" particles
    plot(zep, rowsclass(zep, "Calcite", shuffle=true, n=10), xmax=8.0e3)
"""
function rowsclass(zep::Zeppelin, classname::AbstractString; shuffle=false, n=1000000)
    res = filter(i->zep.data[i,:CLASS]==classname, eachparticle(zep))
    return (shuffle ? Random.shuffle(res) : res)[intersect(eachindex(res),1:n)]
end

rowsclass(zep::Zeppelin, classnames::AbstractVector{<:AbstractString}; shuffle=false, n=1000000) =
    mapreduce(cn->rowsclass(zep, cn, shuffle=shuffle, n=n), union, classnames)

"""
    Base.filter(filt::Function, zep::Zeppelin)

Use a function of the form `filt(row)::Bool` to filter a Zeppelin dataset returning a new Zeppelin dataset
with only the rows for which the function evaluated true.

Example:

    gsr = filter(row->startswith( String(row[:CLASS]), "GSR."), zep)

"""
Base.filter(filt::Function, zep::Zeppelin) = filter(r->filt(zep.data[r,:]), eachparticle(zep))

const MORPH_COLS = ( "NUMBER", "XABS", "YABS", "DAVG", "DMIN", "DMAX", "DPERP", "PERIMETER", "ORIENTATION", "AREA" )
const CLASS_COLS = ( "CLASS", "VERIFIEDCLASS", "IMPORTANCE" )
const COMP_COLS = ( "FIRSTELM", "FIRSTPCT", "SECONDELM", "SECONDPCT", "THIRDELM", "THIRDPCT", "FOURTHELM", "FOURTHPCT", "COUNTS" )
const ALL_ELMS = map(elm->uppercase(elm.symbol), elements[1:95])

const ALL_COMPOSITIONAL_COLUMNS = ( "FIRST", "FIRSTELM", "SECONDELM", "THIRDELM", "FOURTHELM", "FIRSTELM", "SECONDELM",
    "THIRDELM", "FOURTHELM", "COUNTS1", "COUNTS2", "COUNTS3", "COUNTS4", "FIRSTPCT", "SECONDPCT", "THIRDPCT",
    "FOURTHPCT", "TYPE4ET", "COUNTS", "FITQUAL", "COMPHASH" )
const ALL_CLASS_COLS = ( "CLASS", "VERIFIEDCLASS", "IMPORTANCE" )

"""
    translate(zep::Zeppelin, am::AbstractAffineMap)::Zeppelin

The result is a copy of `zep` with the columns `:XABS` and `:YABS` replaced with the DataFtranslated coordinates.
"""
function translate(zep::Zeppelin, am::AbstractAffineMap)::Zeppelin
    res=copy(zep)
    foreach(eachrow(res)) do r
        ( r[:XABS], r[:YABS] ) = am(SVector(r[:XABS], r[:YABS]))
    end
    return res
end


text:
import Base.+
function +(lhs::AFArray, rhs::AFArray)
  out = af_array[0]
  af_add(out, lhs, rhs, true)
  AFArray(out[1])
end
+{T<:Number}(lhs::AFArray, rhs::T) = lhs + AFArray(T[rhs])
+{T<:Number}(lhs::T, rhs::AFArray) = AFArray(T[lhs]) + rhs

import Base.-
function -(lhs::AFArray, rhs::AFArray)
  out = af_array[0]
  af_sub(out, lhs, rhs, true)
  AFArray(out[1])
end
-(in::AFArray) = zeros(in) - in
-{T<:Number}(lhs::AFArray, rhs::T) = lhs - AFArray(T[rhs])
-{T<:Number}(lhs::T, rhs::AFArray) = AFArray(T[lhs]) - rhs

import Base.(.*)
function .*(lhs::AFArray, rhs::AFArray)
  out = af_array[0]
  af_mul(out, lhs, rhs, true)
  AFArray(out[1])
end

import Base.*
*{T<:Number}(lhs::AFArray, rhs::T) = lhs .* AFArray(T[rhs])
*{T<:Number}(lhs::T, rhs::AFArray) = AFArray(T[lhs]) .* rhs

function Base.dot(lhs::AFArray, rhs::AFArray)
  out = af_array[0]
  af_dot(out, lhs, rhs, AF_MAT_NONE, AF_MAT_NONE)
  AFArray(out[1])
end

function matmul(lhs::AFArray, rhs::AFArray, optlhs, optrhs)
  out = af_array[0]
  af_matmul(out, lhs, rhs, optlhs, optrhs)
  AFArray(out[1])
end
*(lhs::AFArray, rhs::AFArray) = matmul(lhs, rhs, AF_MAT_NONE, AF_MAT_NONE)
Base.A_mul_Bt(lhs::AFArray, rhs::AFArray) = matmul(lhs, rhs, AF_MAT_NONE, AF_MAT_TRANS)
Base.At_mul_B(lhs::AFArray, rhs::AFArray) = matmul(lhs, rhs, AF_MAT_TRANS, AF_MAT_NONE)
Base.At_mul_Bt(lhs::AFArray, rhs::AFArray) = matmul(lhs, rhs, AF_MAT_TRANS, AF_MAT_TRANS)

for (fname, afname) in [(:exp, :af_exp),
                        (:expm1, :af_expm1),
                        (:log, :af_log),
                        (:log10, :af_log10),
                        (:log1p, :af_log1p),
                        (:sqrt, :af_sqrt),
                        (:tanh, :af_tanh),
                        (:abs, :af_abs),
                        (:ceil, :af_ceil),
                        (:floor, :af_floor),
                        (:round, :af_round),
                        (:sign, :af_sign),
                        (:trunc, :af_trunc)]
  @eval begin
    function Base.$(fname)(in::AFArray)
      out = af_array[0]
      $afname(out, in)
      AFArray(out[1])
    end
  end
end

import Base: >=, >, <=, <
for (fname, afname) in [(:>=, :af_ge),
                        (:>, :af_gt),
                        (:<=, :af_le),
                        (:<, :af_lt)]
  @eval begin
    function $(fname)(lhs::AFArray, rhs::AFArray)
      out = af_array[0]
      $afname(out, lhs, rhs, true)
      AFArray(out[1])
    end
    $(fname)(lhs::AFArray, rhs::Number) = $(fname)(lhs, AFArray([rhs]))
    $(fname)(lhs::Number, rhs::AFArray) = $(fname)(AFArray([lhs]), rhs)
  end
end


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

function hid_init()
    ccall((:hid_init, hidapi), Cint, ())
end

function hid_exit()
    ccall((:hid_exit, hidapi), Cint, ())
end

function hid_enumerate(vendor_id, product_id)
    ccall((:hid_enumerate, hidapi), Ptr{hid_device_info}, (UInt16, UInt16), vendor_id, product_id)
end

function hid_free_enumeration(devs)
    ccall((:hid_free_enumeration, hidapi), Cvoid, (Ptr{hid_device_info},), devs)
end

function hid_open(vendor_id, product_id, serial_number)
    ccall((:hid_open, hidapi), Ptr{hid_device}, (UInt16, UInt16, Ptr{Cwchar_t}), vendor_id, product_id, serial_number)
end

function hid_open_path(path)
    ccall((:hid_open_path, hidapi), Ptr{hid_device}, (Cstring,), path)
end

function hid_write(dev, data, length)
    ccall((:hid_write, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t), dev, data, length)
end

function hid_read_timeout(dev, data, length, milliseconds)
    ccall((:hid_read_timeout, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t, Cint), dev, data, length, milliseconds)
end

function hid_read(dev, data, length)
    ccall((:hid_read, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t), dev, data, length)
end

function hid_set_nonblocking(dev, nonblock)
    ccall((:hid_set_nonblocking, hidapi), Cint, (Ptr{hid_device}, Cint), dev, nonblock)
end

function hid_send_feature_report(dev, data, length)
    ccall((:hid_send_feature_report, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t), dev, data, length)
end

function hid_get_feature_report(dev, data, length)
    ccall((:hid_get_feature_report, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t), dev, data, length)
end

function hid_get_input_report(dev, data, length)
    ccall((:hid_get_input_report, hidapi), Cint, (Ptr{hid_device}, Ptr{Cuchar}, Csize_t), dev, data, length)
end

function hid_close(dev)
    ccall((:hid_close, hidapi), Cvoid, (Ptr{hid_device},), dev)
end

function hid_get_manufacturer_string(dev, string, maxlen)
    ccall((:hid_get_manufacturer_string, hidapi), Cint, (Ptr{hid_device}, Ptr{Cwchar_t}, Csize_t), dev, string, maxlen)
end

function hid_get_product_string(dev, string, maxlen)
    ccall((:hid_get_product_string, hidapi), Cint, (Ptr{hid_device}, Ptr{Cwchar_t}, Csize_t), dev, string, maxlen)
end

function hid_get_serial_number_string(dev, string, maxlen)
    ccall((:hid_get_serial_number_string, hidapi), Cint, (Ptr{hid_device}, Ptr{Cwchar_t}, Csize_t), dev, string, maxlen)
end

function hid_get_indexed_string(dev, string_index, string, maxlen)
    ccall((:hid_get_indexed_string, hidapi), Cint, (Ptr{hid_device}, Cint, Ptr{Cwchar_t}, Csize_t), dev, string_index, string, maxlen)
end

function hid_error(dev)
    ccall((:hid_error, hidapi), Ptr{Cwchar_t}, (Ptr{hid_device},), dev)
end

function hid_version()
    ccall((:hid_version, hidapi), Ptr{hid_api_version}, ())
end

function hid_version_str()
    ccall((:hid_version_str, hidapi), Cstring, ())
end


text:
function train(m::Model, dtrn::Data; epochs=1, showprogress=true,
               issave=true, optimize=false)
    optimize && hyperoptimization(m, dtrn; showloss=showprogress)
    build(m, dtrn)
    dtrn = minibatch(dtrn.x, dtrn.y, m.params["batchsize"]; shuffle=true)

    showprogress ? progress!(adam(m.model, repeat(
        dtrn, epochs); lr=m.params["lr"])) : adam!(
            m.model, repeat(dtrn, epochs); lr=m.params["lr"])

    issave && savemodel(m)
    m, dtrn
end

function train(m::Model, traindata::Dict{String, Array{T,1} where T};
               epochs=1, showprogress=true, issave=true, optimize=false)
    dtrn = getbatches(m, traindata; batchsize=m.params["batchsize"],
                      showtime=false)
    train(m, dtrn; epochs=epochs, showprogress=showprogress, issave=issave,
          optimize=optimize)
end

function train(m::Model, trainpath::String; epochs=1, showprogress=true,
               issave=true, optimize=true)
    m.datapath = dirname(trainpath)
    extension = splitext(trainpath)[2]
    traindata = in(extension, (".csv", ".tsv")) ? csv2data(trainpath) :
        readtext(trainpath)
    train(m, traindata; epochs=epochs, showprogress=showprogress, issave=issave,
          optimize=optimize)
end

function partialtrain(m::Model, trainpath::String)
    m.datapath = dirname(trainpath)
    traindata = csv2data(trainpath)
    partialtrain(m, traindata)
end

function partialtrain(m::Model, traindata::Dict{String, Array{T,1} where T})
    dtrn = getbatches(m, traindata; batchsize=m.params["batchsize"])
    partialtrain(m, traindata)
end

function partialtrain(m::Model, dtrn::Data)
    train(m, dtrn; epochs=1)
end

function crossvalidate(m::Model, dtrn::Data; k=5, epochs=1, showprogress=false)
    indices = sample(1:size(dtrn.x, 2), size(dtrn.x, 2), replace=false)
    dx, dy = dtrn.x[:, indices], dtrn.y[:, indices]
    xfolds, yfolds = kfolds(dx, k), kfolds(dy, k)
    trainloss, testloss = zeros(2)

    for i=1:k
        foldxtrn = hcat([xfolds[j] for j=1:k if j != i]...)
        foldytrn = hcat([yfolds[j] for j=1:k if j != i]...)

        foldxtst = xfolds[i]
        foldytst = yfolds[i]

        dftrn = minibatch(foldxtrn, foldytrn, m.params["batchsize"];
                          shuffle=true)
        dftst = minibatch(foldxtst, foldytst, m.params["batchsize"];
                          shuffle=true)

        showprogress ? progress!(adam(m.model, repeat(dftrn, epochs))) :
            adam!(m.model, repeat(dftrn, epochs))

        trainloss += sum([m.model(x, y) for (x, y) in dftrn])
        testloss += sum([m.model(x, y) for (x, y) in dftst])
    end

    trainloss /= k
    testloss /= k

    showprogress && println("Train loss: $trainloss")
    showprogress && println("Test loss: $testloss")

    return testloss
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
module SandDance

using Electron, DataValues

import IteratorInterfaceExtensions, TableTraits, IterableTables, JSON, FilePaths

export SandDanceWindow

app = nothing

mutable struct SandDanceWindow
    w::Window

    function SandDanceWindow(data)
        TableTraits.isiterabletable(data) === false && error("'data' is not a table.")

        main_html_path = joinpath(FilePaths.@__PATH__, "..", "sanddanceapp", "build", "index.html")

        global app

        if app == nothing
            app = Application()
        end

        w = Window(app, main_html_path, options=Dict("title" => "SandDance"))

        it = IteratorInterfaceExtensions.getiterator(data)

        data_dict = [Dict{Symbol,Any}(c[1] => isa(c[2], DataValue) ? (isna(c[2]) ? nothing : get(c[2])) : c[2] for c in zip(keys(r), values(r))) for r in it]

        data = JSON.json(data_dict)

        code = "global.startApp($data)"

        run(w, code)

        new(w)
    end
end

end # module


text:
module GMP

if isfile(joinpath(Pkg.dir(string(current_module())),"deps","path.jl"))
    include("../deps/path.jl")
else
    error("GMP not properly installed. Please run Pkg.build(\"GMP\")")
end

include = joinpath(prefix,"include")
lib = joinpath(prefix,"lib")

end # module


text:
##############################################################################
##
## Object constructed by the user
##
##############################################################################
function ife(id::Term, time::Term, rank::Int)
    InteractiveFixedEffectTerm(Symbol(id), Symbol(time), rank)
end

# Object constructed by the user
struct InteractiveFixedEffectTerm <: AbstractTerm
    id::Symbol
    time::Symbol
    rank::Int
end

has_ife(x::InteractiveFixedEffectTerm) = true
has_ife(x::FunctionTerm{typeof(ife)}) = true
has_ife(x::AbstractTerm) = false
function parse_interactivefixedeffect(df::AbstractDataFrame, formula::FormulaTerm)
    m = nothing
    for term in FixedEffectModels.eachterm(formula.rhs)
        if term isa FunctionTerm{typeof(ife)}
            m = InteractiveFixedEffectTerm(term.args_parsed[1].sym, term.args_parsed[2].sym, term.args_parsed[3].n)
        elseif term isa InteractiveFixedEffectTerm
            m = term
        end
    end
    return m, FormulaTerm(formula.lhs, tuple((term for term in FixedEffectModels.eachterm(formula.rhs) if !has_ife(term))...))
end

abstract type AbstractFactorModel{T} end
abstract type AbstractFactorSolution{T} end

# to deprecate
function OldInteractiveFixedEffectFormula(arg)
    arg.head == :tuple && length(arg.args) == 2 || throw("ife does not have a correct syntax")
    arg1 = arg.args[1]
    arg2 = arg.args[2]
    arg1.head == :call && arg1.args[1] == :+ && length(arg1.args) == 3 || throw("ife does not have a correct syntax")
    InteractiveFixedEffectTerm(arg1.args[2], arg1.args[3], arg2)
end

##############################################################################
##
## Factor Model
##
##############################################################################

struct FactorModel{Rank, W, Rid, Rtime} <: AbstractFactorModel{Rank}
    y::Vector{Float64}
    sqrtw::W
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
end

function FactorModel(y::Vector{Float64}, sqrtw::W, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, rank::Int) where {W, Rid, Rtime}
    FactorModel{rank, W, Rid, Rtime}(y, sqrtw, idrefs, timerefs)
end

LinearAlgebra.rank(::FactorModel{Rank}) where {Rank} = Rank

struct FactorSolution{Rank, Tid, Ttime} <: AbstractFactorSolution{Rank}
    idpool::Tid
    timepool::Ttime
end

function FactorSolution(idpool::Tid, timepool::Ttime) where {Tid, Ttime}
    r = size(idpool, 2)
    @assert r == size(timepool, 2)
    FactorSolution{r, Tid, Ttime}(idpool, timepool)
end

function Base.view(f::AbstractFactorSolution, I::Union{AbstractArray,Colon,Int64}...)
    FactorSolution(view(f.idpool, I...), view(f.timepool, I...))
end

## subtract_factor! and subtract_b!
function subtract_factor!(fm::AbstractFactorModel, fs::AbstractFactorSolution)
    for r in 1:rank(fm)
        subtract_factor!(fm, view(fs, :, r))
    end
end

function subtract_factor!(fm::AbstractFactorModel, fs::FactorSolution{1})
    @inbounds @simd for i in 1:length(fm.y)
        fm.y[i] -= fm.sqrtw[i] * fs.idpool[fm.idrefs[i]] * fs.timepool[fm.timerefs[i]]
    end
end

## rescale a factor model
function reverse(m::Matrix{R}) where {R}
    out = similar(m)
    for j in 1:size(m, 2)
        invj = size(m, 2) + 1 - j
        @inbounds @simd for i in 1:size(m, 1)
            out[i, j] = m[i, invj]
        end
    end
    return out
end
function rescale!(fs::FactorSolution{1})
    out = norm(fs.timepool)
    rmul!(fs.idpool, out)
    rmul!(fs.timepool, 1 / out)
end
# normalize factors and loadings so that F'F = Id, Lambda'Lambda diagonal
function rescale!(newfs::AbstractFactorSolution, fs::AbstractFactorSolution)
    U = eigen!(Symmetric(fs.timepool' * fs.timepool))
    sqrtDx = diagm(0 => sqrt.(abs.(U.values)))
    mul!(newfs.idpool,  fs.idpool,  U.vectors * sqrtDx)
    V = eigen!(Symmetric(newfs.idpool' * newfs.idpool))
    mul!(newfs.idpool, fs.idpool, reverse(U.vectors * sqrtDx * V.vectors))
    mul!(newfs.timepool, fs.timepool, reverse(U.vectors * (sqrtDx \ V.vectors)))
    return newfs
end

rescale(fs::FactorSolution) = rescale!(similar(fs), fs)

## Create dataframe from pooledfactors
function getfactors(fp::AbstractFactorModel, fs::AbstractFactorSolution)
    # partial out Y and X with respect to i.id x factors and i.time x loadings
    newfes = FixedEffect[]
    for r in 1:rank(fp)
        idfe = FixedEffect{typeof(fp.idrefs), Vector{Float64}}(fp.idrefs, fs.timepool[fp.timerefs, r], length(fs.idpool))
        push!(newfes, idfe)
        timefe = FixedEffect{typeof(fp.timerefs), Vector{Float64}}(fp.timerefs, fs.idpool[fp.idrefs, r], length(fs.timepool))
        push!(newfes, timefe)
    end
    # obtain the residuals and cross
    return newfes
end

function DataFrames.DataFrame(fp::AbstractFactorModel, fs::AbstractFactorSolution, esample::AbstractVector{Bool})
    df = DataFrame()
    for r in 1:rank(fp)
        # loadings
        df[!, Symbol("loadings$r")] = Vector{Union{Float64, Missing}}(missing, length(esample))
        df[esample, Symbol("loadings$r")] = fs.idpool[:, r][fp.idrefs]

        df[!, Symbol("factors$r")] = Vector{Union{Float64, Missing}}(missing, length(esample))
        df[esample, Symbol("factors$r")] = fs.timepool[:, r][fp.timerefs]
    end
    return df
end

##############################################################################
##
## Interactive Fixed Effect Models
##
##############################################################################

struct InteractiveFixedEffectsModel{Rank, W, Rid, Rtime} <: AbstractFactorModel{Rank}
    y::Vector{Float64}
    sqrtw::W
    X::Matrix{Float64}
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
end

function InteractiveFixedEffectsModel(y::Vector{Float64}, sqrtw::W, X::Matrix{Float64}, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, rank::Int) where {W, Rid, Rtime}
    InteractiveFixedEffectsModel{rank, W, Rid, Rtime}(y, sqrtw, X, idrefs, timerefs)
end

LinearAlgebra.rank(::InteractiveFixedEffectsModel{Rank}) where {Rank} = Rank

function Base.convert(::Type{FactorModel}, f::InteractiveFixedEffectsModel{Rank, W, Rid, Rtime}) where {Rank, W, Rid, Rtime}
    FactorModel{Rank, W, Rid, Rtime}(f.y, f.sqrtw, f.idrefs, f.timerefs)
end

struct InteractiveFixedEffectsSolution{Rank, Tb, Tid, Ttime} <: AbstractFactorSolution{Rank}
    b::Tb
    idpool::Tid
    timepool::Ttime
end
function InteractiveFixedEffectsSolution(b::Tb, idpool::Tid, timepool::Ttime) where {Tb, Tid, Ttime}
    r = size(idpool, 2)
    r == size(timepool, 2) || throw("factors and loadings don't have same dimension")
    InteractiveFixedEffectsSolution{r, Tb, Tid, Ttime}(b, idpool, timepool)
end
Base.convert(::Type{FactorSolution}, f::InteractiveFixedEffectsSolution) = FactorSolution(f.idpool, f.timepool)

struct InteractiveFixedEffectsSolutionT{Rank, Tb, Tid, Ttime} <: AbstractFactorSolution{Rank}
    b::Tb
    idpool::Tid
    timepool::Ttime
end
function InteractiveFixedEffectsSolutionT(b::Tb, idpool::Tid, timepool::Ttime) where {Tb, Tid, Ttime}
    r = size(idpool, 1)
    r == size(timepool, 1) || throw("factors and loadings don't have same dimension")
    InteractiveFixedEffectsSolutionT{r, Tb, Tid, Ttime}(b, idpool, timepool)
end

function rescale(fs::InteractiveFixedEffectsSolution)
    fss = FactorSolution(fs.idpool, fs.timepool)
    newfss = similar(fss)
    rescale!(newfss, fss)
    InteractiveFixedEffectsSolution(fs.b, newfss.idpool, newfss.timepool)
end

struct HalfInteractiveFixedEffectsModel{Rank, W, Rid, Rtime} <: AbstractFactorModel{Rank}
    y::Vector{Float64}
    sqrtw::W
    X::Matrix{Float64}
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
    timepool::Matrix{Float64}
    size::Tuple{Int, Int}
end

function HalfInteractiveFixedEffectsModel(y::Vector{Float64}, sqrtw::W, X::Matrix{Float64}, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, timepool::Matrix{Float64}, size, rank::Int) where {W, Rid, Rtime}
    HalfInteractiveFixedEffectsModel{rank, W, Rid, Rtime}(y, sqrtw, X, idrefs, timerefs, timepool, size)
end

struct HalfInteractiveFixedEffectsSolution{Rank, Tb, Tid} <: AbstractFactorSolution{Rank}
    b::Tb
    idpool::Tid
end

##############################################################################
##
## Results
##
##############################################################################'

struct FactorResult
    esample::BitVector
    augmentdf::DataFrame

    rss::Float64
    iterations::Int64
    converged::Bool
end

# result
struct InteractiveFixedEffectModel <: RegressionModel
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix
    vcov_type::CovarianceEstimator

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::FormulaTerm        # Original formula

    nobs::Int64             # Number of observations
    dof_residual::Int64      # degree of freedoms

    r2::Float64             # R squared
    adjr2::Float64           # R squared adjusted
    r2_within::Float64      # R within

    rss::Float64
    iterations::Int         # Number of iterations
    converged::Bool         # Has the demeaning algorithm converged?

end
StatsBase.coef(x::InteractiveFixedEffectModel) = x.coef
StatsBase.coefnames(x::InteractiveFixedEffectModel) = x.coefnames
StatsBase.vcov(x::InteractiveFixedEffectModel) = x.vcov
StatsBase.nobs(x::InteractiveFixedEffectModel) = x.nobs
StatsBase.dof_residual(x::InteractiveFixedEffectModel) = x.dof_residual
StatsBase.r2(x::InteractiveFixedEffectModel) = x.r2
StatsBase.adjr2(x::InteractiveFixedEffectModel) = x.adjr2
StatsBase.islinear(x::InteractiveFixedEffectModel) = false
StatsBase.rss(x::InteractiveFixedEffectModel) = x.rss
StatsBase.predict(::InteractiveFixedEffectModel, ::AbstractDataFrame) = error("predict is not defined for linear factor models. Use the option save = true")
StatsBase.residuals(::InteractiveFixedEffectModel, ::AbstractDataFrame) = error("residuals is not defined for linear factor models. Use the option save = true")
function StatsBase.confint(x::InteractiveFixedEffectModel)
    scale = quantile(TDist(dof_residual(x)), 1 - (1-0.95)/2)
    se = stderror(x)
    hcat(x.coef -  scale * se, x.coef + scale * se)
end

title(::InteractiveFixedEffectModel) = "Interactive Fixed Effect Model"
top(x::InteractiveFixedEffectModel) = [
            "Number of obs" sprint(show, nobs(x); context=:compact => true);
            "Degree of freedom" sprint(show, nobs(x) - dof_residual(x); context=:compact => true);
            "R2"  @sprintf("%.3f", x.r2);
            "R2 within"  @sprintf("%.3f", x.r2_within);
            "Iterations" sprint(show, x.iterations; context=:compact => true);
            "Converged" sprint(show, x.converged; context=:compact => true)
            ]

function Base.show(io::IO, x::InteractiveFixedEffectModel)
    show(io, coeftable(x))
end
function StatsBase.coeftable(x::InteractiveFixedEffectModel)
    ctitle = title(x)
    ctop = top(x)
    cc = coef(x)
    se = stderror(x)
    coefnms = coefnames(x)
    conf_int = confint(x)
    # put (intercept) last
    if !isempty(coefnms) && ((coefnms[1] == Symbol("(Intercept)")) || (coefnms[1] == "(Intercept)"))
        newindex = vcat(2:length(cc), 1)
        cc = cc[newindex]
        se = se[newindex]
        conf_int = conf_int[newindex, :]
        coefnms = coefnms[newindex]
    end
    tt = cc ./ se
    FixedEffectModels.CoefTable2(
        hcat(cc, se, tt, ccdf.(Ref(FDist(1, dof_residual(x))), abs2.(tt)), conf_int[:, 1:2]),
        ["Estimate","Std.Error","t value", "Pr(>|t|)", "Lower 95%", "Upper 95%" ],
        ["$(coefnms[i])" for i = 1:length(cc)], 4, ctitle, ctop)
end


text:
using JLD2;
using Printf;
using StatsPlots;
using LaTeXStrings;

include("runit.jl"); # for types
include("../experiment_helpers.jl");
include("../thresholds.jl");

name = "experiment_bai2";

@load "$name.dat" dist μ pep srs data δs βs repeats seed

dump_stats(pep, μ, δs, βs, srs, data, repeats);

for i in 1:length(δs)
    plot(boxes(pep, μ, δs[i], βs[i], srs, getindex.(data, i)));
    ylabel!(L"\texttt{Empirical stopping time}")
    savefig("$(name)_$i.pdf");
end


text:
function hello()
    # println("Hello, World!")
    return "Hello, World!"
end

max_stars_repo_name: mauricioriva/learn-julia

max_issues_repo_name: mauricioriva/learn-julia

max_forks_repo_name: mauricioriva/learn-julia


text:
using PyEnum

abstract type AbstractColors <: AbstractIntEnum end
abstract type AbstractPermissions <: AbstractIntFlag end
@pyenum Colors::Int64 begin
    RED = 0
    GREEN = 1
    BLUE = 2
end
@pyenum Permissions::Int64 begin
    R = 1
    W = 2
    X = 16
end
function show()
    color_map = Dict(Colors.RED => "red", Colors.GREEN => "green", Colors.BLUE => "blue")
    a = Colors.GREEN
    if a == Colors.GREEN
        println("green")
    else
        println("Not green")
    end
    b = Permissions.R
    if b == Permissions.R
        println("R")
    else
        println("Not R")
    end
    @assert(length(color_map) != 0)
end

if abspath(PROGRAM_FILE) == @__FILE__
    show()
end


text:
using PyPlot
using Statistics
include("../../../src/extract_planet.jl")
include("../../../src/laplace_wisdom.jl")

function chop_coeff_inner(alpha,j)
# Computes f_1^(j) from equation (10) in Deck & Agol (2015).
beta = j*(1-alpha^1.5)
# Equation (11) in Deck & Agol (2015):
f1 = 2*beta*laplace_wisdom(1//2,j,1,alpha)+
    j*(3+beta^2)*laplace_wisdom(1//2,j,0,alpha)
if j == 1
  f1 -= alpha*(beta^2+2*beta+3)
end
f1 *= alpha/(beta^2*(1-beta^2))
return f1
end

function chop_coeff_outer(alpha,j)
# Computes f_2^(j) from equation (14) in Deck & Agol (2015).
kappa = j*(1/alpha^1.5-1)
# Equation (15) in Deck & Agol (2015):
f2 = (j*(kappa^2+3)+2*kappa)*laplace_wisdom(1//2,j,0,alpha)+
    2*kappa*laplace_wisdom(1//2,j,1,alpha)
if j == 1
  f2 -= (kappa^2-2*kappa+3)/alpha^2
end
f2 /= (kappa^2*(kappa^2-1))
return f2
end

function plot_ttv_diff(data,elements,tt1,count1,nplanet,jmax)

# Now make some plots:
#fig,axes = subplots(4,2,sharex="col",constrained_layout="True")
fig,axes = subplots(7,1,sharex="col",constrained_layout="True",figsize=(8,8))
#fig,axes = subplots(7,1,sharex="col",figsize=(24,24))
#fig,axes = subplots(7,1,sharex="col")
plabel = ["T1b","T1c","T1d","T1e","T1f","T1g","T1h"]
range = [1.5,1.5,1.5,3.0,12.0,8.0,4.0]

# First, go through the planets and compute their
# ephemerides:
coeff_planet = zeros(2,nplanet)
for ip=1:nplanet
  eobs,tobs,sobs,nobs = extract_planet(data,ip)
  fn = zeros(2,nobs); fn[1,:] .= 1.0; fn[2,:] .= eobs
  coeff,cov = regress(fn,tobs,sobs)
  # Now plot TTVFast:
  fn = zeros(Float64,2,count1[ip+1])
  sig = ones(count1[ip+1])
  tti1 = tt1[ip+1,1:count1[ip+1]]
  tt_ref1 = zeros(count1[ip+1])
  epoch = zeros(count1[ip+1])
  for j=1:count1[ip+1]
    fn[1,j] = 1.0
    fn[2,j] = round(Int64,(tti1[j]-elements[ip+1,3])/elements[ip+1,2])
    epoch[j] = round((tti1[j]-coeff[1])/coeff[2])
    tt_ref1[j] = coeff[1]+coeff[2]*epoch[j]
  end
  coeff,cov = regress(fn,tti1,sig)
  coeff_planet[:,ip] .= coeff
end

# Now make some plots:
eoffset = [43,10,75,9,7,3,22]
ord = [5,5,18,30,18,18,18]
for ip=1:nplanet
  ax = axes[ip]
  coeff = zeros(Float64,2)
  # Plot each "Mode"
  eobs,tobs,sobs,nobs = extract_planet(data,ip)
#  println("ip: ",ip," nobs: ",nobs)
  fn = zeros(2,nobs)
  fn[1,:] .= 1.0
  fn[2,:] .= eobs
  coeff,cov = regress(fn,tobs,sobs)
#  println(ip," ",coeff)
  # Now plot TTVFast:
  fn = zeros(Float64,2,count1[ip+1])
  sig = ones(count1[ip+1])
  tti1 = tt1[ip+1,1:count1[ip+1]]
  tt_ref1 = zeros(count1[ip+1])
  epoch = zeros(count1[ip+1])
  for j=1:count1[ip+1]
    fn[1,j] = 1.0
    fn[2,j] = round(Int64,(tti1[j]-elements[ip+1,3])/elements[ip+1,2])
    epoch[j] = round((tti1[j]-coeff[1])/coeff[2])
    tt_ref1[j] = coeff[1]+coeff[2]*epoch[j]
  end
  coeff,cov = regress(fn,tti1,sig)
#  println(ip," ",coeff)
  tt_ref1 = coeff[1] .+coeff[2] .*fn[2,:]
  ttv1 = (tti1 .-tt_ref1) .*(24*60)
#  ax.plot(tti1,ttv1,label=plabel[ip])
  # Now, fit a high-order polynomial to TTVs, and overplot it:
  fn = zeros(ord[ip]+1,count1[ip+1])
  fn[1,:] .= 1.0
  for j=1:ord[ip]
    fn[j+1,:] .= epoch.^j
  end
  coeff_poly,cov = regress(fn,tti1,sig)
  mod_poly = ones(count1[ip+1])*coeff_poly[1]
  for j=1:ord[ip]
     mod_poly .+= fn[j+1,:] .*coeff_poly[j+1]
  end
#  ax.plot(tti1,(mod_poly .-tt_ref1) .*(24*60))
  ttv_mod = (tti1 .- mod_poly) .*(24*60)
  ax.plot(tti1,ttv_mod,label=plabel[ip])
#  ttv_obs = tobs .- coeff_poly[1] .- coeff_poly[2] .*(eobs .+eoffset[ip])
  ttv_obs = tobs .- coeff_poly[1] .- coeff_poly[2] .*eobs
  for j=2:ord[ip]
#    ttv_obs .-= coeff[j+1] .*(eobs .+eoffset[ip]).^j
    ttv_obs .-= coeff_poly[j+1] .*eobs.^j
  end
  ttv_obs .*= 24*60
  sobs_min = sobs .*(24*60)
  if ip < 5
    igood = sobs_min .< 0.5*(maximum(ttv_mod)-minimum(ttv_mod))
  else
    igood = ones(Bool,size(ttv_obs)[1])
  end
  println("Fraction of points: ",sum(igood)/size(sobs)[1])
  ax.errorbar(tobs[igood], ttv_obs[igood], sobs_min[igood],fmt=".")
  if ip == 7
    ax.set_xlabel(L"BJD$_\mathrm{TDB}$-2,450,000")
  end
  if ip == 4
    ax.set_ylabel("TTV - Polynomial [min]")
  end
#  ax.legend(loc="upper right",fontsize=10)
  ax.legend(fontsize=8)
  ax.axis([7200,8900,-range[ip],range[ip]])
#  ax.axis([7960,8860,-range[ip],range[ip]])
  # Now, overplot chopping due to all adjacent planets:
  dt_chop = zeros(length(tti1))
  if ip > 1
    for jp=ip-1:-1:1
      # Compute the observed mean longitudes at the time of transit:
      lam1 = 2*pi*(tti1 .-coeff_planet[1,jp])./coeff_planet[2,jp]
      lam2 = 2*pi*(tti1 .-coeff_planet[1,ip])./coeff_planet[2,ip]
      dlambda = lam1-lam2
      alpha = (coeff_planet[2,jp]/coeff_planet[2,ip])^(2//3)
      for j=1:jmax
        dt_chop .+= elements[1+jp,1]*chop_coeff_outer(alpha,j)*sin.(j*dlambda)
      end
    end
  end
  if ip < nplanet
    for jp = ip+1:nplanet
      # Compute the observed mean longitudes at the time of transit:
      lam1 = 2*pi*(tti1 .-coeff_planet[1,ip])./coeff_planet[2,ip]
      lam2 = 2*pi*(tti1 .-coeff_planet[1,jp])./coeff_planet[2,jp]
      dlambda = lam1-lam2
      alpha = (coeff_planet[2,ip]/coeff_planet[2,jp])^(2//3)
      for j=1:jmax
        dt_chop .+= elements[1+jp,1]*chop_coeff_inner(alpha,j)*sin.(j*dlambda)
      end
    end
  end
  dt_chop .*= coeff_planet[2,ip]/(2pi)
  println(ip," range of chopping: ",minimum(dt_chop*(24*60))," ",maximum(dt_chop*(24*60)))
  # Remove a polynomial:
  coeff_poly,cov = regress(fn,dt_chop,sig)
  fill!(mod_poly,coeff_poly[1])
  for j=1:ord[ip]
     mod_poly .+= fn[j+1,:] .*coeff_poly[j+1]
  end
  ax.plot(tti1,(dt_chop-mod_poly)*(24*60),".",markersize=2.0)
#  ax.plot(tti1,(dt_chop-mod_poly)*(24*60) .- ttv_mod)
end
#ax = axes[8]
#ax.axis("off")
tight_layout()
subplots_adjust(hspace=0)
#println("Ephemerides: ",coeff_planet)
savefig("../T1_chopping.pdf",bbox_inches="tight")
return
end


text:
import JSON

# from https://github.com/JuliaLang/JSON.jl/issues/39
# but without `flatten` or `Nullable`
function deserialize(::Type{T}, json::AbstractString) where T
    deserialize(T, JSON.parse(json))
end

function deserialize(::Type{T}, json::Dict) where T
    T([deserialize(T, f, json) for f in fieldnames(T)]...)
end

function deserialize(::Type{Tf}, field::AbstractString, json::Dict) where Tf
    deserialize(Tf, json[field])
end

function deserialize(::Type{T}, field::Symbol, json::Dict) where T
    deserialize(fieldtype(T,field), string(field), json)
end

function deserialize(::Type{T}, i::Integer) where T<:Integer
    T(i)
end

function deserialize(::Type{T}, i::AbstractString) where T<:Integer
    parse(T, i)
end

function deserialize(::Type{T}, i::AbstractFloat) where T<:AbstractFloat
    T(i)
end

function deserialize(::Type{T}, i::AbstractString) where T<:AbstractFloat
    parse(T, i)
end

# my own extension to vectors
function deserialize(::Type{Vector{T}}, A::Vector{Any}) where T
    T[deserialize(T, a) for a in A]
end

function read_instance(prefix, key::String)
    open(joinpath(prefix, "$(key).instance.json")) do f
        PipeLayout.deserialize(Instance, readstring(f))
    end
end
read_instance(key::String) = read_instance(".", key)

function read_topology(prefix, key::String)
    open(joinpath(prefix, "$(key).topology.json")) do f
        PipeLayout.deserialize(Topology, readstring(f))
    end
end
read_topology(key::String) = read_topology(".", key)


text:
#
# regions.jl -
#
# Implment physical regions.
#

module Regions

export
    Region,
    boundingbox,
    grid2world,
    world2grid

import ...FourierOptics: center, recenter
using ..CoordinateTransforms
using ..Units

"""

A `Region` defines a rectangular sampled piece of physical plane perpendicular
to the direction of propagation.  The sampling step is the same in all
direction and the grid axes are aligned with the `(x,y)` axes of the world
coordinate system (axis `z` is assumed to be the direction of propagation).

A `Region` can be defined by:

```julia
Region(width, height, step=1mm)
```

where `width`, `height` and `step` are the physical dimensions and sampling
step (all lengths in SI units, that is meters) and which yields a region
centered at the origin.  The bounding box (see below) of the resulting region
will be at least of size `width` by `height`.  Another possibility is to use
two instances of `Range` to define the coordinates of the grid of samples along
the two axes:

```julia
Region(X, Y)
```

Finally, the following keyword-only constructor is available:

```julia
Region(; size=(nx,ny), step=stp, center=(xc,yc))
```

where `nx` and `ny` are the number of grid nodes along the dimensions of the
region (both are `100` by default), `stp` is the sampling step size (`1mm` by
default) and `(xc,yc)` are the coordinates of the geometrical center of the
region (`(0mm,0mm)` by default).

Because of rounding errors, the effective grid of samples embedded into the
resulting region may have slightly different coordinates.

A region can be translated or recentered:

```julia
R = Region(X, Y)          # define a region
Rt = R + (tx, ty)         # translate the region
Rc = recenter(R)          # recenter the region (0,0)
Rc = recenter(R, (x0,y0)) # recenter the region at (x0, y0)
```

A number of coordinate transformations can be performed:

```julia
R = Region(X, Y)              # define a region
grid2world(R, (i,j)) -> x, y  # get the world coordinates of a grid node
world2grid(R, (x,y)) -> i, j  # get the grid indices from world coordinates
```

Here, `(x,y)` denotes physical (*world*) coordinates and `(i,j)` denotes,
possibly fractional, grid indices.

Some information can be retrieved from the region instance:

```julia
R = Region(X, Y)  # define a region
length(R)         # yields the number of samples
size(R)           # yields the dimensions (in number of samples) of the region
size(R, i)        # yields the i-th dimension
step(R)           # yields the sampling step (in meters)
grid2world(R)     # yields the grid to world coordinate transformation
world2grid(R)     # yields the world to grid coordinate transformation
extrema(R)        # yields the extreme positions `(xmin,xmax,ymin,ymax)`
boundingbox(R)    # yields the bounding box of the region
center(R)         # yields the physical coordinates of the center of R
```

Compared to `extrema(R)` which yields the extreme world positions `(xmin, xmax,
ymin, ymax)` of the embedded grid nodes, `boundingbox(R)` accounts for an
additional margin of 1/2 sample on all sides.

"""
struct Region
    sz::NTuple{2,Int}        # number of samples in each direction
    i2w::CoordinateTransform # indices to world coordinate transform
    w2i::CoordinateTransform # world to indices coordinate transform
    Region(sz::NTuple{2,Integer}, i2w::CoordinateTransform) =
        (isfinite(step(i2w)) && step(i2w) != 0 ||
         throw(ArgumentError("invalid sampling step"));
         new(sz, i2w, inv(i2w)))
end
Base.size(R::Region) = R.sz
Base.size(R::Region, i) = R.sz[i]
Base.length(R::Region) = prod(size(R))
grid2world(R::Region) = R.i2w
grid2world(R::Region, args...) = R.i2w(args...)
world2grid(R::Region) = R.w2i
world2grid(R::Region, args...) = R.w2i(args...)
Base.step(R::Region) = step(grid2world(R))

Base.show(io::IO, ::MIME"text/plain", R::Region) = show(io, R)

function Base.show(io::IO, R::Region)
    print(io, "Region(size = $(size(R)), step = ")
    showlength(io, step(R))
    print(io, ", center = ")
    showposition(io, center(R))
    print(io, ")")
end

# Keyword-only constructor (to look like the output of the show method).
function Region(;
                size::NTuple{2,Integer} = (100,100),
                step::Real = 1mm,
                center::NTuple{2,Union{Real,Rational,Irrational}} = (0mm, 0mm))
    nx, ny = size
    xc, yc = center
    dx = step*(nx - 1)/2
    dy = step*(ny - 1)/2
    return Region(linspace(xc - dx, xc + dx, nx),
                  linspace(yc - dy, yc + dy, ny))
end

function Region(width::Real, height::Real, step::Real=1mm)
    # Build a large enough region centered at (0,0).
    n1 = ceil(Int, width/step)
    n2 = ceil(Int, height/step)
    x0 = -(n1 + 1)*step/2
    y0 = -(n2 + 1)*step/2
    return Region((n1, n2), CoordinateTransform(step, x0, y0))
end

function Region(X::Range{Tx}, Y::Range{Ty}) where {Tx<:Real,Ty<:Real}
    # Get the sampling step.
    xstp, ystp = step(X), step(Y)
    T = float(promote_type(Tx, Ty))
    if abs(xstp - ystp) > 4*eps(T)*max(abs(xstp), abs(ystp))
        throw(ArgumentError("steps must be (approximately) the same"))
    end
    stp = (xstp + ystp)/2

    # Get the number of samples along each side.
    nx, ny = length(X), length(Y)

    # Compute the offsets so as to keep the central position.
    x0 = ((first(X) - stp) + (last(X) - nx*stp))/2
    y0 = ((first(Y) - stp) + (last(Y) - ny*stp))/2

    # Build the region using the inner constructor.
    return Region((nx, ny), CoordinateTransform(stp, x0, y0))
end

"""
```julia
center(R) -> xc, yc
```

yields the physical coordinates `(xc,yc)` of the geometrical center of the
region `R`.

"""
center(R::Region) = ((n1, n2) = size(R);
                     grid2world(R, (1 + n1)/2, (1 + n2)/2))

"""
```julia
recenter(R, (x0, y0) = (0,0))
```

yields a new physical region of same physical size and sampling step as `R` but
whose center is at `(x0,y0)`.

"""
function recenter(R::Region, c::NTuple{2,Real} = (0,0))
    n1, n2 = size(R)
    stp = step(R)
    x0 = c[1] - (1 + n1)*stp/2
    y0 = c[2] - (1 + n2)*stp/2
    return Region((n1, n2), CoordinateTransform(stp, x0, y0))
end

recenter(R::Region, x0::Real, y0::Real) =
    recenter(R, (x0, y0))

# Translate the region.
Base.:+(R::Region, t::NTuple{2,Real}) =
    ((x0, y0) = origin(grid2world(R));
     Region(size(R), CoordinateTransform(step(R), x0 + t[1], y0 + t[2])))

Base.:-(R::Region, t::NTuple{2,Real}) =
    ((x0, y0) = origin(grid2world(R));
     Region(size(R), CoordinateTransform(step(R), x0 - t[1], y0 - t[2])))

"""
```julia
boundingbox(R) -> xmin, xmax, ymin, ymax
```

yields the physical size of the region `R`.  Compared to `extrema(R)`, an
extra margin of 1/2 sample is taken into account.

"""
function boundingbox(R::Region)
    n1, n2 = size(R)
    x0, y0 = grid2world(R, 0.5, 0.5)
    x1, y1 = grid2world(R, n1 + 0.5, n2 + 0.5)
    return min(x0, x1), max(x0, x1), min(y0, y1), max(y0, y1)
end

function Base.extrema(R::Region)
    x0, y0 = grid2world(R, 1, 1)
    x1, y1 = grid2world(R, size(R))
    return min(x0, x1), max(x0, x1), min(y0, y1), max(y0, y1)
end

end # module


text:
module BayesNets
__precompile__()

using Reexport
using Parameters

include(joinpath("CPDs", "cpds.jl"))
@reexport using BayesNets.CPDs
@reexport using BayesNets.CPDs.ProbabilisticGraphicalModels

import Base: *, /, +, -, normalize
import DataStructures: PriorityQueue, peek
import BayesNets.CPDs.ProbabilisticGraphicalModels: markov_blanket, is_independent, infer
import Iterators: subsets, product
import StatsBase: sample, Weights
import TikzGraphs: plot
import LightGraphs: DiGraph, add_edge!, rem_edge!,
       add_vertex!, rem_vertex!, has_edge,
       edges, topological_sort_by_dfs, in_neighbors,
       out_neighbors, neighbors, is_cyclic, nv, ne,
       outdegree, bfs_tree, dst

export
    BayesNet,
    DAG,

    parents,
    children,
    neighbors,
    descendants,
    has_edge,
    enforce_topological_order!,

    add_edge!,
    has_edge,

    rand_cpd,
    BayesNetSampler,
    gibbs_sample,
    GibbsSampler,

    # tables, formerly known as factors; the DataFrames based approach
    # TODO remove
    Table,
    table,
    sumout,
    normalize,
    # estimate_convergence,
    readxdsl,

    # generate BNs
    rand_discrete_bn,
    rand_bn_inference,
    get_sprinkler_bn,
    get_sat_fail_bn,
    get_asia_bn,

    # Factors
    Factor,
    pattern,
    reducedim!,

    # Sampling
    BayesNetSampler,
    DirectSampler,
    RejectionSampler,
    LikelihoodWeightedSampler,

    get_weighted_sample,
    get_weighted_sample!,

    # Inference
    InferenceState,
    InferenceMethod,
    ExactInference,
    LikelihoodWeightingInference,
    LoopyBelief,
    GibbsSamplingNodewise,
    GibbsSamplingFull,

    infer,

    DirichletPrior,
    UniformPrior,
    BDeuPrior,

    ScoringFunction,
    ScoreComponentCache,
    NegativeBayesianInformationCriterion,
    score_component,
    score_components,

    GraphSearchStrategy,
    K2GraphSearch,
    GreedyHillClimbing,
    GreedyThickThinning,
    ScanGreedyHillClimbing,

    statistics,
    index_data,
    adding_edge_preserves_acyclicity,
    bayesian_score_component,
    bayesian_score_components,
    bayesian_score

include("utils.jl")
include("bayes_nets.jl")
include("io.jl")
include("sampling.jl")
include("gibbs.jl")
include("learning.jl")
include("gen_bayes_nets.jl")

include(joinpath("DiscreteBayesNet", "ndgrid.jl"))
include(joinpath("DiscreteBayesNet", "tables.jl"))
include(joinpath("DiscreteBayesNet", "dirichlet_priors.jl"))
include(joinpath("DiscreteBayesNet", "discrete_bayes_net.jl"))
include(joinpath("DiscreteBayesNet", "structure_scoring.jl"))
include(joinpath("DiscreteBayesNet", "greedy_hill_climbing.jl"))
include(joinpath("DiscreteBayesNet", "scan_greedy_hill_climbing.jl"))
include(joinpath("DiscreteBayesNet", "io.jl"))

include(joinpath("Factors", "factors.jl"))
include(joinpath("Inference", "inference.jl"))
include(joinpath("Acausal", "acausal.jl"))

end # module


text:
# Run forward model

println("Time to load BITEModel.jl")
@time using BITEModel
const BM=BITEModel
const I2=BM.ITMIX2

# Parameters to adjust:
#######################

glaciertype = [:itmix, :synth, :peninsula, :rgi][2]
if !isdefined(:plotyes)
    plotyes = true
end

if glaciertype==:itmix
    # run forward model for ITMIX glaciers
    glaciers = [:Academy, :Aqqutikitsoq, :Austfonna, :Brewster, :Columbia,
                :Devon, :Elbrus, :Freya, :Hellstugubreen, :Kesselwandferner, :Mocho,
                :NorthGlacier, :SouthGlacier, :Starbuck, :Synthetic1, :Synthetic2,
                :Synthetic3, :Tasman, :Unteraar, :Urumqi, :Washmawapta]
    glacier = :Unteraar
    pl_kws = Dict{Symbol,Any}(:use_glogem => false)
    gid = ITMIXGlacier(glacier)
elseif glaciertype==:synth
    pl_kws = Dict{Symbol,Any}()
    if haskey(pl_kws, :use_glogem)
        pop!(pl_kws, :use_glogem)
    end
    gid = SyntheticGlacier(:bench)
    pl_kws = Dict{Symbol,Any}()
elseif glaciertype==:peninsula
    flaskid,starbuckid = BM.PeninsulaGlacier(848), BM.PeninsulaGlacier(900)
    gid = [flaskid,starbuckid][2]
    pl_kws = Dict{Symbol,Any}()
else
    error("not implemented")
end
update_cache = false
verbose = true

# Run it
########
gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose;
                                    update_cache=update_cache, pl_kws...)
pm = BM.MPara(pm, iv_nye=true, iv_dist_exp2=0.01)
print("Time to run forward model:"); gc()
@time fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

qtot1d, qd1d = postproc1d(gb, pp, pm)

if plotyes
    print("Time to plot:")
    @time display(BM.plot_run(fwdsol, reuse=false, plot2d=[:h,:iv,:ivlog][1]))
end


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : false
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 256 - 1
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 256 - 1
###
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    # CUDA.device!(6)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics, LinearAlgebra, MAT

@parallel function smooth!(A2::Data.Array, A::Data.Array, fact::Data.Number)
    @inn(A2) = @inn(A) + 1.0/4.1/fact*(@d2_xi(A) + @d2_yi(A))
    return
end

@parallel function compute_maxloc!(Musτ2::Data.Array, Musτ::Data.Array)
    @inn(Musτ2) = @maxloc(Musτ)
    return
end

macro Mu_eff(ix,iy) esc(:( 1.0/(1.0/Musτ[$ix,$iy] + 1.0/(G*dt)) )) end
@parallel_indices (ix,iy) function compute_iter_params!(dτ_Rho::Data.Array, Gdτ::Data.Array, Musτ::Data.Array, Vpdτ::Data.Number, G::Data.Number, dt::Data.Number, Re::Data.Number, r::Data.Number, max_lxy::Data.Number)
    if (ix<=size(dτ_Rho,1)  && iy<=size(dτ_Rho,2))  dτ_Rho[ix,iy] = Vpdτ*max_lxy/Re/@Mu_eff(ix,iy)  end
    if (ix<=size(Gdτ,1)     && iy<=size(Gdτ,2))     Gdτ[ix,iy]    = Vpdτ^2/dτ_Rho[ix,iy]/(r+2)  end
    return
end

@parallel_indices (ix,iy) function assign_τ!(τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τxy_o::Data.Array)
    if (ix<=size(τxx,1) && iy<=size(τxx,2))  τxx_o[ix,iy] = τxx[ix,iy]  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2))  τyy_o[ix,iy] = τyy[ix,iy]  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2))  τxy_o[ix,iy] = τxy[ix,iy]  end
    return
end

macro av_Gdτ(ix,iy) esc(:( 0.25*(Gdτ[$ix,$iy]+Gdτ[$ix+1,$iy]+Gdτ[$ix,$iy+1]+Gdτ[$ix+1,$iy+1]) )) end
macro Gr(ix,iy)     esc(:( Gdτ[$ix,$iy]/(G*dt) )) end
macro av_Gr(ix,iy)  esc(:( @av_Gdτ($ix,$iy)/(G*dt) )) end
macro av_Mus(ix,iy) esc(:( 0.25*(Mus[$ix,$iy]+Mus[$ix+1,$iy]+Mus[$ix,$iy+1]+Mus[$ix+1,$iy+1]) )) end
@parallel_indices (ix,iy) function compute_Pt_τ!(Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τxy_o::Data.Array, Gdτ::Data.Array, Vx::Data.Array, Vy::Data.Array, Mus::Data.Array, r::Data.Number, G::Data.Number, dt::Data.Number, _dx::Data.Number, _dy::Data.Number)
    if (ix<=size(Pt,1)  && iy<=size(Pt,2))    Pt[ix,iy] = Pt[ix,iy] - r*Gdτ[ix,iy]*( _dx*(Vx[ix+1,iy] - Vx[ix,iy]) + _dy*(Vy[ix,iy+1] - Vy[ix,iy]) )  end
    if (ix<=size(τxx,1) && iy<=size(τxx,2))  τxx[ix,iy] = (τxx[ix,iy] + τxx_o[ix,iy]*   @Gr(ix,iy) + 2.0*Gdτ[ix,iy]*(_dx*(Vx[ix+1,iy] - Vx[ix,iy]))) / (1.0 + Gdτ[ix,iy]/Mus[ix,iy] + @Gr(ix,iy))  end
    if (ix<=size(τyy,1) && iy<=size(τyy,2))  τyy[ix,iy] = (τyy[ix,iy] + τyy_o[ix,iy]*   @Gr(ix,iy) + 2.0*Gdτ[ix,iy]*(_dy*(Vy[ix,iy+1] - Vy[ix,iy]))) / (1.0 + Gdτ[ix,iy]/Mus[ix,iy] + @Gr(ix,iy))  end
    if (ix<=size(τxy,1) && iy<=size(τxy,2))  τxy[ix,iy] = (τxy[ix,iy] + τxy_o[ix,iy]*@av_Gr(ix,iy) + 2.0*@av_Gdτ(ix,iy) * 0.5*(_dy*(Vx[ix+1,iy+1] - Vx[ix+1,iy]) + _dx*(Vy[ix+1,iy+1] - Vy[ix,iy+1])) ) / (1.0 + @av_Gdτ(ix,iy)/@av_Mus(ix,iy) + @av_Gr(ix,iy))  end
    return
end

@parallel_indices (ix,iy) function compute_V!(Vx::Data.Array, Vy::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, dτ_Rho::Data.Array, _dx::Data.Number, _dy::Data.Number, size_innVx_1, size_innVx_2, size_innVy_1, size_innVy_2)
    if (ix<=size_innVx_1 && iy<=size_innVx_2)  Vx[ix+1,iy+1] = Vx[ix+1,iy+1] + (_dx*(τxx[ix+1,iy+1] - τxx[ix,iy+1]) + _dy*(τxy[ix,iy+1] - τxy[ix,iy]) - _dx*(Pt[ix+1,iy+1] - Pt[ix,iy+1])) * (0.5(dτ_Rho[ix,iy+1] + dτ_Rho[ix+1,iy+1]))  end
    if (ix<=size_innVy_1 && iy<=size_innVy_2)  Vy[ix+1,iy+1] = Vy[ix+1,iy+1] + (_dy*(τyy[ix+1,iy+1] - τyy[ix+1,iy]) + _dx*(τxy[ix+1,iy] - τxy[ix,iy]) - _dy*(Pt[ix+1,iy+1] - Pt[ix+1,iy])) * (0.5(dτ_Rho[ix+1,iy] + dτ_Rho[ix+1,iy+1]))  end
    return
end

@parallel_indices (ix,iy) function compute_Res!(∇V::Data.Array, Rx::Data.Array, Ry::Data.Array, Pt::Data.Array, Vx::Data.Array, Vy::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, _dx::Data.Number, _dy::Data.Number)
    if (ix<=size(∇V,1) && iy<=size(∇V,2))  ∇V[ix,iy] = _dx*(Vx[ix+1,iy] - Vx[ix,iy]) + _dy*(Vy[ix,iy+1] - Vy[ix,iy])  end
    if (ix<=size(Rx,1) && iy<=size(Rx,2))  Rx[ix,iy] = _dx*(τxx[ix+1,iy+1] - τxx[ix,iy+1]) + _dy*(τxy[ix,iy+1] - τxy[ix,iy]) - _dx*(Pt[ix+1,iy+1] - Pt[ix,iy+1])  end
    if (ix<=size(Ry,1) && iy<=size(Ry,2))  Ry[ix,iy] = _dy*(τyy[ix+1,iy+1] - τyy[ix+1,iy]) + _dx*(τxy[ix+1,iy] - τxy[ix,iy]) - _dy*(Pt[ix+1,iy+1] - Pt[ix+1,iy])  end
    return
end

@parallel_indices (iy) function bc_x!(A::Data.Array)
    A[1  , iy] = A[2    , iy]
    A[end, iy] = A[end-1, iy]
    return
end

@parallel_indices (ix) function bc_y!(A::Data.Array)
    A[ix, 1  ] = A[ix, 2    ]
    A[ix, end] = A[ix, end-1]
    return
end

@views function Stokes2D_()
    # Physics
    lx, ly    = 10.0, 10.0  # domain extends
    ξ         = 1.0         # Maxwell relaxation time
    μs0       = 1.0         # matrix viscosity
    μsi       = 1e-3        # inclusion viscosity
    G         = 1.0         # elastic shear modulus
    εbg       = 1.0         # background strain-rate
    dt        = μs0/(G*ξ)
    # Numerics
    # nx, ny    = 1*128-1, 1*128-1    # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    nt        = 1#5           # number of time steps
    iterMax   = 100#2e5         # maximum number of pseudo-transient iterations
    nout      = 2000         # error checking frequency
    Re        = 5π          # Reynolds number
    r         = 1.0         # Bulk to shear elastic modulus ratio
    CFL       = 0.8/sqrt(2) # CFL number # DEBUG was 0.9
    ε         = 1e-8        # nonlinear absolute tolerence
    # Derived numerics
    dx, dy    = lx/nx, ly/ny # cell sizes
    max_lxy   = max(lx,ly)
    Vpdτ      = min(dx,dy)*CFL
    _dx, _dy  = 1.0/dx, 1.0/dy
    xc, yc, yv = LinRange(dx/2, lx - dx/2, nx), LinRange(dy/2, ly - dy/2, ny), LinRange(0, ly, ny+1)
    # Array allocations
    Pt        = @zeros(nx  ,ny  )
    dτ_Rho    = @zeros(nx  ,ny  )
    Gdτ       = @zeros(nx  ,ny  )
    ∇V        = @zeros(nx  ,ny  )
    τxx       = @zeros(nx  ,ny  )
    τyy       = @zeros(nx  ,ny  )
    τxy       = @zeros(nx-1,ny-1)
    τxx_o     = @zeros(nx  ,ny  )
    τyy_o     = @zeros(nx  ,ny  )
    τxy_o     = @zeros(nx-1,ny-1)
    Rx        = @zeros(nx-1,ny-2)
    Ry        = @zeros(nx-2,ny-1)
    Mus2      = @zeros(nx  ,ny  )
    Musτ      = @zeros(nx  ,ny  )
    # Initial conditions
    Rad2      =  zeros(nx  ,ny  )
    Vx        =  zeros(nx+1,ny  )
    Vy        =  zeros(nx  ,ny+1)
    Rad2     .= [((ix-1)*dx +0.5*dx -0.5*lx)^2 + ((iy-1)*dy +0.5*dy -0.5*ly)^2 for ix=1:size(Rad2,1), iy=1:size(Rad2,2)]
    Vx        = Data.Array( -εbg.*[((ix-1)*dx -0.5*lx) for ix=1:size(Vx,1), iy=1:size(Vx,2)] )
    Vy        = Data.Array(  εbg.*[((iy-1)*dy -0.5*ly) for ix=1:size(Vy,1), iy=1:size(Vy,2)] )
    Mus       = μs0*ones(nx,ny)
    Mus[Rad2.<1.0] .= μsi
    Mus       = Data.Array( Mus )
    Mus2     .= Mus
    for ism=1:10
        @parallel smooth!(Mus2, Mus, 1.0)
        Mus, Mus2 = Mus2, Mus
    end
    Musτ     .= Mus
    @parallel compute_maxloc!(Musτ, Mus)
    @parallel (1:size(Musτ,2)) bc_x!(Musτ)
    @parallel (1:size(Musτ,1)) bc_y!(Musτ)
    size_innVx_1, size_innVx_2 = size(Vx,1)-2, size(Vx,2)-2
    size_innVy_1, size_innVy_2 = size(Vy,1)-2, size(Vy,2)-2
    # Time loop
    @parallel compute_iter_params!(dτ_Rho, Gdτ, Musτ, Vpdτ, G, dt, Re, r, max_lxy)
    t=0.0; ittot=0; evo_t=[]; evo_τyy=[]; err_evo1=[]; err_evo2=[]; t_tic = 0.0
    for it = 1:nt
        err=2*ε; iter=0
        @parallel assign_τ!(τxx, τyy, τxy, τxx_o, τyy_o, τxy_o)
        # Pseudo-transient iteration
        while err > ε && iter <= iterMax
            if (it==1 && iter==11) t_tic = Base.time()  end
            @parallel compute_Pt_τ!(Pt, τxx, τyy, τxy, τxx_o, τyy_o, τxy_o, Gdτ, Vx, Vy, Mus, r, G, dt, _dx, _dy)
            @parallel compute_V!(Vx, Vy, Pt, τxx, τyy, τxy, dτ_Rho, _dx, _dy, size_innVx_1, size_innVx_2, size_innVy_1, size_innVy_2)
            @parallel (1:size(Vx,1)) bc_y!(Vx)
            @parallel (1:size(Vy,2)) bc_x!(Vy)
            iter += 1
            if iter % nout == 0
                @parallel compute_Res!(∇V, Rx, Ry, Pt, Vx, Vy, τxx, τyy, τxy, _dx, _dy)
                Vmin, Vmax = minimum(Vx), maximum(Vx)
                Pmin, Pmax = minimum(Pt), maximum(Pt)
                norm_Rx    = norm(Rx)/(Pmax-Pmin)*lx/sqrt(length(Rx))
                norm_Ry    = norm(Ry)/(Pmax-Pmin)*lx/sqrt(length(Ry))
                norm_∇V    = norm(∇V)/(Vmax-Vmin)*lx/sqrt(length(∇V))
                # norm_Rx = norm(Rx)/length(Rx); norm_Ry = norm(Ry)/length(Ry); norm_∇V = norm(∇V)/length(∇V)
                err = maximum([norm_Rx, norm_Ry, norm_∇V])
                if isnan(err) error("NaN") end
                push!(err_evo1, maximum([norm_Rx, norm_Ry, norm_∇V])); push!(err_evo2,iter)
                @printf("Step = %d, iter = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_∇V=%1.3e] \n", it, iter, err, norm_Rx, norm_Ry, norm_∇V)
            end
        end
        ittot += iter; t += dt
        # push!(evo_t, t); push!(evo_τyy, maximum(τyy))
    end
    # Performance
    t_toc    = Base.time() - t_tic
    A_eff    = (6*2 + 3*1 + 3)/1e9*nx*ny*sizeof(Data.Number) # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    t_it     = t_toc/(ittot-10)                        # Execution time per iteration [s]
    T_eff    = A_eff/t_it                              # Effective memory throughput [GB/s]
    @printf("Total iters = %d (%d steps), time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", ittot, nt, t_toc, round(T_eff, sigdigits=3))
    # Visualisation
    if do_viz
        p1 = heatmap(xc, yc, Array(τyy)', aspect_ratio=1, xlims=extrema(xc), ylims=extrema(yc), c=:viridis, title="τyy")
        p2 = heatmap(xc, yv, Array(Vy)', aspect_ratio=1, xlims=extrema(xc), ylims=extrema(yv), c=:viridis, title="Vy")
        # p4 = heatmap(X[2:end-1], Yv[2:end-1], log10.(abs.(Array(Ry)')), aspect_ratio=1, xlims=(X[2],X[end-1]), ylims=(Yv[2],Yv[end-1]), c=:viridis, title="log10(Ry)")
        # p5 = scatter(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, framestyle=:box, labels="max(error)", yaxis=:log10)
        p3 = plot(evo_t, evo_τyy , legend=false, xlabel="time", ylabel="max(τyy)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
            #plot!(evo_t, 2.0.*εbg.*μs0.*(1.0.-exp.(.-evo_t.*G./μs0)), linewidth=2.0) # analytical solution
        display(plot(p1, p2, p3))
    end
    if do_save
        !ispath("../output") && mkdir("../output")
        open("../output/out_Stokes2D_ve_perf.txt","a") do io
            println(io, "$(nx) $(ny) $(ittot) $(t_toc) $(A_eff) $(t_it) $(T_eff)")
        end
    end
    if do_save_viz
        !ispath("../out_visu") && mkdir("../out_visu")
        matwrite("../out_visu/Stokes_2D_ve_perf.mat", Dict("Pt_2D"=> Array(Pt), "Mus_2D"=> Array(Mus), "Txy_2D"=> Array(τxy), "Vy_2D"=> Array(Vy), "dx_2D"=> dx, "dy_2D"=> dy); compress = true)
    end
    return xc, yc, Pt
end

if use_return
    xc, yc, P = Stokes2D_();
else
    Stokes2D = begin Stokes2D_(); return; end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    maximalcliques, maximumclique,
    maximalconncliques, maximumconnclique

mutable struct FindCliqueState{T<:UndirectedGraph}
    graph::T
    targetsize::Union{Int,Nothing} # Int

    adjacencies::Dict{Int,Set{Int}}
    expire::Union{UInt64,Nothing} # UInt64, nanoseconds
    Q::Vector{Int}

    cliques::Vector{Set{Int}}
    status::Symbol

    function FindCliqueState{T}(graph; timeout=nothing, targetsize=nothing,
                kwargs...) where {T<:UndirectedGraph}
        if timeout !== nothing
            expire = (time_ns() + timeout * 1_000_000_000)::UInt64
        else
            expire = nothing
        end
        # fast adjacency access
        adj = Dict{Int,Set{Int}}()
        for n in nodeset(graph)
            adj[n] = adjacencies(graph, n)
        end
        new(graph, targetsize, adj, expire, [], [], :ongoing)
    end
end

mutable struct FindConnCliqueState{T<:ModularProduct}
    graph::T
    targetsize::Union{Int,Nothing} # Int

    adjacencies::Dict{Int,Set{Int}}
    connected::Dict{Int,Set{Int}}
    disconn::Dict{Int,Set{Int}}
    expire::Union{UInt64,Nothing} # UInt64, nanoseconds

    cliques::Vector{Set{Int}}
    status::Symbol

    function FindConnCliqueState{T}(graph; timeout=nothing, targetsize=nothing,
                kwargs...) where {T<:ModularProduct}
        if timeout !== nothing
            expire = (time_ns() + timeout * 1_000_000_000)::UInt64
        else
            expire = nothing
        end
        # fast adjacency access
        adj = Dict{Int,Set{Int}}()
        conn = Dict{Int,Set{Int}}()
        disconn = Dict{Int,Set{Int}}()
        for n in nodeset(graph)
            adj[n] = Set{Int}()
            conn[n] = Set{Int}()
            disconn[n] = Set{Int}()
            for (i, a) in neighbors(graph, n)
                push!(adj[n], a)
                if edgeattr(graph, i).hasedge
                    push!(conn[n], a)
                else
                    push!(disconn[n], a)
                end
            end
        end
        new(graph, targetsize, adj, conn, disconn, expire, [], :ongoing)
    end
end

function expand!(state::FindCliqueState, subg, cand)
    (state.status == :timedout || state.status == :targetreached) && return
    if isempty(subg)
        # Report max clique
        push!(state.cliques, Set(state.Q))
        return
    elseif state.expire !== nothing && time_ns() > state.expire
        state.status = :timedout
        return
    elseif state.targetsize !== nothing && length(state.Q) >= state.targetsize
        state.status = :targetreached
        push!(state.cliques, Set(state.Q))
        return
    end
    candnbrcnt(n) = length(intersect(cand, state.adjacencies[n]))
    pivot = sortstablemax(subg, by=candnbrcnt)
    copv = setdiff(cand, state.adjacencies[pivot])

    for q in copv
        push!(state.Q, q)
        qnbrs = state.adjacencies[q]
        subgq = intersect(subg, qnbrs)
        candq = intersect(cand, qnbrs)
        expand!(state, subgq, candq)
        pop!(cand, q)
        pop!(state.Q)
    end
    return
end

function expandconn!(state::FindConnCliqueState, R, P, Q, X, Y)
    (state.status == :timedout || state.status == :targetreached) && return
    if isempty(P) && isempty(X)
        # Report max clique
        push!(state.cliques, copy(R))
        return
    elseif state.expire !== nothing && time_ns() > state.expire
        state.status = :timedout
        return
    elseif state.targetsize !== nothing && length(R) >= state.targetsize
        state.status = :targetreached
        push!(state.cliques, copy(R))
        return
    end
    while !isempty(P)
        n = pop!(P)
        Rnew = union(R, [n])
        Qnew = intersect(Q, state.disconn[n])
        Pnew = union(
            intersect(P, state.adjacencies[n]),
            intersect(Q, state.connected[n]))
        Ynew = intersect(Y, state.disconn[n])
        Xnew = union(
            intersect(X, state.adjacencies[n]),
            intersect(Y, state.connected[n]))
        expandconn!(state, Rnew, Pnew, Qnew, Xnew, Ynew)
        push!(X, n)
    end
    return
end

"""
    maximalcliques(graph::UndirectedGraph; kwargs...
        ) -> Tuple{Vector{Set{Int}}, Symbol}

Return maximal cliques.

# Reference

1. Tomita, E., Tanaka, A., & Takahashi, H. (2006). The worst-case time
   complexity for generating all maximal cliques and computational experiments.
   Theoretical Computer Science, 363(1), 28–42.
   https://doi.org/10.1016/J.TCS.2006.06.015
1. Cazals, F., & Karande, C. (2008). A note on the problem of reporting maximal
   cliques. Theoretical Computer Science, 407(1–3), 564–568.
   https://doi.org/10.1016/j.tcs.2008.05.010
"""
function maximalcliques(graph::T; kwargs...) where {T<:UndirectedGraph}
    state = FindCliqueState{T}(graph; kwargs...)
    expand!(state, nodeset(graph), nodeset(graph))
    if state.status == :ongoing
        state.status = :done
    end
    return (state.cliques, state.status)
end

"""
    maximumclique(graph::UndirectedGraph; kwargs...) -> Tuple{Set{Int}, Symbol}

Return a maximum clique.

"""
function maximumclique(graph::T; kwargs...) where {T<:UndirectedGraph}
    (cliques, status) = maximalcliques(graph; kwargs...)
    return (sortstablemax(cliques, by=length, init=[]), status)
end

"""
    maximalconncliques(graph::ModularProduct; kwargs...
        ) -> Tuple{Vector{Set{Int}}, Symbol}

Return maximal connected cliques.

# Reference

1. Cazals, F., & Karande, C. (2005). An algorithm for reporting maximal
   c-cliques. Theoretical Computer Science, 349(3), 484–490.
   https://doi.org/10.1016/j.tcs.2005.09.038

"""
function maximalconncliques(graph::T; kwargs...) where {T<:UndirectedGraph}
    state = FindConnCliqueState{T}(graph; kwargs...)
    nodes = nodeset(graph)
    done = Set{Int}()
    for n in nodes
        R = Set([n])
        P = intersect(setdiff(nodes, done), state.connected[n])
        Q = intersect(setdiff(nodes, done), state.disconn[n])
        X = intersect(state.connected[n], done)
        Y = intersect(state.disconn[n], done)
        expandconn!(state, R, P, Q, X, Y)
        push!(done, n)
    end
    if state.status == :ongoing
        state.status = :done
    end
    return (state.cliques, state.status)
end

"""
    maximumconnclique(graph::ModularProduct; kwargs...
        ) -> Tuple{Set{Int}, Symbol}

Return a maximum connected clique.

"""
function maximumconnclique(graph::T; kwargs...) where {T<:ModularProduct}
    (cliques, status) = maximalconncliques(graph; kwargs...)
    return (sortstablemax(cliques, by=length, init=[]), status)
end


text:
struct ContextSpecifier
    varName::Symbol
    ctxType::Any
    subtype::Bool
end

macro contextualized(expr::Expr)
    if Meta.isexpr(expr, :function)
        return contextualizedFunction(expr, false)
    elseif Meta.isexpr(expr, :(=))
        return contextualizedFunction(expr, true)
    else
        @error "expected a function definition"
    end
end

function contextualizedFunction(expr::Expr, isshort)
    # Parse complete function
    f = splitdef(expr)
    b = f[:body]

    isinline = (b.args[1] == Expr(:meta,:inline))
    if isshort
        b = b.args[isinline ? 2 : 1 ]
    end

    # Extract context specification
    ctxIdx = isshort ? 1 : 2
    if isinline && !isshort
        ctxIdx += 1
    end
    contextSpec = parseContextSpecifier(b.args[ctxIdx])

    # Remove context specification, i.e. @with macro
    deleteat!(b.args, ctxIdx)

    # Push context and function
    if contextSpec.varName == :(_)
        insert!(f[:args], 1, :(::$(contextSpec.ctxType)))
    else
        insert!(f[:args], 1, :($(contextSpec.varName)::$(contextSpec.ctxType)))
    end
    insert!(f[:args], 2, :(::typeof($(f[:name]))))

    # Change function name
    f[:name] = :($TinyCassette.execute)

    esc(MacroTools.combinedef(f))
end

function parseContextSpecifier(macrocall)
    @assert isexpr(macrocall, :macrocall) "expected a @with context specification"
    @assert macrocall.args[1] == Symbol("@with") "expected a @with specification"

    braces = macrocall.args[3]
    @assert isexpr(braces, :braces) "context specification should be in braces"

    spec = braces.args[1]
    if isexpr(spec, :(::))
        subtype = false
    elseif isexpr(spec, :(<:))
        subtype = true
    else
        @error "context specification should start with ctx::Ctx or ctx<:Ctx"
    end

    varName = spec.args[1]
    ctxType = spec.args[2]
    return ContextSpecifier(varName, ctxType, subtype)
end


text:
# Unique Work Queue

struct WorkQueue{T}
  items::Vector{T}
end

WorkQueue{T}() where T = WorkQueue(T[])

WorkQueue() = WorkQueue{Any}()

WorkQueue(xs) = WorkQueue(collect(xs))

function Base.push!(q::WorkQueue, x)
  i = findfirst(==(x), q.items)
  i === nothing || (deleteat!(q.items, i))
  push!(q.items, x)
  return q
end

Base.pop!(q::WorkQueue) = pop!(q.items)
Base.isempty(q::WorkQueue) = isempty(q.items)

# Simple Graph

struct Graph{T}
  nodes::Dict{T,Int}
  links::Vector{Set{T}}
end

Graph{T}() where T = Graph{T}(Dict(), [])

Base.length(g::Graph) = length(g.links)
Base.keys(g::Graph) = keys(g.nodes)

function Base.getindex(g::Graph, x)
  if !haskey(g.nodes, x)
    push!(g.links, Set())
    g.nodes[x] = length(g)
  end
  return g.links[g.nodes[x]]
end

function connect!(g::Graph, a, b)
  push!(g[a], b)
  push!(g[b], a)
  return g
end