text:
const unix = 1
const linux = 1
const WCSERR_MSG_LENGTH = 160
# Skipping MacroDefinition: ERRLEN(sizeof(structwcserr)/sizeof(int))
# Skipping MacroDefinition: WCSERR_SET(status)err,status,function,__FILE__,__LINE__
# Skipping MacroDefinition: LINLEN(sizeof(structlinprm)/sizeof(int))
const PVN = 30
const PRJX2S_ARGS =
const PRJS2X_ARGS =
# Skipping MacroDefinition: PRJLEN(sizeof(structprjprm)/sizeof(int))
# Skipping MacroDefinition: CELLEN(sizeof(structcelprm)/sizeof(int))
# Skipping MacroDefinition: SPXLEN(sizeof(structspxprm)/sizeof(int))
const SPX_ARGS =
# Skipping MacroDefinition: SPCLEN(sizeof(structspcprm)/sizeof(int))
# Skipping MacroDefinition: TABLEN(sizeof(structtabprm)/sizeof(int))
const WCSSUB_LONGITUDE = 0x1001
const WCSSUB_LATITUDE = 0x1002
const WCSSUB_CUBEFACE = 0x1004
const WCSSUB_CELESTIAL = 0x1007
const WCSSUB_SPECTRAL = 0x1008
const WCSSUB_STOKES = 0x1010
# Skipping MacroDefinition: WCSLEN(sizeof(structwcsprm)/sizeof(int))
# Skipping MacroDefinition: wcscopy(alloc,wcssrc,wcsdst)wcssub(alloc,wcssrc,0x0,0x0,wcsdst)
const WCSHDR_none = 0x00000000
const WCSHDR_all = 0x000FFFFF
const WCSHDR_reject = 0x10000000
const WCSHDR_CROTAia = 0x00000001
const WCSHDR_EPOCHa = 0x00000002
const WCSHDR_VELREFa = 0x00000004
const WCSHDR_CD00i00j = 0x00000008
const WCSHDR_PC00i00j = 0x00000010
const WCSHDR_PROJPn = 0x00000020
const WCSHDR_RADECSYS = 0x00000040
const WCSHDR_VSOURCE = 0x00000080
const WCSHDR_DOBSn = 0x00000100
const WCSHDR_LONGKEY = 0x00000200
const WCSHDR_CNAMn = 0x00000400
const WCSHDR_AUXIMG = 0x00000800
const WCSHDR_ALLIMG = 0x00001000
const WCSHDR_IMGHEAD = 0x00010000
const WCSHDR_BIMGARR = 0x00020000
const WCSHDR_PIXLIST = 0x00040000
const WCSHDO_none = 0x00
const WCSHDO_all = 0xFF
const WCSHDO_safe = 0x0F
const WCSHDO_DOBSn = 0x01
const WCSHDO_TPCn_ka = 0x02
const WCSHDO_PVn_ma = 0x04
const WCSHDO_CRPXna = 0x08
const WCSHDO_CNAMna = 0x10
const WCSHDO_WCSNna = 0x20
immutable Array_160_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
    d5::Uint8
    d6::Uint8
    d7::Uint8
    d8::Uint8
    d9::Uint8
    d10::Uint8
    d11::Uint8
    d12::Uint8
    d13::Uint8
    d14::Uint8
    d15::Uint8
    d16::Uint8
    d17::Uint8
    d18::Uint8
    d19::Uint8
    d20::Uint8
    d21::Uint8
    d22::Uint8
    d23::Uint8
    d24::Uint8
    d25::Uint8
    d26::Uint8
    d27::Uint8
    d28::Uint8
    d29::Uint8
    d30::Uint8
    d31::Uint8
    d32::Uint8
    d33::Uint8
    d34::Uint8
    d35::Uint8
    d36::Uint8
    d37::Uint8
    d38::Uint8
    d39::Uint8
    d40::Uint8
    d41::Uint8
    d42::Uint8
    d43::Uint8
    d44::Uint8
    d45::Uint8
    d46::Uint8
    d47::Uint8
    d48::Uint8
    d49::Uint8
    d50::Uint8
    d51::Uint8
    d52::Uint8
    d53::Uint8
    d54::Uint8
    d55::Uint8
    d56::Uint8
    d57::Uint8
    d58::Uint8
    d59::Uint8
    d60::Uint8
    d61::Uint8
    d62::Uint8
    d63::Uint8
    d64::Uint8
    d65::Uint8
    d66::Uint8
    d67::Uint8
    d68::Uint8
    d69::Uint8
    d70::Uint8
    d71::Uint8
    d72::Uint8
    d73::Uint8
    d74::Uint8
    d75::Uint8
    d76::Uint8
    d77::Uint8
    d78::Uint8
    d79::Uint8
    d80::Uint8
    d81::Uint8
    d82::Uint8
    d83::Uint8
    d84::Uint8
    d85::Uint8
    d86::Uint8
    d87::Uint8
    d88::Uint8
    d89::Uint8
    d90::Uint8
    d91::Uint8
    d92::Uint8
    d93::Uint8
    d94::Uint8
    d95::Uint8
    d96::Uint8
    d97::Uint8
    d98::Uint8
    d99::Uint8
    d100::Uint8
    d101::Uint8
    d102::Uint8
    d103::Uint8
    d104::Uint8
    d105::Uint8
    d106::Uint8
    d107::Uint8
    d108::Uint8
    d109::Uint8
    d110::Uint8
    d111::Uint8
    d112::Uint8
    d113::Uint8
    d114::Uint8
    d115::Uint8
    d116::Uint8
    d117::Uint8
    d118::Uint8
    d119::Uint8
    d120::Uint8
    d121::Uint8
    d122::Uint8
    d123::Uint8
    d124::Uint8
    d125::Uint8
    d126::Uint8
    d127::Uint8
    d128::Uint8
    d129::Uint8
    d130::Uint8
    d131::Uint8
    d132::Uint8
    d133::Uint8
    d134::Uint8
    d135::Uint8
    d136::Uint8
    d137::Uint8
    d138::Uint8
    d139::Uint8
    d140::Uint8
    d141::Uint8
    d142::Uint8
    d143::Uint8
    d144::Uint8
    d145::Uint8
    d146::Uint8
    d147::Uint8
    d148::Uint8
    d149::Uint8
    d150::Uint8
    d151::Uint8
    d152::Uint8
    d153::Uint8
    d154::Uint8
    d155::Uint8
    d156::Uint8
    d157::Uint8
    d158::Uint8
    d159::Uint8
    d160::Uint8
end
immutable wcserr
    status::Cint
    line_no::Cint
    _function::Ptr{Uint8}
    file::Ptr{Uint8}
    msg::Array_160_Uint8
end
# begin enum lin_errmsg_enum
typealias lin_errmsg_enum Uint32
const LINERR_SUCCESS = 0
const LINERR_NULL_POINTER = 1
const LINERR_MEMORY = 2
const LINERR_SINGULAR_MTX = 3
# end enum lin_errmsg_enum
immutable linprm
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    piximg::Ptr{Cdouble}
    imgpix::Ptr{Cdouble}
    unity::Cint
    padding::Cint
    err::Ptr{wcserr}
    i_naxis::Cint
    m_flag::Cint
    m_naxis::Cint
    m_padding::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    padding2::Ptr{Void}
end
# begin enum prj_errmsg_enum
typealias prj_errmsg_enum Uint32
const PRJERR_SUCCESS = 0
const PRJERR_NULL_POINTER = 1
const PRJERR_BAD_PARAM = 2
const PRJERR_BAD_PIX = 3
const PRJERR_BAD_WORLD = 4
# end enum prj_errmsg_enum
immutable Array_4_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
end
immutable Array_30_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
    d7::Cdouble
    d8::Cdouble
    d9::Cdouble
    d10::Cdouble
    d11::Cdouble
    d12::Cdouble
    d13::Cdouble
    d14::Cdouble
    d15::Cdouble
    d16::Cdouble
    d17::Cdouble
    d18::Cdouble
    d19::Cdouble
    d20::Cdouble
    d21::Cdouble
    d22::Cdouble
    d23::Cdouble
    d24::Cdouble
    d25::Cdouble
    d26::Cdouble
    d27::Cdouble
    d28::Cdouble
    d29::Cdouble
    d30::Cdouble
end
immutable Array_40_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
    d5::Uint8
    d6::Uint8
    d7::Uint8
    d8::Uint8
    d9::Uint8
    d10::Uint8
    d11::Uint8
    d12::Uint8
    d13::Uint8
    d14::Uint8
    d15::Uint8
    d16::Uint8
    d17::Uint8
    d18::Uint8
    d19::Uint8
    d20::Uint8
    d21::Uint8
    d22::Uint8
    d23::Uint8
    d24::Uint8
    d25::Uint8
    d26::Uint8
    d27::Uint8
    d28::Uint8
    d29::Uint8
    d30::Uint8
    d31::Uint8
    d32::Uint8
    d33::Uint8
    d34::Uint8
    d35::Uint8
    d36::Uint8
    d37::Uint8
    d38::Uint8
    d39::Uint8
    d40::Uint8
end
immutable Array_10_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
    d7::Cdouble
    d8::Cdouble
    d9::Cdouble
    d10::Cdouble
end
immutable prjprm
    flag::Cint
    code::Array_4_Uint8
    r0::Cdouble
    pv::Array_30_Cdouble
    phi0::Cdouble
    theta0::Cdouble
    bounds::Cint
    name::Array_40_Uint8
    category::Cint
    pvrange::Cint
    simplezen::Cint
    equiareal::Cint
    conformal::Cint
    _global::Cint
    divergent::Cint
    x0::Cdouble
    y0::Cdouble
    err::Ptr{wcserr}
    padding::Ptr{Void}
    w::Array_10_Cdouble
    m::Cint
    n::Cint
    prjx2s::Ptr{Void}
    prjs2x::Ptr{Void}
end
# begin enum cel_errmsg_enum
typealias cel_errmsg_enum Uint32
const CELERR_SUCCESS = 0
const CELERR_NULL_POINTER = 1
const CELERR_BAD_PARAM = 2
const CELERR_BAD_COORD_TRANS = 3
const CELERR_ILL_COORD_TRANS = 4
const CELERR_BAD_PIX = 5
const CELERR_BAD_WORLD = 6
# end enum cel_errmsg_enum
immutable Array_4_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
end
immutable Array_5_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
end
immutable celprm
    flag::Cint
    offset::Cint
    phi0::Cdouble
    theta0::Cdouble
    ref::Array_4_Cdouble
    prj::prjprm
    euler::Array_5_Cdouble
    latpreq::Cint
    isolat::Cint
    err::Ptr{wcserr}
    padding::Ptr{Void}
end
immutable spxprm
    restfrq::Cdouble
    restwav::Cdouble
    wavetype::Cint
    velotype::Cint
    freq::Cdouble
    afrq::Cdouble
    ener::Cdouble
    wavn::Cdouble
    vrad::Cdouble
    wave::Cdouble
    vopt::Cdouble
    zopt::Cdouble
    awav::Cdouble
    velo::Cdouble
    beta::Cdouble
    dfreqafrq::Cdouble
    dafrqfreq::Cdouble
    dfreqener::Cdouble
    denerfreq::Cdouble
    dfreqwavn::Cdouble
    dwavnfreq::Cdouble
    dfreqvrad::Cdouble
    dvradfreq::Cdouble
    dfreqwave::Cdouble
    dwavefreq::Cdouble
    dfreqawav::Cdouble
    dawavfreq::Cdouble
    dfreqvelo::Cdouble
    dvelofreq::Cdouble
    dwavevopt::Cdouble
    dvoptwave::Cdouble
    dwavezopt::Cdouble
    dzoptwave::Cdouble
    dwaveawav::Cdouble
    dawavwave::Cdouble
    dwavevelo::Cdouble
    dvelowave::Cdouble
    dawavvelo::Cdouble
    dveloawav::Cdouble
    dvelobeta::Cdouble
    dbetavelo::Cdouble
    err::Ptr{wcserr}
    padding::Ptr{Void}
end
# begin enum spc_errmsg_enum
typealias spc_errmsg_enum Cint
const SPCERR_NO_CHANGE = -1
const SPCERR_SUCCESS = 0
const SPCERR_NULL_POINTER = 1
const SPCERR_BAD_SPEC_PARAMS = 2
const SPCERR_BAD_X = 3
const SPCERR_BAD_SPEC = 4
# end enum spc_errmsg_enum
immutable Array_8_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
    d5::Uint8
    d6::Uint8
    d7::Uint8
    d8::Uint8
end
immutable Array_7_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
    d7::Cdouble
end
immutable Array_6_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
end
immutable spcprm
    flag::Cint
    _type::Array_8_Uint8
    code::Array_4_Uint8
    crval::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    pv::Array_7_Cdouble
    w::Array_6_Cdouble
    isGrism::Cint
    padding1::Cint
    err::Ptr{wcserr}
    padding2::Ptr{Void}
    spxX2P::Ptr{Void}
    spxP2S::Ptr{Void}
    spxS2P::Ptr{Void}
    spxP2X::Ptr{Void}
end
# begin enum tab_errmsg_enum
typealias tab_errmsg_enum Uint32
const TABERR_SUCCESS = 0
const TABERR_NULL_POINTER = 1
const TABERR_MEMORY = 2
const TABERR_BAD_PARAMS = 3
const TABERR_BAD_X = 4
const TABERR_BAD_WORLD = 5
# end enum tab_errmsg_enum
immutable tabprm
    flag::Cint
    M::Cint
    K::Ptr{Cint}
    map::Ptr{Cint}
    crval::Ptr{Cdouble}
    index::Ptr{Ptr{Cdouble}}
    coord::Ptr{Cdouble}
    nc::Cint
    padding::Cint
    sense::Ptr{Cint}
    p0::Ptr{Cint}
    delta::Ptr{Cdouble}
    extrema::Ptr{Cdouble}
    err::Ptr{wcserr}
    m_flag::Cint
    m_M::Cint
    m_N::Cint
    set_M::Cint
    m_K::Ptr{Cint}
    m_map::Ptr{Cint}
    m_crval::Ptr{Cdouble}
    m_index::Ptr{Ptr{Cdouble}}
    m_indxs::Ptr{Ptr{Cdouble}}
    m_coord::Ptr{Cdouble}
end
# begin enum wcs_errmsg_enum
typealias wcs_errmsg_enum Uint32
const WCSERR_SUCCESS = 0
const WCSERR_NULL_POINTER = 1
const WCSERR_MEMORY = 2
const WCSERR_SINGULAR_MTX = 3
const WCSERR_BAD_CTYPE = 4
const WCSERR_BAD_PARAM = 5
const WCSERR_BAD_COORD_TRANS = 6
const WCSERR_ILL_COORD_TRANS = 7
const WCSERR_BAD_PIX = 8
const WCSERR_BAD_WORLD = 9
const WCSERR_BAD_WORLD_COORD = 10
const WCSERR_NO_SOLUTION = 11
const WCSERR_BAD_SUBIMAGE = 12
const WCSERR_NON_SEPARABLE = 13
# end enum wcs_errmsg_enum
immutable pvcard
    i::Cint
    m::Cint
    value::Cdouble
end
immutable Array_72_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
    d5::Uint8
    d6::Uint8
    d7::Uint8
    d8::Uint8
    d9::Uint8
    d10::Uint8
    d11::Uint8
    d12::Uint8
    d13::Uint8
    d14::Uint8
    d15::Uint8
    d16::Uint8
    d17::Uint8
    d18::Uint8
    d19::Uint8
    d20::Uint8
    d21::Uint8
    d22::Uint8
    d23::Uint8
    d24::Uint8
    d25::Uint8
    d26::Uint8
    d27::Uint8
    d28::Uint8
    d29::Uint8
    d30::Uint8
    d31::Uint8
    d32::Uint8
    d33::Uint8
    d34::Uint8
    d35::Uint8
    d36::Uint8
    d37::Uint8
    d38::Uint8
    d39::Uint8
    d40::Uint8
    d41::Uint8
    d42::Uint8
    d43::Uint8
    d44::Uint8
    d45::Uint8
    d46::Uint8
    d47::Uint8
    d48::Uint8
    d49::Uint8
    d50::Uint8
    d51::Uint8
    d52::Uint8
    d53::Uint8
    d54::Uint8
    d55::Uint8
    d56::Uint8
    d57::Uint8
    d58::Uint8
    d59::Uint8
    d60::Uint8
    d61::Uint8
    d62::Uint8
    d63::Uint8
    d64::Uint8
    d65::Uint8
    d66::Uint8
    d67::Uint8
    d68::Uint8
    d69::Uint8
    d70::Uint8
    d71::Uint8
    d72::Uint8
end
immutable pscard
    i::Cint
    m::Cint
    value::Array_72_Uint8
end
immutable wtbarr
    i::Cint
    m::Cint
    kind::Cint
    extnam::Array_72_Uint8
    extver::Cint
    extlev::Cint
    ttype::Array_72_Uint8
    row::Clong
    ndim::Cint
    dimlen::Ptr{Cint}
    arrayp::Ptr{Ptr{Cdouble}}
end
immutable Array_3_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
end
type wcsprm
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    crval::Ptr{Cdouble}
    cunit::Ptr{Void}
    ctype::Ptr{Void}
    lonpole::Cdouble
    latpole::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    npv::Cint
    npvmax::Cint
    pv::Ptr{pvcard}
    nps::Cint
    npsmax::Cint
    ps::Ptr{pscard}
    cd::Ptr{Cdouble}
    crota::Ptr{Cdouble}
    altlin::Cint
    velref::Cint
    alt::Array_4_Uint8
    colnum::Cint
    colax::Ptr{Cint}
    cname::Ptr{Void}
    crder::Ptr{Cdouble}
    csyer::Ptr{Cdouble}
    dateavg::Array_72_Uint8
    dateobs::Array_72_Uint8
    equinox::Cdouble
    mjdavg::Cdouble
    mjdobs::Cdouble
    obsgeo::Array_3_Cdouble
    radesys::Array_72_Uint8
    specsys::Array_72_Uint8
    ssysobs::Array_72_Uint8
    velosys::Cdouble
    zsource::Cdouble
    ssyssrc::Array_72_Uint8
    velangl::Cdouble
    wcsname::Array_72_Uint8
    ntab::Cint
    nwtb::Cint
    tab::Ptr{tabprm}
    wtb::Ptr{wtbarr}
    lngtyp::Array_8_Uint8
    lattyp::Array_8_Uint8
    lng::Cint
    lat::Cint
    spec::Cint
    cubeface::Cint
    types::Ptr{Cint}
    padding::Ptr{Void}
    lin::linprm
    cel::celprm
    spc::spcprm
    err::Ptr{wcserr}
    m_padding::Ptr{Void}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_crval::Ptr{Cdouble}
    m_cunit::Ptr{Void}
    m_ctype::Ptr{Void}
    m_pv::Ptr{pvcard}
    m_ps::Ptr{pscard}
    m_cd::Ptr{Cdouble}
    m_crota::Ptr{Cdouble}
    m_colax::Ptr{Cint}
    m_cname::Ptr{Void}
    m_crder::Ptr{Cdouble}
    m_csyer::Ptr{Cdouble}
    m_tab::Ptr{tabprm}
    m_wtb::Ptr{wtbarr}
    wcsprm() = new(-1)
end
# begin enum wcshdr_errmsg_enum
typealias wcshdr_errmsg_enum Uint32
const WCSHDRERR_SUCCESS = 0
const WCSHDRERR_NULL_POINTER = 1
const WCSHDRERR_MEMORY = 2
const WCSHDRERR_BAD_COLUMN = 3
const WCSHDRERR_PARSER = 4
const WCSHDRERR_BAD_TABULAR_PARAMS = 5
# end enum wcshdr_errmsg_enum
immutable Array_28_Int16
    d1::Int16
    d2::Int16
    d3::Int16
    d4::Int16
    d5::Int16
    d6::Int16
    d7::Int16
    d8::Int16
    d9::Int16
    d10::Int16
    d11::Int16
    d12::Int16
    d13::Int16
    d14::Int16
    d15::Int16
    d16::Int16
    d17::Int16
    d18::Int16
    d19::Int16
    d20::Int16
    d21::Int16
    d22::Int16
    d23::Int16
    d24::Int16
    d25::Int16
    d26::Int16
    d27::Int16
    d28::Int16
end


text:
#
# Date created: 2022-03-14
# Author: aradclif
#
#
############################################################################################
function sizeblock(N::Int)
    block = Expr(:block)
    for d = 1:N
        ex = Expr(:(=), Symbol(:D_, d), Expr(:call, :size, :A, d))
        push!(block.args, ex)
    end
    block
end

function sizeproductsblock(N::Int)
    block = Expr(:block)
    for k = 3:N
        ex = Expr(:(=), Symbol(:D_, ntuple(identity, k - 1)...),
                  Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)...))
        push!(block.args, ex)
    end
    block
end

function sumprodprecomputed2(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2))
            push!(e.args, ex)
        else
            ex = Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(:i_, k))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodconstant(N::Int)
    Expr(:(=), :D_sp,
         Expr(:call, :+, ntuple(d -> Expr(:call, :*, ntuple(i -> Symbol(:D_, i), d)...), N - 1)...))
end

function bouterloopgen(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] !== Val{1}
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    Expr(:for, block)
end
function binnerloopgen(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] === Val{1}
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    Expr(:for, block)
end

function binnerpost(N::Int, D)
    params = D.parameters
    b = Expr(:ref, :B, ntuple(d -> params[d] === Val{1} ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Val{1} ? 1 : Symbol(:i_, d), N)...)
    Expr(:block, Expr(:(=), b, :m), Expr(:(=), c, :j))
end

function compareblock(f, N::Int)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, Symbol(f), a, :m)) # f should only be > or <
    # mₑ = Expr(:(=), :m, Expr(:if, :y, a, :m))
    # jₑ = Expr(:(=), :j, Expr(:if, :y, d, :j))
    mₑ = Expr(:(=), :m, Expr(:call, :ifelse, :y, a, :m))
    jₑ = Expr(:(=), :j, Expr(:call, :ifelse, :y, d, :j))
    Expr(:block, yₑ, mₑ, jₑ)
end

function compareblock0(f, N::Int, D)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] === Val{1} ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Val{1} ? 1 : Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, Symbol(f), a, b)) # f should only be > or <
    bₑ = Expr(:(=), b, Expr(:call, :ifelse, :y, a, b))
    cₑ = Expr(:(=), c, Expr(:call, :ifelse, :y, d, c))
    Expr(:block, yₑ, bₑ, cₑ)
end

function findmax_quote0(N::Int, D)
    loops = loopgen(N)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    b4 = compareblock0(>, N, D)
    push!(loops.args, b4)
    return quote
        $b1
        $b2
        $b3
        $loops
    end
end
@generated function _bfindmax0!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                                A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D}
    findmax_quote0(N, D)
end
function bfindmax0(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
    Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
    B = fill(typemin(T), Dᴮ)
    C = similar(B, Int)
    _bfindmax0!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end

function findmax_quote(N::Int, D)
    block1 = sizeblock(N)
    block2 = sizeproductsblock(N)
    block3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outerloops = bouterloopgen(N, D)
    block4 = Expr(:block, Expr(:(=), :j, 1), Expr(:(=), :m, Expr(:call, :typemin, :T)))
    innerloops = binnerloopgen(N, D)
    block5 = compareblock(>, N)
    push!(innerloops.args, block5)
    push!(block4.args, innerloops)
    block6 = binnerpost(N, D)
    push!(block4.args, block6.args...)
    push!(outerloops.args, block4)
    return quote
        $block1
        $block2
        $block3
        $outerloops
    end
end
@generated function _bfindmax!(B::AbstractArray{T, N}, C::AbstractArray{Tₒ, N},
                               A::AbstractArray{T, N}, dims::D) where {T, Tₒ, N, D}
    findmax_quote(N, D)
end
function bfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : size(A, d), N)
    Dᴮ′ = ntuple(d -> d ∈ dims ? Val(1) : size(A, d), N)
    B = similar(A, Dᴮ)
    C = similar(B, Int)
    _bfindmax!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end

A = rand(4, 3, 5);
dims = (2,3)
using BenchmarkTools
@benchmark findmax(A, dims=dims)
@benchmark bfindmax(A, dims)
@benchmark bfindmax0(A, dims)
bfindmax0(A, dims) == findmax(A, dims=dims) == bfindmax(A, dims)

for d₂ = 1:ndims(A), d₁ = 1:ndims(A)
    dims = (d₁, d₂)
    @assert bfindmax(A, dims) == findmax(A, dims=dims)
end
DD = typeof(ntuple(d -> d ∈ dims ? Val(1) : size(A, d), ndims(A)))
findmax_quote0(ndims(A), DD)
findmax_quote(ndims(A), DD)
compareblock0(>, ndims(A), DD)

ext: jl


text:
using ROSESoss
import MeasureTheory as MT
import Distributions as Dists

obs = ehtim.obsdata.load_uvfits("/home/ptiede/Research/Projects/RefractiveScatteringThemis/TestDataGen/ER5/lo/hops_3601_M87+netcal.uvfits")
obs.add_scans()
obsavg = obs.avg_coherent(0.0, scan_avg=true)
obsavg.add_logcamp(debias=true, count="min")
obsavg.add_cphase(count="min-cut0bl")
obsavg.add_amp(debias=true)

dlca = ROSESoss.extract_lcamp(obsavg)
dcp = ROSESoss.extract_cphase(obsavg)
dva = ROSESoss.extract_amps(obsavg)

cm = create_joint(ROSESoss.mring(N=1,), dlca, dcp)
dy = DynestyStatic(nlive=200)
tvd, statsd = sample(dy, cm; print_progress=true)

cmva = create_joint(ROSESoss.mring(N=1,), dva, dcp)
tvd, statsd, sample(dy, cmva; print_progress=true)


text:
using ..Data
using MLDatasets
using CSV
using Pkg.Artifacts

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

const twenty_dataset_names = [
        "accidents", "ad", "baudio", "bbc", "bnetflix", "book", "c20ng", "cr52", "cwebkb",
        "dna", "jester", "kdd", "kosarek", "msnbc", "msweb", "nltcs", "plants", "pumsb_star", "tmovie", "tretail",
        "binarized_mnist"
];

const zoo_version = "/Circuit-Model-Zoo-0.1.2";

#####################
# Data loaders
#####################

function dataset(data; do_threshold=false, do_shuffle=false, batch_size=-1)
    shuffled_data = do_shuffle ? shuffle(data) : data
    thresholded_data = do_threshold ? threshold(shuffled_data) : shuffled_data
    batch_size > 0 ? batch(thresholded_data, batch_size) : thresholded_data
end

function mnist()
    # transposing makes slicing by variable much much faster
    # need to take a copy to physically move the data around
    train_x = copy_with_eltype(transpose(MNIST.convert2features(MNIST.traintensor())), Float32)
    test_x  = copy_with_eltype(transpose(MNIST.convert2features(MNIST.testtensor())), Float32)

    train_y::Vector{UInt8} = MNIST.trainlabels()
    test_y::Vector{UInt8}  = MNIST.testlabels()

    train = XYData(XData(train_x),train_y)
    valid = nothing
    test = XYData(XData(test_x),test_y)

    XYDataset(train,valid,test)
end

sampled_mnist() = twenty_datasets("binarized_mnist")

function twenty_datasets(name)
    @assert in(name, twenty_dataset_names)
    data_dir = artifact"density_estimation_datasets"
    function load(type)
        dataframe = CSV.read(data_dir*"/Density-Estimation-Datasets-1.0.1/datasets/$name/$name.$type.data"; header=false,
                 truestrings=["1"], falsestrings=["0"], type=Bool, strict=true)
        XData(BitArray(Base.convert(Matrix{Bool}, dataframe)))
    end
    train = load("train")
    valid = load("valid")
    test = load("test")
    XDataset(train,valid,test)
end

#####################
# Circuit loaders
#####################

zoo_vtree_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/vtrees/$name"

zoo_cnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/cnfs/$name"

zoo_dnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/dnfs/$name"

zoo_lc_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/lcs/$name"

zoo_clt_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/clts/$name"

zoo_psdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/psdds/$name"

zoo_sdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/sdds/$name"

# loaders

zoo_vtree(name) =
    load_vtree(zoo_vtree_file(name))

zoo_dnf(name) =
    load_dnf(zoo_dnf_file(name))

zoo_cnf(name) =
    load_cnf(zoo_cnf_file(name))


text:
using SignalDecomposition, Random, Statistics
m = Sinusoidal([1, 2] ./ 2π)
Random.seed!(151512)
tu = cumsum(rand(L)/4) # non-equispaced timeaxis
pu = periodicf(tu)

@testset "Sinusoidal-2" begin
    for (name, re) in zip(("cos(3t)", "gaussian"), (cos.(3tu), noise))
        s = pu .+ re .- mean(re)
        x, r = decompose(tu, s, m)
        errper = nrmse(pu, x)
        @test errper < 0.1
        errres = nrmse(re, r)
        @test nrmse(r, re) < 0.1
        # println("  "*name*" errper=$errper, errres=$errres")
        errori = nrmse(s, x .+ r)
        @test errori < 1e-15

        # figure()
        # ax1 = subplot(211)
        # plot(pu; alpha = 0.75, label = "x component")
        # plot(x; alpha = 0.75, ls = "dashed", label = "err=$errper")
        # ylabel("x"); legend()
        # subplot(212; sharex = ax1)
        # plot(re; alpha = 0.75, label = "r component")
        # plot(r; alpha = 0.75, ls = "dashed", label = "err=$errres")
        # ylabel("r"); legend()
    end
end

@testset "Sinusoidal-3" begin
    s = pu .+ cos.(3tu) .+ noise
    x, r = decompose(tu, s, Sinusoidal([1, 2, 3]./2π))
    errper = nrmse(pu .+ cos.(3tu), x)
    @test errper < 0.1
    errres = nrmse(noise, r)
    @test errres < 0.1
    errori = nrmse(s, x .+ r)
    @test errori < 1e-15

    # figure()
    # ax1 = subplot(211)
    # plot(pu; alpha = 0.75, label = "x component")
    # plot(x; alpha = 0.75, ls = "dashed", label = "err=$errper")
    # ylabel("x"); legend()
    # subplot(212; sharex = ax1)
    # plot(noise; alpha = 0.75, label = "r component")
    # plot(r; alpha = 0.75, ls = "dashed", label = "err=$errres")
    # ylabel("r"); legend()
end


text:
using CSV, DataFrames, DelimitedFiles, Statistics

const minimum_similarity = 0.8

function setup_data()
  movies = readdlm("top_10_movies_user_rankings.tsv", '\t')
  movies = permutedims(movies, (2,1))
  movies = convert(DataFrame, movies)

  names = convert(Array, movies[1, :])[1,:]

  names!(movies, [Symbol(name) for name in names])
  deleterows!(movies, 1)
  rename!(movies, [Symbol("Movie title") => :User])
end

function movie_similarity(target_movie)
  similarity = Dict{Symbol,Float64}()

  for movie in names(movies[:, 2:end])
    movie == target_movie && continue
    ratings = movies[:, [movie, target_movie]]

    common_users = ratings[(ratings[movie] .>= 0) .& (ratings[target_movie] .> 0), :]

    correlation = try
      cor(common_users[movie], common_users[target_movie])
    catch
      0.0
    end

    similarity[movie] = correlation
  end

  # println("The movie $target_movie is similar to $similarity")
  similarity
end

function recommendations(target_movie)
  recommended = Dict{String,Vector{Tuple{String,Float64}}}()

  # @show target_movie
  # @show movie_similarity(target_movie)

  for (movie, similarity) in movie_similarity(target_movie)
    movie == target_movie && continue
    similarity > minimum_similarity || continue

    # println("Checking to which users we can recommend $movie")

    recommended["$movie"] = Vector{Tuple{String,Float64}}()

    for user_row in eachrow(movies)
      if user_row[target_movie] >= 5
        # println("$(user_row[:User]) has watched $target_movie so we can recommend similar movies")

        if user_row[movie] == 0
          # println("$(user_row[:User]) has not watched $movie so we can recommend it")
          # println("Recommending $(user_row[:User]) the movie $movie")

          push!(recommended["$movie"], (user_row[:User], user_row[target_movie] * similarity))
        end
      end
    end
  end

  recommended
end

const movies = setup_data()
println("Recommendations for users that watched Finding Dory (2016): $(recommendations(Symbol("Finding Dory (2016)")))")

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
"""
Linear interpolation in one dimension

##### Fields

- `breaks::AbstractVector` : A sorted array of grid points on which to interpolate

- `vals::AbstractVector` : The function values associated with each of the grid points

##### Examples

```julia
breaks = cumsum(0.1 .* rand(20))
vals = 0.1 .* sin.(breaks)
li = LinInterp(breaks, vals)

# do interpolation via `call` method on a LinInterp object
li(0.2)

# use broadcasting to evaluate at multiple points
li.([0.1, 0.2, 0.3])
```

"""
struct LinInterp{TV<:AbstractArray,TB<:AbstractVector}
    breaks::TB
    vals::TV
    _n::Int
    _ncol::Int

    function LinInterp{TV,TB}(b::TB, v::TV) where {TB,TV}
        if size(b, 1) != size(v, 1)
            m = "breaks and vals must have same number of elements"
            throw(DimensionMismatch(m))
        end

        if !issorted(b)
            m = "breaks must be sorted"
            throw(ArgumentError(m))
        end
        new{TV,TB}(b, v, length(b), size(v, 2))
    end
end

function Base.:(==)(li1::LinInterp, li2::LinInterp)
    all(getfield(li1, f) == getfield(li2, f) for f in fieldnames(typeof(li1)))
end

function LinInterp(b::TB, v::TV) where {TV<:AbstractArray,TB<:AbstractVector}
    LinInterp{TV,TB}(b, v)
end

function (li::LinInterp{<:AbstractVector})(xp::Number)
    ix = searchsortedfirst(li.breaks, xp)

    # handle corner cases
    @inbounds begin
        ix == 1 && return li.vals[1]
        ix == li._n + 1 && return li.vals[end]

        # now get on to the real work...
        z = (li.breaks[ix] - xp)/(li.breaks[ix] - li.breaks[ix-1])

        return (1-z) * li.vals[ix] + z * li.vals[ix-1]
    end
end

function (li::LinInterp{<:AbstractMatrix})(xp::Number, col::Int)
    ix = searchsortedfirst(li.breaks, xp)
    @boundscheck begin
        if col > li._ncol || col < 1
            msg = "col must be beteween 1 and $(li._ncol), found $col"
            throw(BoundsError(msg))
        end
    end

    @inbounds begin
        # handle corner cases
        ix == 1 && return li.vals[1, col]
        ix == li._n + 1 && return li.vals[end, col]

        # now get on to the real work...
        z = (li.breaks[ix] - xp)/(li.breaks[ix] - li.breaks[ix-1])

        return (1-z) * li.vals[ix, col] + z * li.vals[ix-1, col]
    end
end

_out_eltype(li::LinInterp{TV,TB}) where {TV,TB} = promote_type(eltype(TV), eltype(TB))

function (li::LinInterp{<:AbstractMatrix})(
        xp::Number, cols::AbstractVector{<:Integer}
    )
    ix = searchsortedfirst(li.breaks, xp)
    @boundscheck begin
        for col in cols
            if col > li._ncol || col < 1
                msg = "all cols must be beteween 1 and $(li._ncol), found $col"
                throw(BoundsError(msg))
            end
        end
    end

    out = Vector{_out_eltype(li)}(undef, length(cols))

    @inbounds begin
        # handle corner cases
        if ix == 1
            for (ind, col) in enumerate(cols)
                out[ind] = li.vals[1, col]
            end
            return out
        end

        if ix == li._n + 1
            for (ind, col) in enumerate(cols)
                out[ind] = li.vals[end, col]
            end
            return out
        end

        # now get on to the real work...
        z = (li.breaks[ix] - xp)/(li.breaks[ix] - li.breaks[ix-1])

        for (ind, col) in enumerate(cols)
            out[ind] = (1-z) * li.vals[ix, col] + z * li.vals[ix-1, col]
        end

        return out
    end
end

(li::LinInterp{<:AbstractMatrix})(xp::Number) = li(xp, 1:li._ncol)

"""
    interp(grid::AbstractVector, function_vals::AbstractVector)

Linear interpolation in one dimension

##### Examples

```julia
breaks = cumsum(0.1 .* rand(20))
vals = 0.1 .* sin.(breaks)
li = interp(breaks, vals)

# Do interpolation by treating `li` as a function you can pass scalars to
li(0.2)

# use broadcasting to evaluate at multiple points
li.([0.1, 0.2, 0.3])
```
"""
function interp(grid::AbstractVector, function_vals::AbstractVector)
    if !issorted(grid)
        inds = sortperm(grid)
        return LinInterp(grid[inds], function_vals[inds])
    else
        return LinInterp(grid, function_vals)
    end
end

function interp(grid::AbstractVector, function_vals::AbstractMatrix)
    if !issorted(grid)
        inds = sortperm(grid)
        return LinInterp(grid[inds], function_vals[inds, :])
    else
        return LinInterp(grid, function_vals)
    end
end


text:
using LinearAlgebra, CUDA, NNlib, NNlibCUDA, BatchedBLAS, SymmetricFormats, BenchmarkTools, DataFrames, Gadfly, JLD2

macro belapsed_median(args...)
    return esc(quote
        $BenchmarkTools.time($BenchmarkTools.median($BenchmarkTools.@benchmark $(args...)))/1e9
    end)
end

function doit(L,N)
    x2 = CuArray(rand(L,N));
    x3 = CuArray(rand(L,1,N));
    y2 = CuArray(rand(L,N));
    y3 = CuArray(rand(L,1,N));
    o1 = CuArray(rand(N));
    o3 = CuArray(rand(1,1,N));

    tbgemm = @belapsed_median CUDA.@sync batched_mul!($o3, batched_transpose($x3), $y3)

    tbdot = @belapsed_median CUDA.@sync batched_dot!($o1, $x2, $y2)

    CUDA.unsafe_free!.((x2, x3, y2, y3, o1, o3))
    CUDA.memory_status();  println("\n")

    return tbgemm, tbdot
end

df_b_dot = DataFrame(func = String[], N = Int[], L = Int[], elapsed_time = Float64[])

N=32768
for L=[32, 64, 128, 256]
    times = doit(L, N)
    push!(df_b_dot, ("bgemm", N, L, times[1]))
    push!(df_b_dot, ("bdot",  N, L, times[2]))
end

L=128
for N=[8192, 16384, 65536, 131072]
    times = doit(L, N)
    push!(df_b_dot, ("bgemm", N, L, times[1]))
    push!(df_b_dot, ("bdot",  N, L, times[2]))
end

function doit(L,N)
    _A = rand(L,L,N);
    A3 = CuArray(_A);
    x2 = CuArray(rand(L,N));
    x3 = CuArray(rand(L,1,N));
    y2 = CuArray(rand(L,N));
    y3 = CuArray(rand(L,1,N));

    tbgemm = @belapsed_median CUDA.@sync batched_mul!($y3, $A3, $x3)

    tbgemvn = @belapsed_median CUDA.@sync batched_gemv!('N', 1.0, $A3, $x2, 0.0, $y2)
    tbgemvt = @belapsed_median CUDA.@sync batched_gemv!('T', 1.0, $A3, $x2, 0.0, $y2)

    tbsymvu = @belapsed_median CUDA.@sync batched_symv!('U', 1.0, $A3, $x2, 0.0, $y2)
    tbsymvl = @belapsed_median CUDA.@sync batched_symv!('L', 1.0, $A3, $x2, 0.0, $y2)

    AP = CuArray(hcat([SymmetricPacked(x, :U).tri for x in eachslice(_A, dims=3)]...));
    tbspmvu = @belapsed_median CUDA.@sync batched_spmv!('U', 1.0, $AP, $x2, 0.0, $y2)
    AP = CuArray(hcat([SymmetricPacked(x, :L).tri for x in eachslice(_A, dims=3)]...));
    tbspmvl = @belapsed_median CUDA.@sync batched_spmv!('L', 1.0, $AP, $x2, 0.0, $y2)

    CUDA.unsafe_free!.((A3, AP, x2, x3, y2, y3))
    CUDA.memory_status();  println("\n")

    return tbgemm, tbgemvn, tbgemvt, tbsymvu, tbsymvl, tbspmvu, tbspmvl
end

df_b_mv = DataFrame(func = String[], N = Int[], L = Int[], elapsed_time = Float64[])

N=32768
for L=[32, 64, 128, 256]
    times = doit(L, N)
    push!(df_b_mv, ("bgemm",  N, L, times[1]))
    push!(df_b_mv, ("bgemvn", N, L, times[2]))
    push!(df_b_mv, ("bgemvt", N, L, times[3]))
    push!(df_b_mv, ("bsymvu", N, L, times[4]))
    push!(df_b_mv, ("bsymvl", N, L, times[5]))
    push!(df_b_mv, ("bspmvu", N, L, times[6]))
    push!(df_b_mv, ("bspmvl", N, L, times[7]))
end

L=128
for N=[8192, 16384, 65536, 131072]
    times = doit(L, N)
    push!(df_b_mv, ("bgemm",  N, L, times[1]))
    push!(df_b_mv, ("bgemvn", N, L, times[2]))
    push!(df_b_mv, ("bgemvt", N, L, times[3]))
    push!(df_b_mv, ("bsymvu", N, L, times[4]))
    push!(df_b_mv, ("bsymvl", N, L, times[5]))
    push!(df_b_mv, ("bspmvu", N, L, times[6]))
    push!(df_b_mv, ("bspmvl", N, L, times[7]))
end

function doit(L,N)
    _A = rand(L,L,N);
    A3 = CuArray(_A);
    x2 = CuArray(rand(L,N));
    x3 = CuArray(rand(L,1,N));
    y2 = CuArray(rand(L,N));
    y3 = CuArray(rand(L,1,N));

    tbgemm = @belapsed_median CUDA.@sync batched_mul!($A3, $x3, batched_transpose($x3), -1.0, 1.0)

    tbger = @belapsed_median CUDA.@sync batched_ger!(-1.0, $x2, $y2, $A3)

    tbsyru = @belapsed_median CUDA.@sync batched_syr!('U', -1.0, $x2, $A3)
    tbsyrl = @belapsed_median CUDA.@sync batched_syr!('L', -1.0, $x2, $A3)

    AP = CuArray(hcat([SymmetricPacked(x, :U).tri for x in eachslice(_A, dims=3)]...));
    tbspru = @belapsed_median CUDA.@sync batched_spr!('U', -1.0, $x2, $AP)
    AP = CuArray(hcat([SymmetricPacked(x, :L).tri for x in eachslice(_A, dims=3)]...));
    tbsprl = @belapsed_median CUDA.@sync batched_spr!('L', -1.0, $x2, $AP)

    CUDA.unsafe_free!.((A3, AP, x2, x3, y2, y3))
    CUDA.memory_status();  println("\n")

    return tbgemm, tbger, tbsyru, tbsyrl, tbspru, tbsprl
end

df_b_r = DataFrame(func = String[], N = Int[], L = Int[], elapsed_time = Float64[])

N=32768
for L=[32, 64, 128, 256]
    times = doit(L, N)
    push!(df_b_r, ("bgemm", N, L, times[1]))
    push!(df_b_r, ("bger",  N, L, times[2]))
    push!(df_b_r, ("bsyru", N, L, times[3]))
    push!(df_b_r, ("bsyrl", N, L, times[4]))
    push!(df_b_r, ("bspru", N, L, times[5]))
    push!(df_b_r, ("bsprl", N, L, times[6]))
end

L=128
for N=[8192, 16384, 65536, 131072]
    times = doit(L, N)
    push!(df_b_r, ("bgemm", N, L, times[1]))
    push!(df_b_r, ("bger",  N, L, times[2]))
    push!(df_b_r, ("bsyru", N, L, times[3]))
    push!(df_b_r, ("bsyrl", N, L, times[4]))
    push!(df_b_r, ("bspru", N, L, times[5]))
    push!(df_b_r, ("bsprl", N, L, times[6]))
end

save("bench.jld2", "df_b_dot", df_b_dot, "df_b_mv", df_b_mv, "df_b_r", df_b_r)
#=
df_b_dot = load("bench.jld2", "df_b_dot")
df_b_mv = load("bench.jld2", "df_b_mv")
df_b_r = load("bench.jld2", "df_b_r")
=#

theme = Theme(highlight_width=0mm, point_size=0.6mm, line_style=[:solid, :dot]);

gd = groupby(df_b_dot, [:N, :L])
df_b_dot2 = combine(gd, :func => identity, :elapsed_time => x->x./x[1])

Nbdotplot = plot(df_b_dot2[df_b_dot2.L.==2^7, :],
             x=:N, y=:elapsed_time_function, color=:func_identity,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.colorkey(title=""),
             Guide.title("L=27"),
             Guide.xlabel("N", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

Lbdotplot = plot(df_b_dot2[df_b_dot2.N.==2^15, :],
             x=:L, y=:elapsed_time_function, color=:func_identity,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.xticks(orientation=:horizontal),
             Guide.colorkey(title=""),
             Guide.title("N=215"),
             Guide.xlabel("L", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

df_b_mv[!,:func2] .= [x[end] in ['t', 'u'] ? 2 : 1 for x in df_b_mv[!,:func]]
df_b_mv[!,:func] .=  [rstrip(x, ['u', 'l', 'n', 't']) for x in df_b_mv[!,:func]]

gd = groupby(df_b_mv, [:N, :L])
df_b_mv2 = transform(gd, :elapsed_time => x->x./x[1])

Nbmvplot = plot(df_b_mv2[df_b_mv2.L.==2^7, :],
             x=:N, y=:elapsed_time_function, color=:func, linestyle=:func2,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.colorkey(title=""),
             Guide.title("L=27"),
             Guide.xlabel("N", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

Lbmvplot = plot(df_b_mv2[df_b_mv2.N.==2^15, :],
             x=:L, y=:elapsed_time_function, color=:func, linestyle=:func2,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.xticks(orientation=:horizontal),
             Guide.colorkey(title=""),
             Guide.title("N=215"),
             Guide.xlabel("L", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

df_b_r[!,:func2] .= [x[end]=='u' ? 2 : 1 for x in df_b_r[!,:func]]
df_b_r[!,:func] .=  [rstrip(x, ['u', 'l']) for x in df_b_r[!,:func]]

gd = groupby(df_b_r, [:N, :L])
df_b_r2 = transform(gd, :elapsed_time => x->x./x[1])

Nbrplot = plot(df_b_r2[df_b_r2.L.==2^7, :],
             x=:N, y=:elapsed_time_function, color=:func, linestyle=:func2,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.colorkey(title=""),
             Guide.title("L=27"),
             Guide.xlabel("N", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

Lbrplot = plot(df_b_r2[df_b_r2.N.==2^15, :],
             x=:L, y=:elapsed_time_function, color=:func, linestyle=:func2,
             Scale.x_log2,
             Geom.line, Geom.point,
             Guide.xticks(orientation=:horizontal),
             Guide.colorkey(title=""),
             Guide.title("N=215"),
             Guide.xlabel("L", orientation=:horizontal),
             Guide.ylabel("time re. bgemm", orientation=:vertical),
             theme);

gridstack([Nbdotplot Nbmvplot Nbrplot;
           Lbdotplot Lbmvplot Lbrplot]) |> SVG("bench.svg", 21cm, 12cm)


text:
function run_OnDiskMonitor_tests()
    n = 1000
    m = SimTools.OnDiskMonitor(@__DIR__, "monitor_test.jld2"; dt = 1.0)
    SimTools.reset!(m)
    SimTools.clear_file(m)
    @assert isempty(m.cache)
    @assert !isfile(m.file)

    # Watch task (It can be in another process)
    tests_results = []
    cvals = []
    @async SimTools.watch(m) do ddat
        isempty(ddat) && return
        check = true
        dvals = ddat[:vals]
        for dval in dvals
            check &= dval in cvals
        end
        push!(tests_results, check)
    end

    # Fake job to monitor
    for i in 1:n
        push!(cvals, i)
        SimTools.record!(m) do dat
            dat[:vals] = cvals
        end
        sleep(5.0/n)
    end
    SimTools.reset!(m)
    sleep(1.0)

    @test isempty(m.cache)
    @test length(cvals) == n
    @test all(tests_results)

    SimTools.clear_file(m)
    @test !isfile(m.file)
end
run_OnDiskMonitor_tests()


text:
immutable Trade
    start::DateTime
    finish::DateTime
    open::Float64
    close::Float64
    side::ASCIIString
    quantity::Int
    timeseries::FinancialTimeSeries
    function Trade(start, finish, open, close, side, quantity, timeseries)

                   if timeseries.timestamp[1] > start || timeseries.timestamp[end] < finish
                       error("timeseries doesn't match trade date range")
                   else
                       sd    = findfirst(timeseries.timestamp .== start)
                       fd    = findfirst(timeseries.timestamp .== finish)
                   end
                   new(start, finish, open, close, side, quantity, timeseries[sd:fd])
    end
end

function Trade(b::Blotter, f::FinancialTimeSeries)
    Trade(b.timestamp[1], b.timestamp[length(b.timestamp)], b.values[length(b.timestamp)+1], b.values[end], "long", b.values[1], f)
end

########## get Trades #########

# first get the index (row number) where cumsum == 0, or trades are closed

function closedtrades(b::Blotter)
    # get index of zeros first
    totals  = cumsum(b.values[:,1])
    counter = Int[]
    for d in 1:length(b.timestamp)
        if totals[d] == 0
            push!(counter,d)
        end
    end
    counter
end

function startedtrades(b::Blotter)
    closed  = closedtrades(b)
    counter = [1]
    for n in 1:length(closed) - 1
        push!(counter, closed[n]+1)
    end
    counter
end

function tradearray(b::Blotter, fts::FinancialTimeSeries)
    started  = startedtrades(b)
    closed   = closedtrades(b)
    duration = hcat(started, closed)
    trades   = Trade[]
    for d in 1:size(duration,1)
        s = duration[d,1]
        f = duration[d,2]
        sdate = b.timestamp[s]
        sd    = findfirst(fts.timestamp .== sdate)
        fdate = b.timestamp[f]
        fd    = findfirst(fts.timestamp .== fdate)
        push!(trades, Trade(b[s:f], fts[sd:fd]))
    end
    trades
end

function tradestats(t::Vector{Trade})
#     symbol             = t.[1].timeseries.instrument.ticker
#     txns_total         = length(trades)
#     trades_total       =
#     pnls               = convert(Array{Float64}, ([(t.close-t.open) for t in trades].*100))
#     pos                = pnls[pnls.>0]
#     neg                = pnls[pnls.<0]
#     netprofit          =
#     trade_pl_mean      =
#     med_Trade_PL       =
#     largest_winner     = maximum(pos)
#     largest_loser      = minimum(neg)
#     std_Dev_Trade_PL   =
#     percent_positive   = length(pos)/length(trades)
#     percent_negative   = length(neg)/length(trades)
#     profitfactor       = sum(pos)/abs(sum(neg))
#     avg_Win_Trade      = mean(pos)
#     med_Win_Trade      = median(pos)
#     avg_Losing_Trade   = mean(neg)
#     med_Losing_Trade   = median(neg)
#     avg_Daily_PL       =
#     med_Daily_PL       =
#     std_Dev_Daily_PL   =
#     ann_Sharpe         =
#     max_drawdown       =
#     profit_To_Max_Draw =
#     avg_WinLoss_Ratio  =
#     med_WinLoss_Ratio  =
#     max_equity         =
#     min_equity         =
#     end_equity         =
#     mae                = convert(Array{Float64}, ([t.open - minimum(t.timeseries["Low"].values) for t in trades]))
end


text:
using Test
using JET

using IntervalSets

import ClimaCore: Domains, Meshes, Topologies, Spaces, Fields, Operators
import ClimaCore.Domains: Geometry

# We need to pull these broadcasted expressions out as
# toplevel functions due to how broadcast expressions are
# lowered so JETTest can have a single callsite to analyze.

function opt_InterpolateF2C(face_field)
    I = Operators.InterpolateF2C()
    return I.(identity.(face_field))
end

function opt_WeightedInterpolateF2C(weights, face_field)
    WI = Operators.WeightedInterpolateF2C()
    return identity.(WI.(weights, face_field))
end

function opt_LeftBiasedF2C(face_field)
    LB = Operators.LeftBiasedF2C(left = Operators.SetValue(0.0))
    return LB.(identity.(face_field))
end

function opt_RightBiasedF2C(face_field)
    RB = Operators.RightBiasedF2C(right = Operators.SetValue(0.0))
    return RB.(identity.(face_field))
end

function opt_AdvectionF2F(face_vel, face_field)
    A = Operators.AdvectionF2F()
    return A.(face_vel, identity.(face_field))
end

function opt_FluxCorrectionF2F_Extrapolate(face_vel, face_field)
    FC = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return FC.(face_vel, identity.(face_field))
end

function opt_GradientF2C(face_field)
    ∇ᶜ = Operators.GradientF2C()
    return Geometry.WVector.(∇ᶜ.(sin.(face_field)))
end

function opt_DivergenceF2C(face_field)
    divᶜ = Operators.DivergenceF2C()
    return divᶜ.(Geometry.WVector.(sin.(face_field)))
end

function opt_SetBoundary_SetValue(face_field)
    B = Operators.SetBoundaryOperator(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return B.(sin.(face_field))
end

function opt_InterpolateC2F_SetValue(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return I.(identity.(center_field))
end

function opt_InterpolateC2F_SetGradient(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return I.(identity.(center_field))
end

function opt_InterpolateC2F_Extrapolate(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return I.(identity.(center_field))
end

function opt_WeightedInterpolateC2F_SetValue(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return identity.(WI.(weights, center_field))
end

function opt_WeightedInterpolateC2F_SetGradient(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return identity.(WI.(weights, center_field))
end

function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return identity.(WI.(weights, center_field))
end

function opt_LeftBiasedC2F(center_field)
    LB = Operators.LeftBiasedC2F(left = Operators.SetValue(0.0))
    return LB.(identity.(center_field))
end

function opt_RightBiasedC2F(center_field)
    RB = Operators.RightBiasedC2F(right = Operators.SetValue(0.0))
    return RB.(identity.(center_field))
end

function opt_UpwindBiasedProductC2F_SetValue(face_vel, center_field)
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return UB.(face_vel, identity.(center_field))
end

function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field)
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return UB.(face_vel, identity.(center_field))
end

function opt_AdvectionC2C_SetValue(face_vel, center_field)
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return A.(face_vel, identity.(center_field))
end

function opt_AdvectionC2C_Extrapolate(face_vel, center_field)
    A = Operators.AdvectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return A.(face_vel, identity.(center_field))
end

function opt_FluxCorrectionC2C_Extrapolate(face_vel, center_field)
    FC = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return FC.(face_vel, identity.(center_field))
end

function opt_GradientC2F_SetValue(center_field)
    ∇ᶠ = Operators.GradientC2F(
        left = Operators.SetValue(1.0),
        right = Operators.SetValue(-1.0),
    )
    return Geometry.WVector.(∇ᶠ.(cos.(center_field)))
end

function opt_GradientC2F_SetGradient(center_field)
    ∇ᶠ = Operators.GradientC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return Geometry.WVector.(∇ᶠ.(cos.(center_field)))
end

function opt_DivergenceC2F_SetValue(center_field)
    divᶠ = Operators.DivergenceC2F(
        left = Operators.SetValue(Geometry.WVector(0.0)),
        right = Operators.SetValue(Geometry.WVector(0.0)),
    )
    return divᶠ.(Geometry.WVector.(sin.(center_field)))
end

function opt_DivergenceC2F_SetDivergence(center_field)
    # DivergenceC2F, SetDivergence
    divᶠ = Operators.DivergenceC2F(
        left = Operators.SetDivergence(0.0),
        right = Operators.SetDivergence(0.0),
    )
    return divᶠ.(Geometry.WVector.(cos.(center_field)))
end

function opt_CurlC2F_SetValue(center_field)
    # DivergenceC2F, SetDivergence
    curlᶠ = Operators.CurlC2F(
        left = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
        right = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
    )
    return curlᶠ.(Geometry.Covariant1Vector.(cos.(center_field)))
end

function hspace1d(FT)
    hdomain = Domains.IntervalDomain(
        Geometry.XPoint{FT}(-pi) .. Geometry.XPoint{FT}(pi),
        periodic = true,
    )
    hmesh = Meshes.IntervalMesh(hdomain, nelems = 3)
    htopology = Topologies.IntervalTopology(hmesh)
    Nq = 3
    quad = Spaces.Quadratures.GLL{Nq}()
    return Spaces.SpectralElementSpace1D(htopology, quad)
end

function hspace2d(FT)
    hdomain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-pi) .. Geometry.XPoint{FT}(pi),
        Geometry.YPoint{FT}(-pi) .. Geometry.YPoint{FT}(pi);
        x1periodic = true,
        x2periodic = true,
    )
    Nq = 3
    quad = Spaces.Quadratures.GLL{Nq}()
    hmesh = Meshes.RectilinearMesh(hdomain, 3, 3)
    htopology = Topologies.Topology2D(hmesh)
    return Spaces.SpectralElementSpace2D(htopology, quad)
end

@static if @isdefined(var"@test_opt")
    @testset "Scalar Field ExtrudedFiniteDifferenceSpace" begin
        for FT in (Float64,), hspace in (hspace1d(FT), hspace2d(FT))
            vdomain = Domains.IntervalDomain(
                Geometry.ZPoint{FT}(0.0),
                Geometry.ZPoint{FT}(pi);
                boundary_tags = (:left, :right),
            )
            vmesh = Meshes.IntervalMesh(vdomain; nelems = 16)
            cvspace = Spaces.CenterFiniteDifferenceSpace(vmesh)

            center_space = Spaces.ExtrudedFiniteDifferenceSpace(hspace, cvspace)
            face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(center_space)

            faces = getproperty(Fields.coordinate_field(face_space), :z)
            face_values = ones(FT, face_space)
            face_velocities = Geometry.WVector.(face_values)

            centers = getproperty(Fields.coordinate_field(center_space), :z)
            center_values = ones(FT, center_space)
            center_velocities = Geometry.WVector.(center_values)

            # face space operators
            @test_opt sum(ones(FT, face_space))
            @test_opt sum(sin.(faces))

            @test_opt opt_InterpolateF2C(faces)
            @test_opt opt_WeightedInterpolateF2C(face_values, faces)

            @test_opt opt_LeftBiasedF2C(faces)
            @test_opt opt_RightBiasedF2C(faces)

            @test_opt opt_AdvectionF2F(face_velocities, faces)

            @test_opt opt_FluxCorrectionF2F_Extrapolate(
                center_velocities,
                faces,
            )

            @test_opt opt_GradientF2C(faces)
            @test_opt opt_GradientF2C(faces)
            @test_opt opt_DivergenceF2C(faces)

            @test_opt opt_SetBoundary_SetValue(faces)

            # center space operators
            @test_opt sum(ones(FT, center_space))
            @test_opt sum(sin.(centers))

            @test_opt opt_InterpolateC2F_SetValue(centers)
            @test_opt opt_InterpolateC2F_SetGradient(centers)
            @test_opt opt_InterpolateC2F_Extrapolate(centers)

            @test_opt opt_WeightedInterpolateC2F_SetValue(
                center_values,
                centers,
            )
            @test_opt opt_WeightedInterpolateC2F_SetGradient(
                center_values,
                centers,
            )
            @test_opt opt_WeightedInterpolateC2F_Extrapolate(
                center_values,
                centers,
            )

            @test_opt opt_LeftBiasedC2F(centers)
            @test_opt opt_RightBiasedC2F(centers)

            @test_opt opt_UpwindBiasedProductC2F_SetValue(
                face_velocities,
                centers,
            )
            @test_opt opt_UpwindBiasedProductC2F_Extrapolate(
                face_velocities,
                centers,
            )

            @test_opt opt_AdvectionC2C_SetValue(face_velocities, centers)
            @test_opt opt_AdvectionC2C_Extrapolate(face_velocities, centers)

            @test_opt opt_FluxCorrectionC2C_Extrapolate(
                face_velocities,
                centers,
            )

            @test_opt opt_GradientC2F_SetValue(centers)
            @test_opt opt_GradientC2F_SetGradient(centers)

            @test_opt opt_DivergenceC2F_SetValue(centers)
            @test_opt opt_DivergenceC2F_SetDivergence(centers)
            @test_opt opt_CurlC2F_SetValue(centers)
        end
    end

end


text:
# using Pkg
# Pkg.add(PackageSpec(url="https://github.com/amdeld/LibFEM.jl"))
using LibFEM
# Units system mm/tonne/s/K
# ===============================================================================
# 2D FRAME
# ===============================================================================
# (2)---(3)
#       /
#      /
#     /
#    /
#   /
#  /
# (1)
# ===============================================================================
# The element used here is a linear 2D straight truss with constant cross section
# The degrees of freedom are the u & v displacements
# ===============================================================================
# PARAMETERS
const L = 10000. # length in mm
const A = 100. # cross-sectional area in mm^2
const E = 210000. # modulus of elasticity in MPa [steel]
const FM = 10000. # force modulus in N
# ===============================================PRE-PROCESSING==================
# DEFINING AND DISCRETIZING[MESHING] THE STRUCTURE
# connectivity table
# elt||node_i||node_j
# 1|1|3
# 2|2|3
# grid()
X1pos = 0.;Y1pos = 0.
X2pos = 0.;Y2pos = L
X3pos = L;Y3pos = L
# lengths
L1 = d2_truss_elementlength(X1pos, Y1pos, X3pos, Y3pos) # length of element 1
L2 = d2_truss_elementlength(X2pos, Y2pos, X3pos, Y3pos) # length of element 2
# APPLYING GEOMETRIC&MATERIAL PROPERTIES
A1 = sqrt(2) * A # cross-sectional area of element 1
A2 = A # cross-sectional area of element 2
E1 = E # material of element 1
E2 = E # material of element 2
# writing-defining the element stiffness matrices
K1 = d2_truss_elementstiffness(E1, A1, L1, 45)
println("K1=\r")
display(K1)
K2 = d2_truss_elementstiffness(E2, A2, L2, 0)
println("K2=\r")
display(K2)
# ASSEMBLING THE GLOBAL STIFFNESS MATRIX
# matrices initialization
K = zeros(6, 6);K1P = zeros(6, 6)
# positionning stiffness matrices
K1P = d2_truss_assemble(K, K1, 1, 3)
println("K1P=\r")
display(K1P)
K = zeros(6, 6);K2P = zeros(6, 6)
K2P = d2_truss_assemble(K, K2, 2, 3)
println("K2P=\r")
display(K2P)
# assembling
K = K1P + K2P
println("K=\r")
display(K)
# SOLVING DISPACEMENT EQUATIONS
# extracting displacement submatrix via index vector
K_s = K[5:6,5:6]
# Setting-up the force subvector by applying Load & Boundary Conditions[LBC]]
F_s = [0, -FM]
# solving by gaussian elimination
U_s = K_s \ F_s
# SOLVING FORCE EQUATIONS
# setting-up the global nodal displacement vector
U = [0, 0, 0, 0, U_s[1],U_s[2]]
println("U=\r")
display(U)
# computing the global nodal force vector
F = K * U
println("F=\r")
display(F)
# COMPUTING STRAINS, INTERNAL FORCES, STRESSES
# writing the element nodal displacement vectors
U1 = [U[1], U[2], U[5], U[6]]
U2 = [U[3], U[4], U[5], U[6]]
# computing element strains
ϵ1 = d2_truss_elementstrain(L1, 45, U1)
@show ϵ1
ϵ2 = d2_truss_elementstrain(L2, 0, U2)
@show ϵ2
# computing element forces
f1 = d2_truss_elementforce(E1, A1, L1, 45, U1)
@show f1
f2 = d2_truss_elementforce(E2, A2, L2, 0, U2)
@show f2
# computing element stresses
σ1 = d2_truss_elementstress(E1, L1, 45, U1)
@show σ1
σ2 = d2_truss_elementstress(E2, L2, 0, U2)
@show σ2
# Plotting
using Plots
Xini = [X1pos, X3pos]
Yini = [Y1pos, Y3pos]
plot(Xini,Yini,label=false,marker=5,shape=:sq,color=:black)
Xini = [X2pos, X3pos]
Yini = [Y2pos, Y3pos]
plot!(Xini,Yini,label=false,marker=5, shape=:sq, color=:black)
fampl = 100
Xdef = [(X1pos + fampl * U[1]), (X3pos + fampl * U[5])]
Ydef = [(Y1pos + fampl * U[1]), (Y3pos + fampl * U[6])]
plot!(Xdef,Ydef,label=false,marker=5,shape=:sq,color=:blue)
Xdef = [(X2pos + fampl * U[3]), (X3pos + fampl * U[5])]
Ydef = [(Y2pos + fampl * U[4]), (Y3pos + fampl * U[6])]
plot!(Xdef,Ydef,label=false,title="Undeformed and Deformed Shapes",
xlabel="Horizontal coordinate X",
ylabel="Vertical coordinate Y",marker=5,shape=:sq,color=:blue)


text:
function filtereselection(v0, selection, buttons)
  if !isempty(buttons) && first(buttons) == 0  # if any button is pressed && its the left button
    selection #return diffed index
  else
    Vector2(-1)
  end
end

function haschanged(v0, selection)
  (v0[2] != selection, selection)
end
function edit(style::Style{:Default}, text::String, custumization::Dict{Symbol, Any})
  obj = visualize(text, style; custumization...)
  edit(obj[:text], obj)
end
function edit(style::Style{:Default}, textGPU::Texture{GLGlyph{Uint16}, 4, 2}, obj::RenderObject, custumization::Dict{Symbol, Any})
  screen = custumization[:screen]
  specialkeys = filteritems(screen.inputs[:buttonspressed], [GLFW.KEY_LEFT_CONTROL, GLFW.KEY_ENTER, GLFW.KEY_BACKSPACE], IntSet())

  selectiondata = lift(first, SELECTION[:mouse_hover])
  # Filter out the selected index,
  changed = lift(first, foldl(haschanged, (true, selectiondata.value), selectiondata))

  leftclick_selection = foldl(filtereselection, Vector2(-1), keepwhen(changed, Vector2(-1), selectiondata), screen.inputs[:mousebuttonspressed])
  glypharray = vec(data(textGPU))
  changed    = false
  v00        = (obj, obj.alluniforms[:textlength], textGPU, glypharray, leftclick_selection.value, changed)
  output     = foldl(edit_text, v00, leftclick_selection, screen.inputs[:unicodeinput], specialkeys)
  textsignal = lift(filter(v00, output) do x
    x[end] # <-has changed
  end) do x
    utf16(x[4][1:x[2]]) #x[4] <- glypharray, x[2] <- textlength
  end
  return (obj, textsignal)
end
function pressed{T<:Integer}(keys::Vector{T}, keyset)
    length(keyset) == length(keys) && all(keys) do x
        in(x, keyset)
    end
end
function edit_text(v0, selection1, unicode_keys, special_keys)
    # selection0 tracks, where the carsor is after a new character addition, selection10 tracks the old selection
    obj, textlength, textGPU, glypharray, selection0, changed = v0
    selected_object, selected_index = selection0
    selected_index += 1
    chars_added = 0
    edit_inbound = selected_index<=textlength+1 && selected_index>0 && selected_object == obj.id
    if selection1 != Vector2(-1) # if a new index was selected update selection
        return (obj, textlength, textGPU, glypharray, selection1, false)
    elseif edit_inbound && !(!isempty(unicode_keys) && IntSet(GLFW.MOD_SHIFT)==special_keys) && (!isempty(special_keys) || !isempty(unicode_keys))# something will get edited
        if !isempty(special_keys) && isempty(unicode_keys)
            if in(GLFW.KEY_BACKSPACE, special_keys) && selected_index>1
                splice!(glypharray, selected_index-1)
                chars_added = -1
            elseif in(GLFW.KEY_ENTER, special_keys)
                insert!(glypharray, selected_index, GLGlyph('\n', 0,0,0))
                chars_added = 1
            elseif pressed([GLFW.KEY_LEFT_CONTROL, GLFW.KEY_V], special_keys)
                p = clipboard()
                pasted = [GLGlyph(c,0,0,0) for c in p]
                if selected_index == 1
                  glypharray = [pasted, glypharray]
                elseif selected_index == textlength+1
                  glypharray = [glypharray, pasted]
                elseif selected_index>1 && selected_index>=textlength
                  glypharray = [sub(glypharray, 1:selected_index), pasted, sub(glypharray, selected_index:length(glypharray))]
                else
                  error("invalid text cursor index. Index: ", selected_index)
                end
                chars_added = length(p)
            end
        elseif edit_inbound && !isempty(unicode_keys) && (isempty(special_keys) || IntSet(GLFW.MOD_SHIFT)==special_keys)# else unicode input must have occured
            insert!(glypharray, selected_index, GLGlyph(first(unicode_keys), 0,0,0))
            chars_added = 1
        end
        if chars_added != 0
            textlength  += chars_added
            newselection = chars_added + selected_index
            selection0   = Vector2(selected_object, max(min(newselection, textlength+1), 1)-1)
            update_glyphpositions!(glypharray)
            if textlength > length(textGPU)
              resize!(textGPU, [size(textGPU,1), size(textGPU,2)*2])
            end
            remaining = div(textlength, 1024)
            if remaining < 1
                textGPU[1:textlength, 1:1] = reshape(glypharray[1:textlength], textlength)
            else
                textGPU[1:end, 1:remaining] = reshape(glypharray[1:1024*remaining], 1024, remaining)
            end
            obj[:postrender, renderinstanced] = (obj.vertexarray, textlength)
        end
        return (obj, textlength, textGPU, glypharray, selection0, true)
    end
    return (obj, textlength, textGPU, glypharray, selection0, false)
end

# Filters a signal. If any of the items is in the signal, the signal is returned.
# Otherwise default is returned
function filteritems{T}(a::Signal{T}, items, default::T)
  lift(a) do signal
    if any(item-> in(item, signal), items)
      signal
    else
      default
    end
  end
end

function Base.delete!(s::Array{GLGlyph{Uint16}, 1}, Index::Integer)
  if Index == 0
    return s
  elseif Index == length(s)
    return s[1:end-1]
  end
  return [s[1:max(Index-1, 0)], s[min(Index+1,length(s)):end]]
end

addchar(s::Array{GLGlyph{Uint16}, 1}, glyph::Char, Index::Integer) = addchar(s, GLGlyph(glyph, 0, 0, 0), int(Index))
function addchar(s::Array{GLGlyph{Uint16}, 1}, glyph::GLGlyph{Uint16}, i::Integer)
  if i == 0
    return [glyph, s]
  elseif i == length(s)
    return [s, glyph]
  elseif i > length(s) || i < 0
    return s
  end
  return [s[1:i], glyph, s[i+1:end]]
end


text:
module VectorValuedNonLinearFEOperatorsTests

using Test
using Gridap
import Gridap: ∇

# Material parameters
const λ = 30.0
const μ = 40.0

# Constitutive law (St. Venant–Kirchhoff Material)
S(E) = λ*tr(E)*one(E) + 2*μ*E

# Green strain
E(F) = 0.5*( F'*F - one(F) )
dE(F,dF) = 0.5*( dF'*F + F'*dF)

# Operations to be performed at integration points
@law σ(x,∇u) = ∇u * S(E(∇u))
@law dσ(x,∇du,∇u) = ∇du * S(E(∇u)) + ∇u * S(dE(∇u,∇du))

# Define manufactured functions
ufun(x) = VectorValue(x[1] + x[2],x[1])
ufun_grad(x) = TensorValue(1.0,1.0,1.0,0.0)
∇(::typeof(ufun)) = ufun_grad
bfun(x) = VectorValue(0.0,0.0)

# Construct the discrete model
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

# Construct the FEspace
order = 2
diritag = "boundary"
T = VectorValue{2,Float64}
fespace = H1ConformingFESpace(T,model,order,diritag)

# Define test and trial
V = TestFESpace(fespace)
U = TrialFESpace(fespace,ufun)

# Setup integration on the volume
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

# Transform shape function according to the model

# Terms in the volume
bfield = CellField(trian,bfun)
res(u,v) = inner( ∇(v), σ(∇(u)) ) - inner(v,bfield)
jac(u,v,du) = inner(∇(v), dσ(∇(du),∇(u)) )
t_Ω = NonLinearFETerm(res,jac,trian,quad)

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define the FEOperator
op = NonLinearFEOperator(V,U,assem,t_Ω)

# Define the FESolver
ls = LUSolver()
tol = 1.e-10
maxiters = 20
nls = NewtonRaphsonSolver(ls,tol,maxiters)
solver = NonLinearFESolver(nls)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
u = CellField(trian,ufun)
e = u - uh

# Define norms to measure the error
l2(u) = inner(u,u)
sh1(u) = inner(∇(u),∇(u))
h1(u) = sh1(u) + l2(u)

# Compute errors
el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

@test el2 < 1.e-8
@test eh1 < 1.e-8

# Further tests
op = NonLinearFEOperator(V,U,t_Ω)
op = NonLinearFEOperator(V,U,t_Ω,t_Ω)

end # module


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

half_pi{T}(::Type{T}) = pi_interval(T) / 2
half_pi{T<:AbstractFloat}(x::T) = half_pi(T)

two_pi{T}(::Type{T})  = pi_interval(T) * 2

range_atan2{T<:Real}(::Type{T}) = Interval(-(pi_interval(T).hi), pi_interval(T).hi)
half_range_atan2{T}(::Type{T}) = (temp = half_pi(T); Interval(-(temp.hi), temp.hi) )
pos_range_atan2{T<:Real}(::Type{T}) = Interval(zero(T), pi_interval(T).hi)

doc"""Finds the quadrant(s) corresponding to a given floating-point
number. The quadrants are labelled as 0 for x ∈ [0, π/2], etc.
For numbers very near a boundary of the quadrant, a tuple of two quadrants
is returned. The minimum or maximum must then be chosen appropriately.

This is a rather indirect way to determine if π/2 and 3π/2 are contained
in the interval; cf. the formula for sine of an interval in
Tucker, *Validated Numerics*."""

function find_quadrants(x::AbstractFloat)
    temp = x / half_pi(x)
    (floor(Int, temp.lo), floor(Int, temp.hi))
end

function sin{T}(a::Interval{T})
    isempty(a) && return a

    whole_range = Interval(-one(T), one(T))

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return Interval(-one(T), one(T))
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range  # in same quadrant but separated by almost 2pi
        lo = Interval(sin(a.lo, RoundDown), sin(a.lo, RoundUp))
        hi = Interval(sin(a.hi, RoundDown), sin(a.hi, RoundUp))
        return hull(lo, hi)

    elseif lo_quadrant==3 && hi_quadrant==0
        return Interval(sin(a.lo, RoundDown), sin(a.hi, RoundUp))

    elseif lo_quadrant==1 && hi_quadrant==2
        return Interval(sin(a.hi, RoundDown), sin(a.lo, RoundUp))

    elseif ( lo_quadrant == 0 || lo_quadrant==3 ) && ( hi_quadrant==1 || hi_quadrant==2 )
        return Interval(min(sin(a.lo, RoundDown), sin(a.hi, RoundDown)), one(T))

    elseif ( lo_quadrant == 1 || lo_quadrant==2 ) && ( hi_quadrant==3 || hi_quadrant==0 )
        return Interval(-one(T), max(sin(a.lo, RoundUp), sin(a.hi, RoundUp)))

    else#if( lo_quadrant == 0 && hi_quadrant==3 ) || ( lo_quadrant == 2 && hi_quadrant==1 )
        return whole_range
    end
end

function cos{T}(a::Interval{T})
    isempty(a) && return a

    whole_range = Interval(-one(T), one(T))

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return Interval(-one(T), one(T))
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant # Interval limits in the same quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range
        lo = Interval(cos(a.lo, RoundDown), cos(a.lo, RoundUp))
        hi = Interval(cos(a.hi, RoundDown), cos(a.hi, RoundUp))
        return hull(lo, hi)

    elseif lo_quadrant == 2 && hi_quadrant==3
        return Interval(cos(a.lo, RoundDown), cos(a.hi, RoundUp))

    elseif lo_quadrant == 0 && hi_quadrant==1
        return Interval(cos(a.hi, RoundDown), cos(a.lo, RoundUp))

    elseif ( lo_quadrant == 2 || lo_quadrant==3 ) && ( hi_quadrant==0 || hi_quadrant==1 )
        return Interval(min(cos(a.lo, RoundDown), cos(a.hi, RoundDown)), one(T))

    elseif ( lo_quadrant == 0 || lo_quadrant==1 ) && ( hi_quadrant==2 || hi_quadrant==3 )
        return Interval(-one(T), max(cos(a.lo, RoundUp), cos(a.hi, RoundUp)))

    else#if ( lo_quadrant == 3 && hi_quadrant==2 ) || ( lo_quadrant == 1 && hi_quadrant==0 )
        return whole_range
    end
end

function tan{T}(a::Interval{T})
    isempty(a) && return a

    diam(a) > pi_interval(T).lo && return entireinterval(a)

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    lo_quadrant_mod = mod(lo_quadrant, 2)
    hi_quadrant_mod = mod(hi_quadrant, 2)

    if lo_quadrant_mod == 0 && hi_quadrant_mod == 1
        (half_pi(T) ⊆ a || -half_pi(T) ⊆ a) && return entireinterval(a)
    elseif lo_quadrant_mod == hi_quadrant_mod && hi_quadrant > lo_quadrant
        hi_quadrant == lo_quadrant+2 && return entireinterval(a)
    end

    Interval(tan(a.lo, RoundDown), tan(a.hi, RoundUp))
end

function asin{T}(a::Interval{T})

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    Interval(asin(a.lo, RoundDown), asin(a.hi, RoundUp))
end

function acos{T}(a::Interval{T})

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    Interval(acos(a.hi, RoundDown), acos(a.lo, RoundUp))
end

function atan{T}(a::Interval{T})
    isempty(a) && return a

    Interval(atan(a.lo, RoundDown), atan(a.hi, RoundUp))
end

#atan2{T<:Real, S<:Real}(y::Interval{T}, x::Interval{S}) = atan2(promote(y, x)...)

function atan2(y::Interval{Float64}, x::Interval{Float64})
    (isempty(y) || isempty(x)) && return emptyinterval(Float64)

    float(atan2(big53(y), big53(x)))
end

function atan2(y::Interval{BigFloat}, x::Interval{BigFloat})
    (isempty(y) || isempty(x)) && return emptyinterval(BigFloat)

    T = BigFloat

    # Prevents some non-sense results whenever y has a zero signed
    if y.lo == zero(T)
        y = Interval(zero(T), y.hi)
    end
    if y.hi == zero(T)
        y = Interval(y.lo, zero(T))
    end

    # Check cases based on x
    if x == zero(x)

        y == zero(y) && return emptyinterval(T)
        y.lo ≥ zero(T) && return half_pi(T)
        y.hi ≤ zero(T) && return -half_pi(T)
        return half_range_atan2(T)

    elseif x.lo > zero(T)

        y == zero(y) && return y
        y.lo ≥ zero(T) &&
            return @round(T, atan2(y.lo, x.hi), atan2(y.hi, x.lo)) # refinement lo bound
        y.hi ≤ zero(T) &&
            return @round(T, atan2(y.lo, x.lo), atan2(y.hi, x.hi))
        return @round(T, atan2(y.lo, x.lo), atan2(y.hi, x.lo))

    elseif x.hi < zero(T)

        y == zero(y) && return pi_interval(T)
        y.lo ≥ zero(T) &&
            return @round(T, atan2(y.hi, x.hi), atan2(y.lo, x.lo))
        y.hi < zero(T) &&
            return @round(T, atan2(y.hi, x.lo), atan2(y.lo, x.hi))
        return range_atan2(T)

    else # zero(T) ∈ x

        if x.lo == zero(T)
            y == zero(y) && return y
            y.lo ≥ zero(T) &&
                return @round(T, atan2(y.lo, x.hi), half_range_atan2(T).hi)
            y.hi ≤ zero(T) &&
                return @round(T, half_range_atan2(T).lo, atan2(y.hi, x.hi))
            return half_range_atan2(T)
        elseif x.hi == zero(T)
            y == zero(y) && return pi_interval(T)
            y.lo ≥ zero(T) &&
                return @round(T, half_pi(T).lo, atan2(y.lo, x.lo))
            y.hi < zero(T) &&
                return @round(T, atan2(y.hi, x.lo), -(half_pi(T).lo))
            return range_atan2(T)
        else
            y.lo ≥ zero(T) &&
                return @round(T, atan2(y.lo, x.hi), atan2(y.lo, x.lo))
            y.hi < zero(T) &&
                return @round(T, atan2(y.hi, x.lo), atan2(y.hi, x.hi))
            return range_atan2(T)
        end

    end
end


text:
include("static_mala_hmc.jl")
include("static_map_optimize.jl")
include("static_mh.jl")
include("static_mh_collapsed.jl")
include("quickstart.jl")
include("dynamic_map_optimize_gibbs.jl")
include("dynamic_mh.jl")
include("dynamic_resimulation_mh.jl")


text:
function gemm_kernel!(C, A, B)
    @avx for n ∈ 1:size(A, 1), m ∈ 1:size(B, 2)
        Cₙₘ = zero(eltype(C))
        for k ∈ 1:size(A, 2)
            Cₙₘ += A[n,k] * B[k,m]
        end
        C[n,m] = Cₙₘ
    end
end

function add_gemm_kernel!(C, A, B, ::Val{factor}) where {factor}
    if factor == 1
        _add_gemm_kernel!(C, A, B)
    elseif factor == -1
        _sub_gemm_kernel!(C, A, B)
    else
        _add_gemm_kernel!(C, C, B, Val(factor))
    end
end

function _add_gemm_kernel!(C, A, B)
    @avx for n ∈ 1:size(A, 1), m ∈ 1:size(B, 2)
        Cₙₘ = zero(eltype(C))
        for k ∈ 1:size(A, 2)
            Cₙₘ += A[n,k] * B[k,m]
        end
        C[n,m] += Cₙₘ
    end
end

function _sub_gemm_kernel!(C, A, B)
    @avx for n ∈ 1:size(A, 1), m ∈ 1:size(B, 2)
        Cₙₘ = zero(eltype(C))
        for k ∈ 1:size(A, 2)
            Cₙₘ -= A[n,k] * B[k,m]
        end
        C[n,m] += Cₙₘ
    end
end

_add_gemm_kernel!(C, A, B, ::Val{1}) = _add_gemm_kernel!(C, A, B)
_add_gemm_kernel!(C, A, B, ::Val{-1}) = _sub_gemm_kernel!(C, A, B)


text:
x = 0
c = Clock(tick=10.0) do count, origin, delta
    global x
    x = x + 1
    if count>10 return -1 end
end

sleep(1)

@test x == 11


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

export
    SmilesParser,
    SmartsParser,
    smilestomol,
    smartstomol

mutable struct SmartsParserState{N<:AbstractNode,E<:UndirectedEdge}
    input::String
    allow_disconnected::Bool

    pos::Int
    done::Bool
    node::Int # No. of current node
    branch::Int # No. of node at the current branch root
    root::Int # No. of node at the current tree root
    ringlabel::Dict{Int,Tuple{Int,Union{E,Symbol,Nothing}}} # TODO: strange union

    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{N}
    edgeattrs::Vector{E}
    connectivity::Vector{Vector{Int}}

    function SmartsParserState{N,E}(str, disconn
            ) where {N<:AbstractNode,E<:UndirectedEdge}
        new(str, disconn, 1, false, 0, 1, 1, Dict(), [], [], [], [])
    end
end

SmilesParser = SmartsParserState{SmilesAtom,SmilesBond}
SmartsParser = SmartsParserState{SmartsAtom,SmartsBond}

function Base.parse(::Type{SMILES}, str::AbstractString)
    state = SmilesParser(str, true)
    fragment!(state)
    return graphmol(state.edges, state.nodeattrs, state.edgeattrs)
end

function Base.parse(::Type{SMARTS}, str::AbstractString)
    if occursin('.', str)
        state = SmartsParser(str, true)
        componentquery!(state)
    else
        state = SmartsParser(str, false)
        fragment!(state)
    end
    return querymol(
        state.edges, state.nodeattrs, state.edgeattrs, state.connectivity)
end

"""
    smilestomol(smiles::AbstractString) -> GraphMol{SmilesAtom,SmilesBond}

Parse SMILES string into `GraphMol` object.
"""
function smilestomol(smiles::AbstractString)
    mol = parse(SMILES, smiles)
    kekulize!(mol)
    setdiastereo!(mol)
    return mol
end

"""
    smartstomol(smarts::AbstractString) -> QueryMol{SmartsAtom,SmartsBond}

Parse SMARTS string into `QueryMol` object.
"""
smartstomol(smarts::AbstractString) = parse(SMARTS, smarts)

function lookahead(state::SmartsParserState, pos::Int)
    # Negative pos can be used
    newpos = state.pos + pos
    if newpos > length(state.input) || newpos < 1
        return Char(0)
    else
        c = state.input[newpos]
        if isascii(c)
            return state.input[newpos]
        else
            throw(ErrorException("invalid charactor $(c)"))
        end
    end
end

Base.read(state) = lookahead(state, 0)  # TODO: rename

function forward!(state::SmartsParserState, num::Int)
    # Negative num can be used
    state.pos += num
    if state.pos > length(state.input)
        state.done = true
    elseif state.done
        throw(ErrorException("charactors in the buffer were used up"))
    elseif state.pos < 1
        throw(ErrorException("no more backtracking!"))
    end
end

forward!(state) = forward!(state, 1)
backtrack!(state, num) = forward!(state, -num)
backtrack!(state) = backtrack!(state, 1)


text:
using GPUifyLoops

f1(x) = sin(x)
f(x) = 1 + f1(x)

function kernel!(A, B, h)
    @inbounds @loop for i in (1:size(A,1); threadIdx().x)
        A[i] = h(B[i])
    end
    nothing
end

data = rand(Float32, 1024)
fdata = similar(data)
kernel!(fdata, data, f)

@assert f.(data) ≈ fdata

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    function kernel!(A::CuArray, B::CuArray)
        @launch CUDA() threads=length(A) kernel!(A, B, f)
    end

    cudata = CuArray(data)
    cufdata = similar(cudata)
    kernel!(cufdata, cudata)

    @assert f.(data) ≈ cufdata
end