text:
######################################################################
# binned_diromni.jl
#
# like diromni but discretized
######################################################################

export BinnedDirOmni
struct BinnedDirOmni <: Sensor
    means::Vector{Float64}
    stds::Vector{Float64}

    bin_range::UnitRange{Int64}

    function BinnedDirOmni(file::AbstractString)
        means = vec(readcsv(file)[:,2])
        stds = 2*ones(360)
        bin_range = -30:30
        return new(means, stds, bin_range)
    end
end

function observe(m::SearchDomain, s::BinnedDirOmni, p::Pose)
    # determine the relative bearing
    rel_bearing = p[3] - true_bearing(p, m.theta)
    if rel_bearing < 0.0
        rel_bearing += 360.0
    end
    rel_int = round(Int, rel_bearing, RoundDown) + 1

    float_obs = s.means[rel_int] + s.stds[rel_int]*randn()
    return obs2bin(float_obs, s)
end

function O(s::BinnedDirOmni, theta::LocTuple, xp::Pose, o)
    #rel_bearing = x.heading - true_bearing(xp, theta)
    rel_bearing = xp[3] - true_bearing(xp, theta)
    #println("rel_bearing = ", rel_bearing)
    if rel_bearing < 0.0
        rel_bearing += 360.0
    end
    # had to do this because rel_bearing was sometimes -2.8e-14
    # when adding this to 360.0, julia returned 360.0
    if rel_bearing == 360.0
        rel_bearing -= 1.0
    end
    rel_int = round(Int, rel_bearing, RoundDown) + 1

    low_val = floor(o)
    high_val = low_val + 1
    d = Normal(s.means[rel_int], s.stds[rel_int])
    p = cdf(d, high_val) - cdf(d, low_val)
    return p
end

# here, num_bins isn't too important; we just bin to nearest integer
obs2bin(o::Float64, s::BinnedDirOmni) = round(Int, o, RoundDown)


text:
module adapters_postgresql_subquery_test

using Test # @test
using Octo.Adapters.PostgreSQL # from to_sql Schema.model SELECT FROM WHERE ORDER BY LIMIT INNER JOIN LATERAL ON TRUE

struct Game
end
Schema.model(Game, table_name="games")

struct GameSold
end
Schema.model(GameSold, table_name="games_sold")

g0 = from(Game, :g0)
gs = from(GameSold, :gs)
f1 = from([SELECT * FROM gs WHERE gs.game_id == g0.id ORDER BY gs.sold_on LIMIT 2], :f1)
@test to_sql([SELECT (g0.name, f1.sold_on) FROM g0 INNER JOIN LATERAL f1 ON TRUE]) ==
             "SELECT g0.name, f1.sold_on FROM games g0 INNER JOIN LATERAL (SELECT * FROM games_sold gs WHERE gs.game_id = g0.id ORDER BY gs.sold_on LIMIT 2) AS f1 ON TRUE"

end # module adapters_postgresql_subquery_test

ext: jl


text:
using Weave.Dates

test_doctypes = filter(first.(Weave.list_out_formats())) do doctype
    # don't test doctypes which need external programs
    doctype ∉ ("pandoc2html", "pandoc2pdf", "md2pdf")
end

function test_func(body)
    @test !isempty(body)
    date_str = string(Date(now()))
    @test occursin(date_str, body)
end

# julia markdown
julia_markdown_body = """
# doc chunk

this is text with `j :inline` code

code chunk:
```julia
using Dates
Date(now())
```
"""

for doctype in test_doctypes
    test_mock_weave(test_func, julia_markdown_body; informat = "markdown", doctype = doctype)
end

# TODO: test noweb format

# julia script
julia_script_body = """
#' # doc chunk
#'
#' this is text with `j :inline` code
#'
#' code chunk:
#+

using Dates
Date(now())
"""
for doctype in test_doctypes
    test_mock_weave(test_func, julia_script_body; informat = "script", doctype = doctype)
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03kd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 6
    NMAX = 50
    LDA1 = NMAX
    LDA2 = NMAX
    LDQ1 = NMAX
    LDQ2 = NMAX
    LDWORK = max( 42*KMAX + NMAX, 80*KMAX - 48, KMAX + max( 2*NMAX, 8*KMAX ) )
    LIWORK = max( 4*KMAX, 2*KMAX + NMAX )
    HUND = 1.0e2
    ZERO = 0.0e0
    QIND = Array{BlasInt,1}(undef, KMAX)
    S = Array{BlasInt,1}(undef, KMAX)
    A = Array{Float64,3}(undef, LDA1,LDA2,KMAX)
    Q = Array{Float64,3}(undef, LDQ1,LDQ2,KMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX)
    ALPHAI = Array{Float64,1}(undef, NMAX)
    BETA = Array{Float64,1}(undef, NMAX)
    SCAL = Array{BlasInt,1}(undef, NMAX)
    ND = Array{BlasInt,1}(undef, KMAX)
    NI = Array{BlasInt,1}(undef, KMAX)
    SELECT = Array{Bool,1}(undef, NMAX)
    T = Array{Float64,1}(undef, NMAX*NMAX*KMAX)
    IXT = Array{BlasInt,1}(undef, KMAX)
    QK = Array{Float64,1}(undef, NMAX*NMAX*KMAX)
    IXQ = Array{BlasInt,1}(undef, KMAX)
# WARNING: desperate attempt to initialize M
    M = 0
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    LDQ = Array{BlasInt,1}(undef, KMAX)
    LDT = Array{BlasInt,1}(undef, KMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    DEFL =  vs[2][1]
    COMPQ =  vs[3][1]
    STRONG =  vs[4][1]
    K = parse(BlasInt, vs[5])
    N = parse(BlasInt, vs[6])
    H = parse(BlasInt, vs[7])
    ILO = parse(BlasInt, vs[8])
    IHI = parse(BlasInt, vs[9])
    if ( N<0 || N>NMAX )
        @error "illegal N"
    end
    TOL = HUND

    vs = String[]
    _isz = K
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    S[1:_isz] .= parsex.(BlasInt, vs)

    vs = String[]
    _isz,_jsz,_ksz = (N,N,K)
    while length(vs) < _isz*_jsz*_ksz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for k in 1:_ksz
      for i in 1:_isz
        _i0 = (i-1)*_jsz + (k-1)*_jsz*_isz
        A[i,1:_jsz,k] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
      end
    end

    if ( LSAME( COMPQ, 'U' ) )
        vs = String[]
        _isz,_jsz,_ksz = (N,N,K)
        while length(vs) < _isz*_jsz*_ksz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        for k in 1:_ksz
            for i in 1:_isz
                _i0 = (i-1)*_jsz + (k-1)*_jsz*_isz
                Q[i,1:_jsz,k] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
            end
        end
    end

    if ( LSAME( COMPQ, 'P' ) )

        vs = String[]
        _isz = K
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        QIND[1:_isz] .= parsex.(BlasInt, vs)
        vs = String[]
        _isz,_jsz,_ksz = (N,N,K)
        while length(vs) < _isz*_jsz*_ksz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        for k in 1:_ksz
            if QIND[k] > 0
                for i in 1:_isz
                    _i0 = (i-1)*_jsz + (k-1)*_jsz*_isz
                    Q[i,1:_jsz,QIND[k]] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
                end
            end
        end
    end # if
    close(f)
    if LSAME(JOB,'E')
        JOB = 'S'
    end
    println(io, "JOB = $JOB")
    println(io, "COMPQ = $COMPQ")
# interp call 1

    INFO, IWARN = SLICOT.mb03bd!(JOB, DEFL, COMPQ, QIND, K, N, H, ILO, IHI, S, A, Q, ALPHAR, ALPHAI, BETA, SCAL, LIWORK, LDWORK)
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")

    if ( IWARN!=0 )
        @warn "mb03bd returned IWARN=$IWARN"
    end
    verbose = false
    if verbose
        ## added for diagnostic
        if ( LSAME( JOB, 'S' ) || LSAME( JOB, 'T' ) )
            # interp output 1
            println(io, "after Schur (HessIdx=$H) A:")
            _nc = N
            _nr = N
            _nk = K
            show(io, "text/plain", A[1:_nr,1:_nc,1:_nk])
            println(io)
        end # if
        if ( LSAME( COMPQ, 'U' ) || LSAME( COMPQ, 'I' ) )
            println(io, "Q:")
            _nc = N
            _nr = N
            _nk = K
            show(io, "text/plain", Q[1:_nr,1:_nc,1:_nk])
            println(io)
        elseif ( LSAME( COMPQ, 'P' ) )
            println(io, "TODO: show Q for COMPQ='P'")
            for L in 1:K
            #if ( QIND[L]>0 )
                #println(io, "Q:")
                #show(io, "text/plain", Q)
                #println(io)
                #nd # if
            end # for
        end # if
    end ## diagnostic

    # prepare data for calling MB03KD with screwy data structures and factor ordering
    for L in 1:K
        ND[L] = max(1,N)
        NI[L] = 0
        LDT[L] = max(1,N)
        IXT[L] = (L-1)*LDT[L]*N + 1
        LDQ[L] = max(1,N)
        IXQ[L] = IXT[L]
        if ( L<= K÷2 )
            i = S[ K - L + 1 ]
            S[K-L+1] = S[L]
            S[L] = i
        end # if
    end
    for L in 1:K
        # interp call 2
        #    FOREIGN.dlacpy!( 'Full', N, N, A( 1, 1, K-L+1 ), LDA1, T( IXT( L ) ), LDT( L ) )
        ir1 = IXT[L]
        for ic in 1:N
            T[ir1:ir1+N-1] .= A[1:N,ic,K-L+1]
            ir1 += LDT[L]
        end
    end

    if verbose
        println(io, "stored T:")
        for k in 1:K
            println(io, "K=$k")
            Ttmp = reshape(T[IXT[k]:IXT[k]+LDT[k]*N-1],N,N)
            show(io, "text/plain", Ttmp)
            println(io)
        end
    end # diagnostic

    if ( LSAME( COMPQ, 'U' ) || LSAME( COMPQ, 'I' ) )
        COMPQ = 'U'

        for L in 1:K
            # interp call 3
            # FOREIGN.dlacpy!( 'Full', N, N, Q( 1, 1, K-L+1 ), LDQ1, QK( IXQ( L ) ), LDQ( L ) )
            ir1 = IXQ[L]
            for ic in 1:N
                QK[ir1:ir1+N-1] .= Q[1:N,ic,K-L+1]
                ir1 += LDQ[L]
            end
        end
    elseif ( LSAME( COMPQ, 'P' ) )
        COMPQ = 'W'
        for L in 1:K
            if QIND[L] < 0
                QIND[L] = 2
            end
            P = QIND[L]
            if P != 0
                ir1 = IXQ[P]
                Qk[ir1:ir1+N-1] .= Q[1:N,ic,K-P+1]
                ir1 += LDQ[P]
            end
        end
    end # if

    if verbose    # diag
        println(io, "stored Q:")
        for k in 1:K
            println(io, "K=$k")
            Ttmp = reshape(QK[IXQ[k]:IXQ[k]+LDQ[k]*N-1],N,N)
            show(io, "text/plain", Ttmp)
            println(io)
        end
    end
    SELECT .= false
    for i in 1:N
        SELECT[i] = ALPHAR[i] < 0
    end
        # interp output 1
    println(io, "ALPHAR:")
    _nr = N
    show(io, "text/plain", ALPHAR[1:_nr])
    println(io)

    # interp output 2
    println(io, "ALPHAI:")
    _nr = N
    show(io, "text/plain", ALPHAI[1:_nr])
    println(io)

    # interp output 3
    println(io, "BETA:")
    _nr = N
    show(io, "text/plain", BETA[1:_nr])
    println(io)

    # interp output 4
    println(io, "SCAL:")
    _nr = N
    show(io, "text/plain", SCAL[1:_nr])
    println(io)

    # interp call 4
    select = convert(Vector{BlasBool}, SELECT)
    println(io, "select:"); show(io, "text/plain", select[1:K]); println()
    M, INFO = SLICOT.mb03kd!(COMPQ, QIND, STRONG, K, N, H, ND, NI, S, select, T, LDT, IXT, QK, LDQ, IXQ, TOL)
    @test INFO == 0
    if INFO != 0
        println(io, "mb03kd returned INFO = $INFO")
        return
    end

    if ( INFO!=0 )
        #@error "mb03kd returned INFO=$INFO"
    end
    Tf = zeros(N,N,K)
    for L in 1:K
        P = K - L + 1
        ir1 = IXT[P]
        for ic in 1:N
            Tf[1:N,ic,L] .= T[ir1:ir1+N-1]
            ir1 += LDT[P]
        end
    end
    # unable to translate write loop:
    #  write ( IXT( P ) + I - 1 + ( J - 1 )*LDT( P ) ), J = 1, N
    # interp output 5
    println(io, "T:")
    show(io, "text/plain", Tf)
    println(io)

    if ( LSAME( COMPQ, 'U' ) || LSAME( COMPQ, 'I' ) )
        Qf = zeros(N,N,K)
        for L in 1:K
            P = K - L + 1
            ir1 = IXQ[P]
            for ic in 1:N
                Qf[1:N,ic,L] .= QK[ir1:ir1+N-1]
                ir1 += LDQ[P]
            end
        end
        # unable to translate write loop:
        #  write ( IXQ( P ) + I - 1 + ( J - 1 )*LDQ( P ) ), J = 1, N
        # interp output 6
        println(io, "QK:")
        show(io, "text/plain", Qf)
        println(io)

    elseif ( LSAME( COMPQ, 'W' ) )
        for L in 1:K
            if ( QIND[L]>0 )
                println(io, "factor $(QIND[L]):")
                Qf = zeros(N,N)
                P = K - QIND[L] + 1
                ir1 = IXQ[P]
                for ic in 1:N
                    Qf[1:N,ic] .= QK[ir1:ir1+N-1]
                    ir1 += LDQ[P]
                end
                show(io, "text/plain", Qf)
                println(io)
            end
        end
            # unable to translate write loop:
            #  write ( IXQ( P ) + I - 1 + ( J - 1 )*LDQ( P ) ), J = 1, N
            # interp output 7
    end # if
end # module


text:
function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/hover")}}, params)
    return TextDocumentPositionParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/hover")},TextDocumentPositionParams}, server)
    if !haskey(server.documents, URI2(r.params.textDocument.uri))
        send(JSONRPC.Response(r.id, CancelParams(r.id)), server)
        return
    end
    documentation = Any[]
    doc = server.documents[URI2(r.params.textDocument.uri)]
    x = get_expr(getcst(doc), get_offset(doc, r.params.position))

    get_hover(x, documentation, server)

    send(JSONRPC.Response(r.id, Hover(unique(documentation))), server)
end

function get_hover(x, documentation, server) end

function get_hover(x::EXPR, documentation, server)
    if x.parent isa EXPR  && (x.kind === CSTParser.Tokens.END || x.kind === CSTParser.Tokens.RPAREN || x.kind === CSTParser.Tokens.RBRACE || x.kind === CSTParser.Tokens.RSQUARE)
        push!(documentation, MarkedString("Closes $(x.parent.typ) expression."))
    elseif CSTParser.isidentifier(x) && StaticLint.hasref(x)
        if x.ref isa StaticLint.Binding
            get_hover(x.ref, documentation, server)
        elseif x.ref isa SymbolServer.SymStore
            append!(documentation, split_docs(x.ref.doc))
        end
    end
end

function get_hover(b::CSTParser.Binding, documentation, server)
    if b.val isa EXPR
        if CSTParser.defines_function(b.val)
            while true
                if b.val isa EXPR
                    if CSTParser.defines_function(b.val)
                        pushfirst!(documentation, MarkedString(Expr(CSTParser.get_sig(b.val))))
                    elseif CSTParser.defines_datatype(b.val)
                        pushfirst!(documentation, MarkedString(Expr(b.val)))
                    end
                elseif b.val isa SymbolServer.SymStore
                    push!(documentation, b.val.doc)
                else
                    break
                end
                if b.overwrites isa CSTParser.Binding && b.overwrites != b && (b.overwrites.t == getsymbolserver(server)["Core"].vals["Function"] || b.overwrites.t == getsymbolserver(server)["Core"].vals["DataType"])
                    b = b.overwrites
                else
                    break
                end
            end
        else
            push!(documentation, MarkedString(Expr(b.val)))
        end
    elseif b.val isa SymbolServer.SymStore
        append!(documentation, split_docs(b.val.doc))
    elseif b.val isa CSTParser.Binding
        get_hover(b.val, documentation, server)
    end
end

"""
    split_docs(s::String)

Returns an array of Union{String,MarkedString} by separating code blocks (denoted by ```sometext```) within s.
"""
function split_docs(s::String)
    out = Any[]
    locs = Int[]
    i = 1
    while i < length(s)
        m = match(r"```", s, i)
        if m isa Nothing
            isempty(out) && push!(out, s)
            break
        else
            push!(locs, m.offset)
            i = m.offset + 3
        end
    end
    if length(locs) > 0 && iseven(length(locs))
        if locs[1] !== 1
            push!(out, s[1:locs[1]-1])
        end
        for i = 1:2:length(locs)
            push!(out, LanguageServer.MarkedString("julia", replace(s[locs[i]+3:locs[i+1]-1], "jldoctest"=>"")))
            if i + 1 == length(locs)
                push!(out, s[locs[i+1]+3:end])
            else
                push!(out, s[locs[i+1]+3:locs[i+2]-1])
            end
        end
    end
    return out
end


text:
using PDMats
import Statistics: mean, var
using LinearAlgebra

"""
   `calc_periodogram_orig(times, y_obs, σ_ys; options...)`
Computes a periodogram for times, observations and uncertainties.

Inputs:
- times: obseration times
- y_obs: observed values at each time
- σ_ys:  measurement uncertainty at each time (may be either a Vector or positive definite matrix)
Optional parameters:

- period_min (2.0)
- period_max (4*span of times)
- num_periods (4000)

Output (as NamedTuple):
- period_grid: orbital periods searched
- periodogram: χ² for best fit given each orbital period
- period_best_fit: period searched resulting to minimum χ²
- coeff_best_fit: best-fit coefficients for period_best_fit
- phase_best_fit: orbital phase for period_best_fit at time = 0
- predict: predictions of best-fit model at each time
- rms: root mean square of predictions (of best-fit model) minus observations
"""
function calc_periodogram_orig  end

function calc_periodogram_orig(t::Vector, y_obs::Vector, σ_ys::Vector; kwargs...)
 	calc_periodogram_orig(t, y_obs, PDiagMat(σ_ys.^2); kwargs...)
end

function calc_periodogram_orig(t::Vector, y_obs::Vector, covar_mat::MT;
	    		period_min::Real = 2.0, period_max::Real = 4*(maximum(t)-minimum(t)), num_periods::Integer = 4000) where {
						T3<:Real, MT<:AbstractPDMat{T3} }
	period_grid =  1.0 ./ range(1.0/period_max, stop=1.0/period_min, length=num_periods)
	periodogram = map(p->-0.5*calc_χ²_general_linear_least_squares_orig(calc_design_matrix_circ(p,t),covar_mat,y_obs),period_grid)
	period_fit = period_grid[argmax(periodogram)]
	design_matrix_fit = calc_design_matrix_circ(period_fit,t)
	coeff_fit = fit_general_linear_least_squares_orig(design_matrix_fit,covar_mat,y_obs)
	phase_fit = atan(coeff_fit[1],coeff_fit[2])
	pred = design_matrix_fit * coeff_fit
	rms = sqrt(mean((y_obs.-pred).^2))
	return (;period_grid=period_grid, periodogram=periodogram, period_best_fit = period_fit, coeff_best_fit=coeff_fit, phase_best_fit=phase_fit, predict=pred, rms=rms )
end

"""
   `calc_periodogram(times, y_obs, σ_ys; options...)`
Computes a periodogram for times, observations and uncertainties.

Inputs:
- times: obseration times
- y_obs: observed values at each time
- σ_ys:  measurement uncertainty at each time (may be either a Vector or positive definite matrix)
Optional parameters:

- period_min (2.0)
- period_max (4*span of times)
- num_periods (4000)

Output (as NamedTuple):
- period_grid: orbital periods searched
- periodogram: χ² for best fit given each orbital period
- period_best_fit: period searched resulting to minimum χ²
- coeff_best_fit: best-fit coefficients for period_best_fit
- phase_best_fit: orbital phase for period_best_fit at time = 0
- predict: predictions of best-fit model at each time
- rms: root mean square of predictions (of best-fit model) minus observations
"""
function calc_periodogram end

function calc_periodogram(t::Vector, y_obs::Vector, σ_ys::Vector; kwargs...)
 	calc_periodogram(t, y_obs, PDiagMat(σ_ys.^2), kwargs...)
end

function calc_periodogram(t::Vector, y_obs::Vector, covar_mat::MT;
				period_min::Real = 2.0, period_max::Real = 4*(maximum(t)-minimum(t)), num_periods::Integer = 4000) where {
						MT<:AbstractPDMat }
	period_grid =  1.0 ./ range(1.0/period_max, stop=1.0/period_min, length=num_periods)
	design_matrix_workspace = Array{eltype(t)}(undef,length(t),3)
	calc_χ²_workspace = prealloc_workspace(design_matrix_workspace,covar_mat,y_obs)
	periodogram = map(p->-0.5*calc_χ²_general_linear_least_squares(calc_design_matrix_circ!(design_matrix_workspace,p,t),covar_mat,y_obs,workspace=calc_χ²_workspace),period_grid)
	period_fit = period_grid[argmax(periodogram)]
	design_matrix_fit = calc_design_matrix_circ!(design_matrix_workspace,period_fit,t)
	coeff_fit = fit_general_linear_least_squares(design_matrix_fit,covar_mat,y_obs)
	phase_fit = atan(coeff_fit[1],coeff_fit[2])
	pred = design_matrix_fit * coeff_fit
	rms = sqrt(mean((y_obs.-pred).^2))
	return (;period_grid=period_grid, periodogram=periodogram, period_best_fit = period_fit, coeff_best_fit=coeff_fit, phase_best_fit=phase_fit, predict=pred, rms=rms )
end

# Internal routines called by calc_periodogram[_orig]
function calc_design_matrix_circ!(result::AM, period::Real, times::AV) where { R1<:Real, AM<:AbstractMatrix{R1}, AV<:AbstractVector{R1} }
	n = length(times)
	@assert size(result,1) == n
	@assert 2<= size(result,2) <=3
	for i in 1:n
		( result[i,1], result[i,2] ) = sincos(2π/period .* times[i])
	end
	if size(result,2) == 3
		result[:,3] .= one(eltype(result))
	end
	return result
end

function calc_design_matrix_circ(period::Real, times::AV) where { R1<:Real, AV<:AbstractVector{R1} }
	n = length(times)
	dm = Array{promote_type(typeof(period),eltype(times))}(undef,n,3)
	calc_design_matrix_circ!(dm,period,times)
	return dm
end

function prealloc_workspace( design_mat::ADM, covar_mat::APD, obs::AA) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	(;	invA_X = Array{Float64,2}(undef,size(covar_mat,1),size(design_mat,2)),
		Xt_invA_X =  Array{Float64,2}(undef,size(design_mat,2),size(design_mat,2)),
		invA_y =  Array{Float64,1}(undef,length(obs)),
		X_invA_y =  Array{Float64,1}(undef,size(design_mat,2)),
		AB_hat = Array{Float64,1}(undef,size(design_mat,2) ),
		predict =  Array{Float64,1}(undef,length(obs))  )
end

function fit_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA;
			workspace = prealloc_workspace(design_mat,covar_mat,obs)
	 			) where { ADM<:AbstractMatrix, APD<:PDiagMat, AA<:AbstractArray }
	#workspace.invA_X .= covar_mat \ design_mat
	workspace.invA_X .= covar_mat.inv_diag .* design_mat
	mul!(workspace.Xt_invA_X,design_mat', workspace.invA_X)
	#workspace.inv A_y .= covar_mat \ obs
	workspace.invA_y .= covar_mat.inv_diag .* obs
	mul!(workspace.X_invA_y,  design_mat', workspace.invA_y)
	#workspace.AB_hat .= workspace.Xt_invA_X \ workspace.X_invA_y
	lufac = lu!( workspace.Xt_invA_X, check=false)
	workspace.AB_hat .= lufac \ workspace.X_invA_y
end

function fit_general_linear_least_squares_orig( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	Xt_inv_covar_X = Xt_invA_X(covar_mat,design_mat)
	X_inv_covar_y =  design_mat' * (covar_mat \ obs)
	AB_hat = Xt_inv_covar_X \ X_inv_covar_y                            # standard GLS
end

function predict_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA;
			workspace = prealloc_workspace(design_mat,covar_mat,obs)
			 	) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	#param =
	fit_general_linear_least_squares(design_mat,covar_mat,obs,workspace=workspace)
	mul!(workspace.predict, design_mat, workspace.AB_hat)
	return workspace.predict
end

function predict_general_linear_least_squares_orig( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	param = fit_general_linear_least_squares_orig(design_mat,covar_mat,obs)
	design_mat * param
end

function calc_χ²_general_linear_least_squares_orig( design_mat::ADM, covar_mat::APD, obs::AA ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	pred = predict_general_linear_least_squares_orig(design_mat,covar_mat,obs)
	invquad(covar_mat,obs-pred)
end

function calc_χ²_general_linear_least_squares( design_mat::ADM, covar_mat::APD, obs::AA;
				workspace = prealloc_workspace(design_mat,covar_mat,obs)
				 ) where { ADM<:AbstractMatrix, APD<:AbstractPDMat, AA<:AbstractArray }
	#workspace.predict .=
	predict_general_linear_least_squares(design_mat,covar_mat,obs,workspace=workspace)
	workspace.predict .-= obs
	invquad(covar_mat,workspace.predict)
end


text:
function pool{T}(x::CuArray{T}; handle=cudnnhandle(), alpha=1,
                 o...) # window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    y = similar(x, pdims(x; o...))
    beta = 0
    @cuda(cudnn, cudnnPoolingForward,
          (Cptr, Cptr, Ptr{T}, Cptr, Ptr{T}, Ptr{T}, Cptr,Ptr{T}),
          handle, PD(x;o...), Ref(T(alpha)), TD(x), x.ptr, Ref(T(beta)), TD(y), y.ptr)
    return y
end

function pool_grad{T}(x::CuArray{T}, y::CuArray{T}, dy::CuArray{T};
                      handle=cudnnhandle(), alpha=1, mode=0,
                      o...) # window=2, padding=0, stride=window, maxpoolingNanOpt=0
    if alpha!=1 && mode==0; error("Gradient of pool(alpha!=1,mode=0) broken in CUDNN"); end
    dx = similar(x)
    beta = 0
    @cuda(cudnn,cudnnPoolingBackward,
          (Cptr, Cptr, Ptr{T}, Cptr, Ptr{T}, Cptr, Ptr{T}, Cptr, Ptr{T}, Ptr{T}, Cptr, Ptr{T}),
          handle, PD(x; mode=mode, o...), Ref(T(alpha)), TD(y), y.ptr,
          TD(dy), dy.ptr, TD(x), x.ptr, Ref(T(beta)), TD(dx), dx.ptr)
    return dx
end

function unpool(x; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,x))
    y = similar(x,updims(x; window=window, o...))
    poolx(y,x,x.*w; o..., window=window, mode=1, alpha=1/alpha)
end

function unpool_grad(dy; window=2, alpha=1, o...) # padding=0, stride=window, mode=0,
                                                    # maxpoolingNanOpt=0
    w = prod(psize(window,dy))
    pool(dy.ptr; o..., window=window, mode=1, alpha=1/alpha) * w
end


text:
@marginalrule typeof(*)(:A_in) (m_out::NormalDistributionsFamily, m_A::PointMass, m_in::F, meta::Any) where { F <: NormalDistributionsFamily } = begin
    # A = mean(m_A)
    # xi_out, W_out = weightedmean_precision(m_out)
    # W = A' * W_out * A
    # b_in = convert(promote_variate_type(F, NormalWeightedMeanPrecision), A' * xi_out, W)
    b_in = @call_rule typeof(*)(:in, Marginalisation) (m_out = m_out, m_A = m_A, meta = meta)
    q_in = prod(ProdAnalytical(), b_in, m_in)
    return (A = m_A, in = q_in)
end


text:
using KernelAbstractions, Test

@kernel function mul2(A)
    I = @index(Global, Cartesian)
    A[I] = 2 * A[I]
end

@kernel function add3(A, B, C)
    I = @index(Global, Cartesian)
    A[I] = B[I] + C[I]
end

@kernel function addi(A, C, i)
    I = @index(Global, Cartesian)
    A[I] = i + C[I]
end

function test_typed_kernel_dynamic(backend, ArrayT)
    A = ArrayT(ones(1024, 1024))
    kernel = mul2(backend())
    res = if backend == CPU
        @ka_code_typed kernel(A, ndrange=size(A), workgroupsize=16)
    else
        @ka_code_typed kernel(A, ndrange=size(A), workgroupsize=(32, 32))
    end
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})
end

function test_typed_kernel_dynamic_no_info(backend, ArrayT)
    A = ArrayT(ones(1024, 1024))
    B = similar(A)
    C = similar(A)
    kernel = add3(backend())
    res = @ka_code_typed kernel(A, B, C, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})
end

function test_typed_kernel_static(backend, ArrayT)
    A = ArrayT(ones(1024, 1024))
    kernel = if backend == CPU
        mul2(backend(), 16)
    else
        mul2(backend(), (32, 32))
    end
    res = @ka_code_typed kernel(A, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})
end

function test_typed_kernel_no_optimize(backend, ArrayT)
    A = ArrayT(ones(1024, 1024))
    B = similar(A)
    C = similar(A)
    kernel = if backend == CPU
        add3(backend(), 16)
    else
        add3(backend(), (32, 32))
    end
    res = @ka_code_typed optimize=false kernel(A, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo,Core.TypeofBottom})
    res_opt = @ka_code_typed kernel(A, ndrange=size(A))
    @test size(res[1].code) < size(res_opt[1].code)
end

function test_expr_kernel(backend, ArrayT)
    A = ArrayT(ones(1024, 1024))
    C = similar(A)
    kernel = if backend == CPU
        addi(backend())
    else
        addi(backend(), (32, 32))
    end
    res = @ka_code_typed kernel(A, C, 1+2, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})
end

function reflection_testsuite(backend, ArrayT)
    test_typed_kernel_dynamic(backend, ArrayT)
    test_typed_kernel_dynamic_no_info(backend, ArrayT)
    test_typed_kernel_static(backend, ArrayT)
    test_typed_kernel_no_optimize(backend, ArrayT)
    test_expr_kernel(backend, ArrayT)
end


text:
module ImpactModeling

import Colors

using DataFrames
using Distributions
using ProgressMeter
using Random: shuffle
using StatsBase
using VegaLite

import JSON

function summarize_vec(vec::Vector{T} where T <: Real)
    res = DataFrame(quantile(vec, [0.0, 0.25, 0.5, 0.75, 1.0])', [:min, :LQ, :median, :UQ, :max])
    res[!, :mean] .= mean(vec)
    return res
end

function summarize_sample(distr::UnivariateDistribution, n_objects::Int, n_samples::Int=10000)
    sample = vcat(mapslices(summarize_vec, rand(distr, (n_samples, n_objects)), dims=2)...)
    res = aggregate(sample, mean)
    names!(res, names(sample))
    return res
end

function best_additional_value(n_samples::Int, n_applicants::Int, distr::UnivariateDistribution)::Vector{Float64}
    samples = mapslices(sort, rand(distr, (n_samples, n_applicants)), dims=2)
    return samples[:, end] - samples[:, end - 1]
end

function estimate_applicant_impact_summary(n_samples::Int, n_applicants::Vector{Int}, distr::UnivariateDistribution)
    best_values = best_additional_value.(n_samples, n_applicants, distr);
    applicant_impact_df = vcat(summarize_vec.(best_values)...)
    applicant_impact_df[!, :n_applicants] = n_applicants;
    return applicant_impact_df
end

function organization_impact_fraction(n_samples::Int, n_organizations::Int, organization_dist::UnivariateDistribution)
    samples = mapslices(x -> sort(x, rev=true), rand(organization_dist, (n_samples, n_organizations)), dims=2)
    frac_of_top = samples ./ samples[:, 1]
    frac_of_total = samples ./ sum(samples, dims=2)
    cum_frac_of_total = mapslices(cumsum, frac_of_total, dims=2)

    res = vcat(mapslices(summarize_vec, frac_of_top, dims=1)...)
    names!(res, [Symbol("top_$s") for s in names(res)])

    res_total = vcat(mapslices(summarize_vec, frac_of_total, dims=1)...)
    names!(res_total, [Symbol("total_$s") for s in names(res_total)])

    res_cum = vcat(mapslices(summarize_vec, cum_frac_of_total, dims=1)...)
    names!(res_cum, [Symbol("cum_total_$s") for s in names(res_cum)])

    res = hcat(res, res_total, res_cum)
    res[!, :n_organizations] .= n_organizations
    res[!, :organization_id] .= 1:size(res, 1)

    return res
end

# total_impact_per_company(n_samples::Int, n_applicants::Vector{Int}, n_organizations::Vector{Int},
#     applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution) =
#         vcat(total_impact_per_company.(n_samples, n_applicants, n_organizations, applicant_dist, organization_dist)...)

function aggregate_impact(applicant_impact::Union{Matrix{T}, Vector{T}, T}, company_impact::Union{Matrix{T}, Vector{T}, T}; method::Symbol=:prod) where T <: Real
    if method == :prod
        return applicant_impact .* company_impact
    elseif method == :geometric
        return sqrt.(applicant_impact .* company_impact)
    else
        error("Unknown method: $method")
    end
end

function total_impact_per_company(personal_impacts::Matrix{Float64}, company_impacts::Matrix{Float64},
        personal_fit_dist::UnivariateDistribution; matching::Symbol=:exact, aggregation::Symbol=:prod)
    n_applicants = size(personal_impacts, 2)
    total_impacts = zeros(size(personal_impacts, 1))
    company_impacts = mapslices(x -> sort(x, rev=true), company_impacts, dims=2)
    for si in 1:size(personal_impacts, 1)
        personal_fits = rand(personal_fit_dist, (n_applicants, size(company_impacts, 2)));
        pers_impacts_adj = personal_fits .* personal_impacts[si,:]
        for ci in 1:size(company_impacts, 2)
            id = 0
            if matching == :exact
                ai = findmax(pers_impacts_adj[:, 1])[2]
            elseif matching == :weighted
                ai = sample(1:n_applicants, FrequencyWeights(pers_impacts_adj[:, 1]), 1)
            elseif matching == :random
                ai = sample(1:n_applicants, 1)
            else
                error("Unknown matching type: $matching")
            end

            total_impacts[si] += aggregate_impact(pers_impacts_adj[ai, 1], company_impacts[si, ci], method=aggregation)
            pers_impacts_adj = pers_impacts_adj[:, 2:end]
            pers_impacts_adj[ai, :] .= 0
        end
    end

    return total_impacts
end

function arrange_personal_impacts(personal_impacts::Matrix{Float64}; matching::Symbol=:exact, n_exact::Int=0)
    if n_exact > 0 && matching != :exact
        personal_impacts = mapslices(x -> sort(x, rev=true), personal_impacts, dims=2)
    end

    if matching == :exact
        personal_impacts = mapslices(x -> sort(x, rev=true), personal_impacts, dims=2)
    elseif matching == :weighted
        personal_impacts[:, (n_exact + 1):end] = mapslices(row -> sample(row, FrequencyWeights(row), length(row); replace=false),
            personal_impacts[:, (n_exact + 1):end], dims=2)
    elseif matching == :random
        personal_impacts[:, (n_exact + 1):end] = mapslices(shuffle, personal_impacts[:, (n_exact + 1):end], dims=2)
    else
        error("Unknown matching type: $matching")
    end

    return personal_impacts
end

function total_impact_per_company(personal_impacts::Matrix{Float64}, company_impacts::Matrix{Float64};
        matching::Symbol=:exact, n_exact::Int=0, aggregation::Symbol=:prod)
    n_min = min(size(personal_impacts, 2), size(company_impacts, 2))
    company_impacts = mapslices(x -> sort(x, rev=true), company_impacts, dims=2)[:, 1:n_min]
    personal_impacts = arrange_personal_impacts(personal_impacts, matching=matching, n_exact=n_exact)[:, 1:n_min]
    return sum(aggregate_impact(personal_impacts, company_impacts, method=aggregation), dims=2) |> vec
end

function total_impact_per_company(n_samples::Int, n_applicants::Int, n_organizations::Int,
        applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution;
        personal_fit_dist::Union{UnivariateDistribution, Nothing}=nothing, matching::Symbol=:exact, aggregation::Symbol=:prod)
    personal_impacts = rand(applicant_dist, (n_samples, n_applicants))
    company_impacts = rand(organization_dist, (n_samples, n_organizations))

    total_impact = personal_fit_dist === nothing ?
        total_impact_per_company(personal_impacts, company_impacts; matching=matching, aggregation=aggregation) :
        total_impact_per_company(personal_impacts, company_impacts, personal_fit_dist; matching=matching, aggregation=aggregation)

    res = summarize_vec(total_impact)
    res[!, :n_applicants] .= n_applicants
    res[!, :n_organizations] .= n_organizations
    return res
end

function impact_loss_from_matching(n_samples::Int, n_applicants::Int, n_organizations::Int,
        applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution;
        matching1::Symbol=:exact, matching2::Symbol=:weighted, n_exact::Int=0, aggregation::Symbol=:prod)

    personal_impacts = rand(applicant_dist, (n_samples, n_applicants))
    company_impacts = rand(organization_dist, (n_samples, n_organizations))

    total_impact = total_impact_per_company(personal_impacts, company_impacts; matching=matching1, n_exact=n_exact, aggregation=aggregation)
    total_impact2 = total_impact_per_company(personal_impacts, company_impacts; matching=matching2, n_exact=n_exact, aggregation=aggregation)

    res = summarize_vec((total_impact .- total_impact2) ./ total_impact)
    res[!, :n_applicants] .= n_applicants
    res[!, :n_organizations] .= n_organizations

    return res
end

function prepare_wrong_choice_info(n_samples::Int, n_applicants::Int, n_organizations::Int,
    applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution)
    personal_impacts = rand(applicant_dist, (n_samples, n_applicants));
    company_impacts = rand(organization_dist, (n_samples, n_organizations));

    n_min = min(n_applicants, n_organizations)
    personal_impacts = mapslices(x -> sort(x, rev=true), personal_impacts, dims=2)[:, 1:n_min]
    company_impacts = mapslices(x -> sort(x, rev=true), company_impacts, dims=2)[:, 1:n_min];

    return personal_impacts, company_impacts
end

function wrong_choice_impact_loss(personal_impacts::Matrix{Float64}, company_impacts::Matrix{Float64},
        shift::Int, real_id::Int; aggregation::Symbol=:prod)
    applied_id = real_id + shift
    personal_impacts = deepcopy(personal_impacts)
    pers_imp_real = deepcopy(personal_impacts[:, real_id])

    original_impacts = pers_imp_real .* company_impacts[:, real_id];
    original_impacts_total = sum(personal_impacts .* company_impacts, dims=2);

    personal_impacts[:, real_id:(applied_id-1)] .= personal_impacts[:, (real_id + 1):applied_id];
    personal_impacts[:, applied_id] .= pers_imp_real

    final_impacts = aggregate_impact(pers_imp_real, company_impacts[:, applied_id], method=aggregation);
    final_impacts_total = sum(aggregate_impact(personal_impacts, company_impacts, method=aggregation), dims=2);

    res = summarize_vec(vec((original_impacts_total .- final_impacts_total) ./ original_impacts))
    rename!(res, [Symbol("personal_$s") for s in names(res)])

    res_total = summarize_vec(vec((original_impacts_total .- final_impacts_total) ./ original_impacts_total))
    rename!(res_total, [Symbol("total_$s") for s in names(res_total)])

    res = hcat(res, res_total)

    res[!, :real_id] .= real_id
    res[!, :shift] .= shift

    return res
end

function wrong_choice_impact_loss(n_samples::Int, n_applicants::Int, n_organizations::Int,
        applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution;
        real_ids::Vector{Int}, shifts::Vector{Int}, aggregation::Symbol=:prod)
    personal_impacts, company_impacts = prepare_wrong_choice_info(n_samples, n_applicants, n_organizations, applicant_dist, organization_dist)

    # @assert (maximum(real_ids) + maximum(shifts)) <= size(personal_impacts, 2)
    return vcat([wrong_choice_impact_loss(personal_impacts, company_impacts, s, ri, aggregation=aggregation)
            for s in shifts for ri in real_ids if s + ri <= size(personal_impacts, 2)]...)
end

function wrong_choice_impact_loss_frac(n_samples::Int, n_applicants::Int, n_organizations::Int,
    applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution;
    real_quants::Vector{Float64}, shift_quants::Vector{Float64}, base::Symbol=:min)
    @assert (maximum(real_quants) <= 1) & (maximum(shift_quants) <= 1) & (minimum(real_quants) >= 0) & (minimum(shift_quants) >= 0)
    base = (base == :applicants) ? n_applicants : ((base == :organizations) ? n_organizations : min(n_applicants, n_organizations))

    real_ids, shift_ids = [unique(round.(Int, (base - 1) .* v .+ 1)) for v in [real_quants, shift_quants]]
    return wrong_choice_impact_loss(n_samples, n_applicants, n_organizations, applicant_dist, organization_dist;
        real_ids=real_ids, shifts=shift_ids)
end

function bin_search_n_apps_per_new_org(n_samples::Int, n_applicants::Int, n_organizations::Int,
        applicant_dist::UnivariateDistribution, organization_dist::UnivariateDistribution;
        max_iters=100, quartile_frac=0.05)
    n_a_upper = n_applicants
    impact_target = total_impact_per_company(n_samples, n_applicants, n_organizations + 1, applicant_dist, organization_dist);

    converged = false
    for i in 1:max_iters
        n_a_upper *= 2
        impact_cur = total_impact_per_company(n_samples, n_a_upper, n_organizations, applicant_dist, organization_dist);
        if impact_cur.median > impact_target.median
            converged = true
            break
        end
    end

    if !converged
        error("First phase didn't converge ($n_applicants, $n_organizations)")
    end

    converged = false
    n_a_mean = n_a_lower = n_applicants + 1
    for i in 1:max_iters
        n_a_mean = round(Int, (n_a_lower + n_a_upper) / 2)
        impact_cur = total_impact_per_company(n_samples, n_a_mean, n_organizations, applicant_dist, organization_dist);
        if impact_cur.median > ((1 - quartile_frac) * impact_target.median + quartile_frac * impact_target.UQ)
            n_a_upper = n_a_mean
        elseif impact_cur.median < ((1 - quartile_frac) * impact_target.median + quartile_frac * impact_target.LQ)
            n_a_lower = n_a_mean
        else
            converged = true
            break
        end
    end

    if !converged
        @warn "Second phase didn't converge ($n_applicants, $n_organizations)"
    end

    n_a_mean -= n_applicants
    return DataFrame(:n_applicants => n_applicants, :n_organizations => n_organizations,
        :n_additional_applicants => n_a_mean, :frac_additional_applicants => n_a_mean / n_applicants)
end

## Utils

function ints_to_colors(vals::Vector{Int}, colormap::String="Oranges"; min_val=1, rev=true)
    res = Colors.colormap(colormap)[round.(Int, range(min_val, 100; length=length(unique(vals))))]
    if rev
        return reverse(res)
    else
        return res
    end
end

rotate_df(df_row::DataFrame; other_cols...) =
    DataFrame(:values => collect(df_row[1,:]), :type => names(df_row), other_cols...)

## Visualization utils

function savehtml(filepath::AbstractString, spec::VegaLite.AbstractVegaSpec)
    tmp_path = VegaLite.writehtml_full(spec)
    return cp(tmp_path, filepath, force=true)
end

function savehtml2(filepath::AbstractString, spec::VegaLite.AbstractVegaSpec)
    spec_str = VegaLite.getparams(spec) |> JSON.json
    content = """

    <title>VegaLite plot</title>
    <meta charset="UTF-8">

    """
    open(filepath, "w") do f
        print(f, content)
    end
end

function optimize_df(df::DataFrame, quant_vars::Vector{Symbol}, id_vars::Vector{Symbol})
    df = deepcopy(df[:, vcat(quant_vars, id_vars)])
    df[:, quant_vars] .= round.(df[:, quant_vars], sigdigits=4)
    return df
end

function optimize_df(df::DataFrame, quant_vars::Vector{Pair{Symbol, Symbol}}, id_vars::Vector{Pair{Symbol, Symbol}})
    quant_vars = Dict(quant_vars)
    id_vars = Dict(id_vars)
    df = deepcopy(df[:, vcat(keys(id_vars)..., keys(quant_vars)...)])
    df[:, collect(keys(quant_vars))] .= round.(df[:, collect(keys(quant_vars))], sigdigits=4)
    rename!(df, quant_vars..., id_vars...)
    return df
end

## Plots

function vl_plot_base(x::Pair{Symbol, String}, y::Pair{Symbol, String}, color::Pair{Symbol, String};
    extent::Int=40, width::Int=400, height::Int=300)
    return @vlplot(
        config={axisY={minExtent=extent, maxExtent=extent}},
        x={x[1], title=x[2], type="quantitative"},
        y={y[1], title=y[2], type="quantitative"},
        color={color[1], scale={scheme = "category20b"}, type="nominal", title=color[2]},
        opacity = {condition = {selection = "legend", value = 1}, value = 0.1},
        transform = [{filter = {selection = :Selectors}}],
        width=width, height=height
    )
end

vl_errorbar(y::Symbol=:UQ, y2::Symbol=:LQ) =
    @vlplot(mark={:errorbar, ticks=true}, y={y, type="quantitative", title=""}, y2={y2})

end

max_stars_repo_name: VPetukhov/EaImpactModeling

max_issues_repo_name: VPetukhov/EaImpactModeling

max_forks_repo_name: VPetukhov/EaImpactModeling


text:
using Flux: Chain, train!, params
using Flux.Losses:mse
using Flux.Data:DataLoader
using Flux.Optimise:ADAM
using IterTools:ncycle
using Statistics:mean

"""
    AE(encoder= Chain(),
       decoder = Chain(),
       batchsize= 32,
       epochs = 1,
       shuffle = false,
       partial = true,
       opt = ADAM(),
       loss = mse)

Calculate the anomaly score of an instance based on the reconstruction loss of an autoencoder, see [1] for an
explanation of auto encoders.

Parameters
----------
$_ae_params

$_ae_loss

Examples
--------
$(_score_unsupervised("AE"))

References
----------
[1] Aggarwal, Charu C. (2017): Outlier Analysis.
"""
MMI.@mlj_model struct AE <: UnsupervisedDetector
    encoder::Chain = Chain()
    decoder::Chain = Chain()
    batchsize::Integer = 32
    epochs::Integer = 1
    shuffle::Bool = false
    partial::Bool = true
    opt::Any = ADAM()
    loss::Function = mse
end

struct AEModel <: Model
    chain::Chain
end

function fit(detector::AE, X::Data)::Fit
    loader = DataLoader(X, batchsize=detector.batchsize, shuffle=detector.shuffle, partial=detector.partial)

    # Create the autoencoder
    model = Chain(detector.encoder, detector.decoder)

    # train the neural network model
    train!(x -> detector.loss(model(x), x), params(model), ncycle(loader, detector.epochs), detector.opt)

    scores = detector.loss(model(X), X, agg=instance_mean)
    Fit(AEModel(model), scores)
end

@score function score(detector::AE, model::Fit, X::Data)::Result
    detector.loss(model.chain(X), X, agg=instance_mean)
end


text:
# GetProperty
struct GetProperty{name} <: Function
end
@inline GetProperty(name::Symbol) = GetProperty{name}()

"""
    getproperty(name::Symbol)

Return a curried function equivalent to `x -> getproperty(x, name)`.

Internally, `name` is stored as a type parameter of a `GetProperty` object for the purpose
of constant propagation. See also `getproperties`.

# Example

Extract property `b` from a `NamedTuple`.

```julia
julia> nt = (a = 1, b = 2.0, c = false)
(a = 1, b = 2.0, c = false)

julia> getproperty(:b)(nt)
2.0
```
"""
@inline function Base.getproperty(name::Symbol)
    return GetProperty(name)
end

@inline (::GetProperty{name})(x) where {name} = getproperty(x, name)

# GetProperties
struct GetProperties{names} <: Function
end

@inline GetProperties(names::Tuple{Vararg{Symbol}}) = GetProperties{names}()

"""
    getproperties(names::Symbol...)

Return a function that extracts a set of properties with the given `names` from an object,
returning a new object with just those properties.

Internally, the `names` are stored as a type parameter of a `GetProperties` object for the
purpose of constant propagation. You may overload the `propertytype` function to control
the type of the object that is returned, which by default will be a `NamedTuple`. See also
`getproperty`.

# Example

Extract properties `a` and `c` from a `NamedTuple`.

```julia
julia> nt = (a = 1, b = 2.0, c = false)
(a = 1, b = 2.0, c = false)

julia> getproperties(:a, :c)(nt)
(a = 1, c = false)
```
"""
@inline function getproperties(names::Symbol...)
    return GetProperties(names)
end

@generated function (::GetProperties{names})(x) where {names}
    exprs = [:($n = getproperty(x, $(QuoteNode(n)))) for n in names]
    return Expr(:call, Expr(:call, :propertytype, :x), Expr(:tuple, exprs...))
end

"""
    propertytype(x)

Return a constructor for an object similar to `x` that can accept a `NamedTuple` with
arbitrary properties and support `getproperty`. Used for determining the return type of a
`getproperties` function. The defaults return type is `NamedTuple`.
"""
propertytype(x) = identity # the input is always a `NamedTuple`, and so this is the default propertytype

struct Compute{names, F} <: Function
    f::F
end

(::Type{Compute{names}})(f) where {names} = Compute{names, typeof(f)}(f)

@generated function (c::Compute{names})(x) where {names}
    exprs = [:(getproperty(x, $(QuoteNode(n)))) for n in names]
    return Expr(:call, :(c.f), exprs...)
end

const UNARY_OPERATORS = [:!, :-]
const BINARY_OPERATORS = [:+, :-, :*, :/, :\, ://, :÷, :×, :∈, :∉, :⊂, :⊆, :⊃, :⊇, :⊄, :⊈, :⊅, :⊉, :(==), :(===), :<, :<=, :>, :>=, :!=, :!==, :\, :≠, :≡, :≢, :≯ ,:≱, :≮, :≰]

"""
    @Compute(...)

The `@Compute` macro returns a function which performs a calculation on the properties of an
object, such as a `NamedTuple`.

The input expression is standard Julia code, with `\$` prepended to property names. For
example. if you want to refer to a property named `a` then use `\$a` in the expression.

# Example

```julia
julia> nt = (a = 1, b = 2.0, c = false)
(a = 1, b = 2.0, c = false)

julia> @Compute(\$a + \$b)(nt)
3.0
```
"""
macro Compute(expr)
    # We traverse the expression tree and locate the $ nodes, which we note the name of and
    # replace with getproperty calls
    names = Symbol[]

    # First we check if we can simplify the expression to make it more introsepctable
    if expr isa Expr
        # Simple column extraction
        if expr.head == :$ && length(expr.args) == 1 && expr.args[1] isa Symbol
            return :($(GetProperty{expr.args[1]}()))
        end

        if expr.head == :call && countdollars(expr.args[1]) == 0
            # Single argument function
            if length(expr.args) == 2 && expr.args[2] isa Expr && expr.args[2].head == :$ && length(expr.args[2].args) == 1 && expr.args[2].args[1] isa Symbol
                f = expr.args[1]
                name = expr.args[2].args[1]
                return :(Compute{($(QuoteNode(name)),)}($f))
            end

            # Two argument functions
            if length(expr.args) == 3
                # Binary function
                if expr.args[2] isa Expr && expr.args[2].head == :$ && length(expr.args[2].args) == 1 && expr.args[2].args[1] isa Symbol && expr.args[3] isa Expr && expr.args[3].head == :$ && length(expr.args[3].args) == 1 && expr.args[3].args[1] isa Symbol
                    f = expr.args[1]
                    name1 = expr.args[2].args[1]
                    name2 = expr.args[3].args[1]
                    return :(Compute{($(QuoteNode(name1)),$(QuoteNode(name2)))}($f))
                end

                # Fix1
                if countdollars(expr.args[2]) == 0 && expr.args[3] isa Expr && expr.args[3].head == :$ && length(expr.args[3].args) == 1 && expr.args[3].args[1] isa Symbol
                    f = expr.args[1]
                    x = expr.args[2]
                    name = expr.args[3].args[1]
                    return :(Compute{($(QuoteNode(name)),)}(Base.Fix1($f, $x)))
                end

                # Fix2
                if expr.args[2] isa Expr && expr.args[2].head == :$ && length(expr.args[2].args) == 1 && expr.args[2].args[1] isa Symbol && countdollars(expr.args[3]) == 0
                    f = expr.args[1]
                    name = expr.args[2].args[1]
                    x = expr.args[3]
                    return :(Compute{($(QuoteNode(name)),)}(Base.Fix2($f, $x)))
                end
            end
        end

        # Unary operator
        if expr.head ∈ UNARY_OPERATORS && length(expr.args) == 1 && expr.args[1].head == :$ && length(expr.args[1].args) == 1 && expr.args[1].args[1] isa Symbol
            f = expr.head
            name = expr.args[1].args[1]
            return :(Compute{($(QuoteNode(name)),)}($f))
        end

        # Binary operator
        if expr.head ∈ BINARY_OPERATORS && length(expr.args) == 2
            # Binary function
            if expr.args[1].head == :$ && length(expr.args[1].args) == 1 && expr.args[1].args[1] isa Symbol && expr.args[2].head == :$ && length(expr.args[2].args) == 1 && expr.args[2].args[1] isa Symbol
                f = expr.head
                name1 = expr.args[1].args[1]
                name2 = expr.args[2].args[1]
                return :(Compute{($(QuoteNode(name1)),$(QuoteNode(name2)))}($f))
            end

            # Fix1
            if countdollars(expr.args[1]) == 0 && expr.args[2].head == :$ && length(expr.args[2].args) == 1 && expr.args[2].args[1] isa Symbol
                f = expr.head
                x = expr.args[1]
                name = expr.args[2].args[1]
                return :(Compute{($(QuoteNode(name)),)}(Base.Fix1($f, $x)))
            end

            # Fix2
            if expr.args[1].head == :$ && length(expr.args[1].args) == 1 && expr.args[1].args[1] isa Symbol && countdollars(expr.args[2]) == 0
                f = expr.head
                name = expr.args[1].args[1]
                x = expr.args[2]
                return :(Compute{($(QuoteNode(name)),)}(Base.Fix2($f, $x)))
            end
        end
    end

    # In the general case, we make a new closure
    compute_expr!(names, expr)
    return :(Compute{$(tuple(names...))}((x...) -> $expr))
end

countdollars(expr) = countdollars(expr, 0)
countdollars(expr, n::Int) = n
function countdollars(expr::Expr, n::Int)
    if expr.head == :$
        @assert length(expr.args) == 1 && expr.args[1] isa Symbol
        n += 1
    else
        n += sum(map(countdollars, expr.args))
    end
    return n
end

function compute_expr!(names::Vector{Symbol}, expr::Expr)
    if expr.head == :$
        @assert length(expr.args) == 1
        n = expr.args[1]::Symbol

        i = 0
        for j = 1:length(names)
            if names[j] == n
                i = j
                break
            end
        end
        if i == 0
            push!(names, n)
            i = length(names)
        end

        expr.head = :ref
        expr.args = [:x, i]
    else
        compute_expr!(names, expr.head)
        foreach(ex -> compute_expr!(names, ex), expr.args)
    end

    return nothing
end

compute_expr!(names::Vector{Symbol}, expr::Any) = nothing

struct Select{names, Fs <: Tuple} <: Function
    fs::Fs
end

(::Type{Select{names}})(fs::Tuple) where {names} = Select{names, typeof(fs)}(fs)

@generated function (s::Select{names})(x) where {names}
    exprs = [:($(names[i]) = (s.fs[$i])(x)) for i in 1:length(names)]
    return Expr(:call, Expr(:call, :propertytype, :x), Expr(:tuple, exprs...))
end

"""
    @Select(...)

The `@Select` macro returns a function which performs an arbitrary transformation of the
properties of an object, such as a `NamedTuple`.

The input expression is a comma-seperated list of `lhs = rhs` pairs. The `lhs` is the name
of the new property to calculate. The `rhs is  standard Julia code, with `\$` prepended to
input property names. For example. if you want to rename an input property `a` to be called
`b`, use `@Select(b = \$a)`.

As a special case, if a property is to be simply replicated the `= rhs` part can be dropped,
for example `@Select(a)` is synomous with `@Select(a = \$a)`.

# Example

```julia
julia> nt = (a = 1, b = 2.0, c = false)
(a = 1, b = 2.0, c = false)

julia> @Select(a, sum_a_b = \$a + \$b)(nt)
(a = 1, sum_a_b = 3.0)
```
"""
macro Select(exprs...)
    # For each express we extract the output property name and the associated `Compute`.
    names = Symbol[]
    fs = Expr[]

    if exprs isa Tuple{Vararg{Symbol}}
        return :($(GetProperties{exprs}()))
    end

    for expr in exprs
        if expr isa Symbol
            push!(names, expr)
            push!(fs, :(GetProperty{$(QuoteNode(expr))}()))
        elseif expr isa Expr && expr.head == :(=)
            @assert length(expr.args) == 2
            push!(names, expr.args[1])
            internal_names = Symbol[]
            compute_expr!(internal_names, expr.args[2])
            push!(fs, :(Compute{$(tuple(internal_names...))}((x...) -> $(expr.args[2]))))
        else
            error("Bad input to @select")
        end
    end

    @assert length(unique(names)) == length(names)

    return :(Select{$(tuple(names...))}(tuple($(fs...))))
end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test

println("Scalar g")
A = [-1.0 0.0; 0.0 -0.5]
u0 = [1.0, 1.0]; tspan = (0.0,1.0)
_f = (u,p,t) -> t*(A*u)
_g = (u,p,t) -> 1.0
prob = SDEProblem(SDEFunction(_f, _g), _g, u0, tspan)
integrator = init(prob, SKenCarp(); adaptive=false, dt=0.01)
step!(integrator)
@test_broken solve(prob, SOSRI(); adaptive=false, dt=0.01) isa RODESolution
solve(prob, SOSRA(); adaptive=false, dt=0.01)
solve(prob, EM(); adaptive=false, dt=0.01)
@test_broken solve(prob, RKMil(); adaptive=false, dt=0.01) isa RODESolution
solve(prob, SRIW1(); adaptive=false, dt=0.01)
@test_broken solve(prob, SOSRI2(); adaptive=false, dt=0.01) isa RODESolution
solve(prob, SOSRA2(); adaptive=false, dt=0.01)

println("Vector g")
_g = (u,p,t) -> [1.0, 1.0]
prob = SDEProblem(SDEFunction(_f, _g), _g, u0, tspan)
println("Implicit EM")
integrator = init(prob, ImplicitEM(); adaptive=false, dt=0.01)
step!(integrator)
println("SKenCarp")
integrator = init(prob, SKenCarp(); adaptive=false, dt=0.01)
step!(integrator)
solve(prob, SOSRI(); adaptive=false, dt=0.01)
solve(prob, SOSRA(); adaptive=false, dt=0.01)
solve(prob, EM(); adaptive=false, dt=0.01)
solve(prob, RKMil(); adaptive=false, dt=0.01)
solve(prob, SOSRI2(); adaptive=false, dt=0.01)
solve(prob, SOSRA2(); adaptive=false, dt=0.01)


text:
# This file contains functions for creating lattices

"""
    reggeomlattice(n, m)
Create a regular geometry-lattice of ``n`` rows and ``m`` columns.

# Optional Arguments
- `gtype=:polygon`: geometry. :polygon or :point.
- `direction=:rightdown`: direction. :rightup, :rightdown, :leftdown or :leftup.
"""
function reggeomlattice(n::Int, m::Int; gtype::Symbol = :polygon, direction::Symbol = :rightdown)
    n > 0 || throw(ArgumentError("The number of rows must be greater than 0"))
    m > 0 || throw(ArgumentError("The number of columns must be greater than 0"))

    (direction == :rightup) || (direction == :rightdown) ||
    (direction == :leftdown) || (direction == :leftup) ||
        throw(ArgumentError("`direction` must be :rightup, :rightdown, :leftdown or :leftup"))

    (gtype == :polygon) || (gtype == :point) || throw(ArgumentError("`gtype` must be :polygon or :point"))

    if direction == :rightup
        iiter = 1:n
        jiter = 1:m
    elseif direction == :rightdown
        iiter = reverse(1:n)
        jiter = 1:m
    elseif direction == :leftdown
        iiter = reverse(1:n)
        jiter = reverse(1:m)
    elseif direction == :leftup
        iiter = 1:n
        jiter = reverse(1:m)
    end

    if gtype == :polygon
        geometries = MultiPolygon[]
    elseif gtype == :point
        geometries = Point[]
    end

    for i in iiter
        for j in jiter
            if gtype == :polygon
                push!(geometries,
                    MultiPolygon(
                        Polygon([
                            [[0.0 + j, 0.0 + i],
                            [1.0 + j, 0.0 + i],
                            [1.0 + j, 1.0 + i],
                            [0.0 + j, 1.0 + i] ]])))
            elseif gtype == :point
                push!(geometries, Point([0.0 + j, 0.0 + i]))
            end
        end
    end

    return geometries
end


text:
# calculate probability by
# gaussian distribution pdf
# plot as 2d-contour
# 3 different contours

module MultiGaussDist
    using Plots, Base.Iterators, LinearAlgebra
    pyplot()

    function pdf(x, mu, cov_mat)
        # determinant
        det_mat = det(cov_mat)

        diff_x = x - mu
        return exp(-((diff_x'/cov_mat)*diff_x)/2)/(2*pi*sqrt(det_mat))
    end

    function main(is_test=false)
        # x-y grid
        vx = 0:200
        vy = 0:100

        # mu, covariance matrix
        # contour a
        mu_a = [50; 50]
        cov_a = [50 0; 0 100]
        # contour b
        mu_b = [100; 50]
        cov_b = [125 0; 0 25]
        # contour c
        mu_c = [150; 50]
        cov_c = [100 -25*sqrt(3); -25*sqrt(3) 50]

        # probability density
        # contour a
        z_a = [pdf([x; y], mu_a, cov_a) for x in vx, y in vy]
        # contour b
        z_b = [pdf([x; y], mu_b, cov_b) for x in vx, y in vy]
        # contour c
        z_c = [pdf([x; y], mu_c, cov_c) for x in vx, y in vy]

        # plot contours
        contour(vx, vy, z_a', c=:haline, aspect_ratio=:equal)
        contour!(vx, vy, z_b', c=:haline, aspect_ratio=:equal)
        contour!(vx, vy, z_c', c=:haline, aspect_ratio=:equal)

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/multi_dim_gauss_dist/multiple_gauss_dist/multiple_gauss_dist.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module MyPkgDemo

export add_nums
export cat_str

include("add_nums.jl")
include("cat_str.jl")

end

max_stars_repo_name: junedwx/MyPkgDemo

max_issues_repo_name: junedwx/MyPkgDemo


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ d6bc6cde-7686-11eb-0a34-f1e7c9841275
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ f030f70c-7686-11eb-3c72-c59f85af6c83
md"""# Iteration

This chapter is about iteration, which is the ability to run a block of statements repeatedly. We saw a kind of iteration, using recursion, in “Recursion”. We saw another kind, using a `for` loop, in “Simple Repetition”. In this chapter we’ll see yet another kind, using a `while` statement. But first I want to say a little more about variable assignment.
"""

# ╔═╡ b464d026-7687-11eb-2e68-99bb4a60f392
md"""## Reassignment

As you may have discovered, it is legal to make more than one assignment to the same variable. A new assignment makes an existing variable refer to a new value (and stop referring to the old value):

```julia
julia> x = 5
5
julia> x = 7
7
```

The first time we display `x`, its value is `5`; the second time, its value is `7`.

Figure 7-1 shows what *reassignment* looks like in a state diagram.
"""

# ╔═╡ f001d996-7687-11eb-1215-5940f40d3fed
Drawing(width=720, height=50) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
    rect(x=310, y=10, width=100, height=40, fill="rgb(242, 242, 242)", stroke="black")
	text(x=330, y=35, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("x")
	end
	text(x=390, y=25, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("5")
	end
	text(x=390, y=45, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("7")
	end
	line(x1=335, y1=30, x2=370, y2=23, stroke="black", stroke_dasharray="5,5", marker_end="url(#arrow)")
	line(x1=335, y1=32, x2=370, y2=37, stroke="black", marker_end="url(#arrow)")
end

# ╔═╡ 55c7f344-768d-11eb-2f35-e52d557f4c61
md"*7-1. State diagram.*"

# ╔═╡ 67411830-768d-11eb-1cf0-f9011222c68d
md"""At this point I want to address a common source of confusion. Because Julia uses the equals sign (`=`) for assignment, it is tempting to interpret a statement like `a = b` as a mathematical proposition of equality; that is, the claim that `a` and `b` are equal. But this interpretation is wrong.

First, equality is a symmetric relationship and assignment is not. For example, in mathematics, if ``a = 7`` then ``7 = a``. But in Julia, the statement `a = 7` is legal and `7 = a` is not.

Also, in mathematics, a proposition of equality is either true or false for all time. If ``a = b`` now, then ``a`` will always equal ``b``. In Julia, an assignment statement can make two variables equal, but they don’t have to stay that way:

```julia
julia> a = 5
5
julia> b = a # a and b are now equal
5
julia> a = 3 # a and b are no longer equal
3
julia> b
5
```

The third line changes the value of `a `but does not change the value of `b`, so they are no longer equal.

!!! danger
    Reassigning variables is often useful, but you should use it with caution. If the values of variables change frequently, it can make the code difficult to read and debug.

    It is illegal to define a function that has the same name as a previously defined variable.
"""

# ╔═╡ 13602200-768e-11eb-389b-fb726834e96f
md"""## Updating Variables

A common kind of reassignment is an update, where the new value of the variable depends on the old:

```julia
julia> x = x + 1
8
```

or

```julia
julia> x += 1
8
```

This means “get the current value of `x`, add 1, and then update `x` with the new value.”

If you try to update a variable that doesn’t exist, you get an error, because Julia evaluates the right side before it assigns a value to x:

```julia
julia> y = y + 1
ERROR: UndefVarError: y not defined
```

Before you can update a variable, you have to initialize it, usually with a simple assignment:

```julia
julia> y = 0
0
julia> y = y + 1
1
```

Updating a variable by adding 1 is called an *increment*; subtracting 1 is called a *decrement*.
"""

# ╔═╡ 658792ba-768f-11eb-0e59-6d3a8c392c59
md"""## The `while` Statement

Computers are often used to automate repetitive tasks. Repeating identical or similar tasks without making errors is something that computers do well and people do poorly. In a computer program, repetition is also called *iteration*.

We have already seen two functions, `countdown` and `printn`, that iterate using recursion. Because iteration is so common, Julia provides language features to make it easier. One is the `for` statement we saw in “Simple Repetition”. We’ll get back to that later.

Another is the *`while` statement*. Here is a version of `countdown` that uses a `while` statement:

```julia
function countdown(n)
	while n > 0
		print(n, " ")
		n -= 1
	end
	println("Blastoff!")
end
```

You can almost read the `while` statement as if it were English. It means, “While `n` is greater than `0`, display the value of `n` and then decrement `n`. When you get to `0`, display the word `"Blastoff!"`”

More formally, here is the flow of execution for a `while` statement:

1. Determine whether the condition is true or false.
2. If false, exit the `while` statement and continue execution at the next statement.
3. If the condition is true, run the body and then go back to step 1.

This type of flow is called a loop because the third step loops back around to the top.

The body of the loop should change the value of one or more variables so that the condition becomes false eventually and the loop terminates. Otherwise, the loop will repeat forever, which is called an *infinite loop*. An endless source of amusement for computer scientists is the observation that the directions on shampoo bottles, “Lather, rinse, repeat,” are an infinite loop.

In the case of `countdown`, we can prove that the loop terminates: if `n` is `0` or negative, the loop never runs. Otherwise, `n` gets smaller each time through the loop, so eventually we have to get to `0`.

For some other loops, it is not so easy to tell. For example:

```julia
function seq(n)
	while n != 1
		println(n)
		if n % 2 == 0   # n is even
			n /= 2
		else
			n = 3n + 1  # n is odd
		end
	end
end
```

The condition for this loop is `n != 1`, so the loop will continue until `n` is `1`, which makes the condition false.

Each time through the loop, the program outputs the value of `n` and then checks whether it is even or odd. If it is even, `n` is divided by 2. If it is odd, the value ofnis replaced with `3n + 1`. For example, if the argument passed to `seq` is `3`, the resulting values of `n` are `3`, `10`, `5.0`, `16.0`, `8.0`, `4.0`, `2.0`, `1.0`.

Since `n` sometimes increases and sometimes decreases, there is no obvious proof that `n` will ever reach `1`, or that the program terminates. For some particular values of `n`, we can prove termination. For example, if the starting value is a power of two, `n` will be even every time through the loop until it reaches `1`. The previous example ends with such a sequence, starting with `16`.

The hard question is whether we can prove that this program terminates for all positive values of `n`. So far, no one has been able to prove it or disprove it!
"""

# ╔═╡ f51d8d5c-7690-11eb-3eea-e5a92768f7e6
md"""#### Exercise 7-1

Rewrite the function `printn` from “Recursion” using iteration instead of recursion.
"""

# ╔═╡ 87bc44b0-76a4-11eb-3e5a-13df6d96717d
md"""!!! languages
    Both Python and the C programming language have a `while` statement:

    ```python
    def countdown(n):
        while n > 0:
            print(n, end=' ')
            n = n - 1

        print('Blastoff!')
    ```

    The conventions for each language are always similar.

    ```c
    void countdown(int n) {
        while (n > 0) {
            printf("%i ", n);
            n -= 1;
        }
        printf("Blastoff!\n");
    }
    ```
"""

# ╔═╡ 5adbadaa-76a6-11eb-34b0-0b02b1cd9f95
md"""## `break`

Sometimes you don’t know it’s time to end a loop until you get halfway through the body. In that case you can use the *`break` statement* to jump out of the loop.

For example, suppose you want to take input from the user until they type **`done`**. You could write:

```julia
while true
    print("> ")
    line = readline()
    if line == "done"
		break
	end
	println(line)
end
println("Done!")
```

The loop condition is `true`, which is always true, so the loop runs until it hits the
`break` statement.

Each time through, it prompts the user with an angle bracket. If the user types **`done`**, the `break` statement exits the loop. Otherwise, the program echoes whatever the user types and goes back to the top of the loop. Here’s a sample run:

```
> not done
not done
> done
Done!
```

This way of writing `while` loops is common because you can check the condition anywhere in the loop (not just at the top) and you can express the stop condition affirmatively (“stop when this happens”) rather than negatively (“keep going until that happens”).
"""

# ╔═╡ 0aa7430a-76a7-11eb-1626-05ef9338e919
md"""## `continue`

The `break` statement exits the loop. When a *`continue` statement* is encountered inside a loop, control jumps to the beginning of the loop for the next iteration, skipping the execution of statements inside the body of the loop for the current iteration. For example, this

```julia
for i in 1:10
	if i % 3 == 0
		continue
	end
	print(i, " ")
end
```

outputs

```
1 2 4 5 7 8 10
```

If `i` is divisible by `3`, the continue statement stops the current iteration and the nex iteration starts. Only the numbers in the range 1 to 10 not divisible by 3 are printed.
"""

# ╔═╡ c6aaba90-76a8-11eb-1e6f-495fb7759e5d
md"""!!! languages
    Both Python and the C programming language have `break` and `continue` statements.
"""

# ╔═╡ e74d8976-76a8-11eb-199d-c33a925a4d7f
md"""## Square Roots

Loops are often used in programs that compute numerical results by starting with an approximate answer and iteratively improving it.

For example, one way of computing square roots is Newton’s method. Suppose that you want to know the square root of a. If you start with almost any estimate, x, you can compute a better estimate with the following formula:

```math
y=\frac{1}{2}\left(x+\frac{a}{x}\right)
```

For example, if ``a`` is ``4`` and ``x`` is ``3``:

```julia
julia> a = 4
4
julia> x = 3
3
julia> y = 0.5(x + a/x)
2.1666666666666665
```

The result is closer to the correct answer (``\sqrt 4 = 2``). If we repeat the process with the new estimate, it gets even closer:

```julia
julia> x = y
2.1666666666666665

julia> y = 0.5(x + a/x)
2.0064102564102564
```

After a few more updates, the estimate is almost exact:

```julia
julia> x = y
2.0064102564102564

julia> y = 0.5(x + a/x)
2.0000102400262145

julia> x = y
2.0000102400262145

julia> y = 0.5(x + a/x)
2.0000000000262146
```

In general we don’t know ahead of time how many steps it takes to get to the right answer, but we know when we get there because the estimate stops changing:

```julia
julia> x = y
2.0000000000262146

julia> y = 0.5(x + a/x)
2.0

julia> x = y
2.0

julia> y = 0.5(x + a/x)
2.0
```

When `y == x`, we can stop. Here is a loop that starts with an initial estimate, `x`, and improves it until it stops changing:

```julia
while true
	println(x)
	y = 0.5(x + a/x)
	if y == x
		break
	end
	x = y
end
```

For most values of a this works fine, but in general it is dangerous to test float equality. Floating-point values are only approximately right: most rational numbers, like ``\frac{1}{3}``, and irrational numbers, like ``\sqrt 2``, can’t be represented exactly with a `Float64`.

Rather than checking whether `x` and `y` are exactly equal, it is safer to use the built-in function `abs` to compute the absolute value, or magnitude, of the difference between them:

```julia
if abs(y - x) < ε
	break
end
```

where `ε` (**`\varepsilon TAB`**) has a value like `0.0000001` that determines how close is close enough.
"""

# ╔═╡ c5327386-76af-11eb-1311-0f0f5cd2c507
md"""## Algorithms

Newton’s method is an example of an *algorithm*: it is a mechanical process for solving a category of problems (in this case, computing square roots).

To understand what an algorithm is, it might help to start with something that is not an algorithm. When you learned to multiply single-digit numbers, you probably memorized the multiplication table. In effect, you memorized 100 specific solutions. That kind of knowledge is not algorithmic.

But if you were “lazy”, you might have learned a few tricks. For example, to find the product of ``n`` and ``9``, you can write ``n−1`` as the first digit and ``10−n`` as the second digit. This trick is a general solution for multiplying any single-digit number by ``9``. That’s an algorithm!

Similarly, the techniques you learned for addition with carrying, subtraction with borrowing, and long division are all algorithms. One of the characteristics of algo‐ rithms is that they do not require any intelligence to carry out. They are mechanical processes where each step follows from the last according to a simple set of rules.

Executing algorithms is boring, but designing them is interesting, intellectually challenging, and a central part of computer science.

Some of the things that people do naturally, without difficulty or conscious thought, are the hardest to express algorithmically. Understanding natural language is a good example. We all do it, but so far no one has been able to explain how we do it, at least not in the form of an algorithm.
"""

# ╔═╡ 0c60279e-76b0-11eb-0687-bf8d44d9cd47
md"""## Debugging
As you start writing bigger programs, you might find yourself spending more time debugging. More code means more chances to make an error and more places for bugs to hide.

One way to cut your debugging time is “debugging by bisection.” For example, if there are 100 lines in your program and you check them one at a time, it will take 100 steps.

Instead, try to break the problem in half. Look at the middle of the program, or near it, for an intermediate value you can check. Add a print statement (or something else that has a verifiable effect) and run the program.

If the midpoint check is incorrect, there must be a problem in the first half of the pro‐ gram. If it is correct, the problem is in the second half.

Every time you perform a check like this, you halve the number of lines you have to search. In theory, after six steps (which is significantly fewer than 100), you’ll be down to one or two lines of code.

In practice, it is not always clear what the “middle of the program” is and not always possible to check it. It doesn’t make sense to count lines and find the exact midpoint. Instead, think about places in the program where there might be errors and places where it is easy to put a check. Then choose a spot where you think the chances are about the same that the bug is before or after the check.
"""

# ╔═╡ 2b84de3a-76b0-11eb-184e-fbf49560e122
md"""## Glossary

*reassignment*:
Assigning a new value to a variable that already exists.

*update*:
An assignment where the new value of the variable depends on the old.

*initialization*:
An assignment that gives an initial value to a variable that will be updated.

*increment*:
An update that increases the value of a variable (often by 1).

*decrement*:
An update that decreases the value of a variable.

*iteration*:
Repeated execution of a set of statements using either a recursive function call or a loop.

*`while` statement*:
A statement that allows iterations controlled by a condition.

*infinite loop*:
A loop in which the terminating condition is never satisfied.

*`break` statement*:
A statement allowing you to jump out of a loop.

*`continue` statement*:
A statement inside a loop that jumps to the beginning of the loop for the next iteration.

*algorithm*:
A general process for solving a category of problems.
"""

# ╔═╡ 7e6a42ea-76b0-11eb-0587-0b2e3932b601
md"""## Exercises

#### Exercise 7-2

Copy the loop from “Square Roots” and encapsulate it in a function called `mysqrt` that takes `a` as a parameter, chooses a reasonable value of `x`, and returns an estimate of the square root of `a`.

To test it, write a function named `testsquareroot` that prints a table like this:

```
a   mysqrt             sqrt               diff
-   ------             ----               ----
1.0 1.0                1.0                0.0
2.0 1.414213562373095  1.4142135623730951 2.220446049250313e-16
3.0 1.7320508075688772 1.7320508075688772 0.0
4.0 2.0                2.0                0.0
5.0 2.23606797749979   2.23606797749979   0.0
6.0 2.449489742783178  2.449489742783178  0.0
7.0 2.6457513110645907 2.6457513110645907 0.0
8.0 2.82842712474619   2.8284271247       4.440892098500626e-16
9.0 3.0                3.0                0.0
```

The first column is a number, `a`; the second column is the square root of a computed with `mysqrt`; the third column is the square root computed by `sqrt`; and the fourth column is the absolute value of the difference between the two estimates.
"""

# ╔═╡ f24e401a-76b0-11eb-1530-930c3df061fc
md"""#### Exercise 7-3

The built-in function `Meta.parse` takes a string and transforms it into an expression. This expression can be evaluated in Julia with the function `Core.eval`. For example:

```julia
julia> expr = Meta.parse("1+2*3")
:(1 + 2 * 3)
julia> eval(expr)
7
julia> expr = Meta.parse("sqrt(π)")
:(sqrt(π))
julia> eval(expr)
1.7724538509055159
```

Write a function called `evalloop` that iteratively prompts the user, takes the resulting input and evaluates it using `eval`, and prints the result. It should continue until the user enters **`done`**, and then return the value of the last expression it evaluated.
"""

# ╔═╡ 2c14de1c-76b1-11eb-03f7-b73b33ae0d24
md"""## Exercise 7-4

The mathematician Srinivasa Ramanujan found an infinite series that can be used to generate a numerical approximation of ``\frac{1}{π}``:

```math
\frac{1}{π}=\frac{2\sqrt 2}{9801}\sum_{k=0}^\infty\frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}
```

Write a function called `estimatepi` that uses this formula to compute and return an estimate of ``π``. It should use a while loop to compute terms of the summation until the last term is smaller than `1e-15` (which is Julia notation for ``10^{−15}``). You can check the result by comparing it to `π`.
"""

# ╔═╡ Cell order:
# ╟─d6bc6cde-7686-11eb-0a34-f1e7c9841275
# ╟─f030f70c-7686-11eb-3c72-c59f85af6c83
# ╟─b464d026-7687-11eb-2e68-99bb4a60f392
# ╠═f001d996-7687-11eb-1215-5940f40d3fed
# ╟─55c7f344-768d-11eb-2f35-e52d557f4c61
# ╟─67411830-768d-11eb-1cf0-f9011222c68d
# ╟─13602200-768e-11eb-389b-fb726834e96f
# ╟─658792ba-768f-11eb-0e59-6d3a8c392c59
# ╟─f51d8d5c-7690-11eb-3eea-e5a92768f7e6
# ╟─87bc44b0-76a4-11eb-3e5a-13df6d96717d
# ╟─5adbadaa-76a6-11eb-34b0-0b02b1cd9f95
# ╟─0aa7430a-76a7-11eb-1626-05ef9338e919
# ╟─c6aaba90-76a8-11eb-1e6f-495fb7759e5d
# ╟─e74d8976-76a8-11eb-199d-c33a925a4d7f
# ╟─c5327386-76af-11eb-1311-0f0f5cd2c507
# ╟─0c60279e-76b0-11eb-0687-bf8d44d9cd47
# ╟─2b84de3a-76b0-11eb-184e-fbf49560e122
# ╟─7e6a42ea-76b0-11eb-0587-0b2e3932b601
# ╟─f24e401a-76b0-11eb-1530-930c3df061fc
# ╟─2c14de1c-76b1-11eb-03f7-b73b33ae0d24


text:
@enum AggOperations begin

    MEAN = 0
    MEDIAN = 1
    EACH = 2

end

@enum FilterTypes begin

    BUTTERWORTH = 0
    CHEBYSHEV_TYPE_1 = 1
    BESSEL = 2

end

function perform_lowpass(data, sampling_rate, cutoff, order, filter_type, ripple)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_lowpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    elseif Sys.isapple()
        ec = ccall((:perform_lowpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    else
        ec = ccall((:perform_lowpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_lowpass ", ec), ec))
    end
    return
end

function perform_highpass(data, sampling_rate, cutoff, order, filter_type, ripple)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_highpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    elseif Sys.isapple()
        ec = ccall((:perform_highpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    else
        ec = ccall((:perform_highpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, cutoff, order, filter_type, ripple)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_highpass ", ec), ec))
    end
    return
end

function perform_bandpass(data, sampling_rate, center_freq, band_width, order, filter_type, ripple)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_bandpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    elseif Sys.isapple()
        ec = ccall((:perform_bandpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    else
        ec = ccall((:perform_bandpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_bandpass ", ec), ec))
    end
    return
end

function perform_bandstop(data, sampling_rate, center_freq, band_width, order, filter_type, ripple)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_bandstop, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    elseif Sys.isapple()
        ec = ccall((:perform_bandstop, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    else
        ec = ccall((:perform_bandstop, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
        	data, length(data), sampling_rate, center_freq, band_width, order, filter_type, ripple)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_bandstop ", ec), ec))
    end
    return
end

function perform_rolling_filter(data, period, operation)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_rolling_filter, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
        	data, length(data), period, operation)
    elseif Sys.isapple()
        ec = ccall((:perform_rolling_filter, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
        	data, length(data), period, operation)
    else
        ec = ccall((:perform_rolling_filter, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
        	data, length(data), period, operation)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_rolling_filter ", ec), ec))
    end
    return
end

function perform_wavelet_denoising(data, wavelet, decomposition_level)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_wavelet_denoising, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
        	data, length(data), wavelet, decomposition_level)
    elseif Sys.isapple()
        ec = ccall((:perform_wavelet_denoising, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
        	data, length(data), wavelet, decomposition_level)
    else
        ec = ccall((:perform_wavelet_denoising, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
        	data, length(data), wavelet, decomposition_level)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_wavelet_denoising ", ec), ec))
    end
    return
end

function perform_downsampling(data, period, operation)
	len = Integer(floor(length(data) / period))
	downsampled_data = Vector{Float64}(undef, len)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_downsampling, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
        	data, length(data), period, operation, downsampled_data)
    elseif Sys.isapple()
        ec = ccall((:perform_downsampling, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
        	data, length(data), period, operation, downsampled_data)
    else
        ec = ccall((:perform_downsampling, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
        	data, length(data), period, operation, downsampled_data)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_downsampling ", ec), ec))
    end
    downsampled_data
end

function write_file(data, file_name, file_mode)
	shape = size(data)
	flatten = reshape(transpose(data), (1, shape[1] * shape[2]))
	flatten = copy(flatten)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:write_file, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
        	flatten, shape[1], shape[2], file_name, file_mode)
    elseif Sys.isapple()
        ec = ccall((:write_file, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
        	flatten, shape[1], shape[2], file_name, file_mode)
    else
        ec = ccall((:write_file, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
        	flatten, shape[1], shape[2], file_name, file_mode)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in write_file ", ec), ec))
    end
    return
end

function read_file(file_name)
	num_elements = Vector{Cint}(undef, 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_num_elements_in_file, "DataHandler.dll"), Cint, (Ptr{UInt8}, Ptr{Cint}),
        		file_name, num_elements)
    elseif Sys.isapple()
        ec = ccall((:get_num_elements_in_file, "libDataHandler.dylib"), Cint, (Ptr{UInt8}, Ptr{Cint}),
        		file_name, num_elements)
    else
        ec = ccall((:get_num_elements_in_file, "libDataHandler.so"), Cint, (Ptr{UInt8}, Ptr{Cint}),
        		file_name, num_elements)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_num_elements_in_file ", ec), ec))
    end

    data = Vector{Float64}(undef, num_elements[1])
    num_cols = Vector{Cint}(undef, 1)
    num_rows = Vector{Cint}(undef, 1)

    if Sys.iswindows()
        ec = ccall((:read_file, "DataHandler.dll"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
        		data, num_rows, num_cols, file_name, num_elements[1])
    elseif Sys.isapple()
        ec = ccall((:read_file, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
        		data, num_rows, num_cols, file_name, num_elements[1])
    else
        ec = ccall((:read_file, "libDataHandler.so"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
        		data, num_rows, num_cols, file_name, num_elements[1])
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in read_file ", ec), ec))
    end
    value = transpose(reshape(data, (num_cols[1], num_rows[1])))
    value
end

function perform_wavelet_transform(data, wavelet, decomposition_level)
	wavelet_coeffs = Vector{Float64}(undef, length(data) + 2 * (40 + 1))
	lengths = Vector{Cint}(undef, decomposition_level + 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_wavelet_transform, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
        	data, length(data), wavelet, decomposition_level, wavelet_coeffs, lengths)
    elseif Sys.isapple()
        ec = ccall((:perform_wavelet_transform, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
        	data, length(data), wavelet, decomposition_level, wavelet_coeffs, lengths)
    else
        ec = ccall((:perform_wavelet_transform, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
        	data, length(data), wavelet, decomposition_level, wavelet_coeffs, lengths)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_wavelet_transform ", ec), ec))
    end
    wavelet_coeffs[1:sum(lengths)], lengths
end

function perform_inverse_wavelet_transform(wavelet_output, original_data_len, wavelet, decomposition_level)
	original_data = Vector{Float64}(undef, original_data_len)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_inverse_wavelet_transform, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
        	wavelet_output[1], original_data_len, wavelet, decomposition_level, wavelet_output[2], original_data)
    elseif Sys.isapple()
        ec = ccall((:perform_inverse_wavelet_transform, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
        	wavelet_output[1], original_data_len, wavelet, decomposition_level, wavelet_output[2], original_data)
    else
        ec = ccall((:perform_inverse_wavelet_transform, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
        	wavelet_output[1], original_data_len, wavelet, decomposition_level, wavelet_output[2], original_data)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_inverse_wavelet_transform ", ec), ec))
    end
    original_data
end

function perform_fft(data)

	function is_power_of_two(value)
		(value != 0) && (value & (value - 1) == 0)
	end

	if !is_power_of_two(length(data))
		throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
	end

	temp_re = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
	temp_im = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
	res = Vector{Complex}(undef, Integer(length(data) / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_fft, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}),
        	data, length(data), temp_re, temp_im)
    elseif Sys.isapple()
        ec = ccall((:perform_fft, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}),
        	data, length(data), temp_re, temp_im)
    else
        ec = ccall((:perform_fft, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}),
        	data, length(data), temp_re, temp_im)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_fft ", ec), ec))
    end
    for i in 1:Integer(length(data) / 2) + 1
    	res[i] = Complex(temp_re[i], temp_im[i])
    end
    res
end

function perform_ifft(data)

	temp_re = Vector{Float64}(undef, length(data))
	temp_im = Vector{Float64}(undef, length(data))
	res = Vector{Float64}(undef, 2 * (length(data) - 1))

	for i in 1:length(data)
		temp_re[i] = data[i].re
		temp_im[i] = data[i].im
    end
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_ifft, "DataHandler.dll"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
        	temp_re, temp_im, length(res), res)
    elseif Sys.isapple()
        ec = ccall((:perform_ifft, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
        	temp_re, temp_im, length(res), res)
    else
        ec = ccall((:perform_ifft, "libDataHandler.so"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
        	temp_re, temp_im, length(res), res)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_ifft ", ec), ec))
    end
    res
end

max_stars_repo_name: kusumikakd/Brainflow

max_issues_repo_name: kusumikakd/Brainflow

max_forks_repo_name: kusumikakd/Brainflow


text:
# ALL TESTS THAT CONCERN UNITCELLS

# include
using LatticePhysics
using Base.Test

# begin the Unitcell testblock
plaquette_testset = @testset "Plaquettes of lattices" begin

################################################################################
#
#    FINDING PLAQUETTES IN LATTICE OBJECTS
#
################################################################################

# begin the testset
@testset "Finding Plaquettes" begin

    # First, 2D
    @testset "Plaquettes in 2D (honeycomb)" begin

        # get a test unitcell
        unitcell = getUnitcellHoneycomb()
        # get a test lattice
        lattice = getLatticeByBondDistance(unitcell, 5)

        @testset "Finding plaquettes of one site" for testsite in [1,2,3]
            # does finding give the correct type
            @test typeof(getPlaquettesOfSite(lattice, testsite, 6)) == Array{Array{Int64, 1}, 1}
            # does finding without length give the correct type
            @test typeof(getPlaquettesOfSite(lattice, testsite)) == Array{Array{Int64, 1}, 1}
            # correct number of plaquettes found? 3 hexagons per site in honeycomb but nothing smaller or slightly bigger
            @test length(getPlaquettesOfSite(lattice, testsite)) == 3
            @test length(getPlaquettesOfSite(lattice, testsite, 6)) == 3
            @test sum([length(getPlaquettesOfSite(lattice, testsite, l)) for l in [3,4,5,7,8,9]]) == 0
        end;

        @testset "Finding plaquettes of complete lattice" begin
            # does finding give the correct type
            @test typeof(getPlaquettesOfLattice(lattice, 6)) == Array{Array{Int64, 1}, 1}
            # does finding without length give the correct type
            @test typeof(getPlaquettesOfLattice(lattice)) == Array{Array{Int64, 1}, 1}
            # correct number of plaquettes found? 3 hexagons per site in honeycomb but nothing smaller or slightly bigger
            @test length(getPlaquettesOfLattice(lattice)) > 0
            @test length(getPlaquettesOfLattice(lattice, 6)) > 0
            @test sum([length(getPlaquettesOfLattice(lattice, l)) for l in [3,4,5,7,8,9]]) == 0
        end;

    # end the 2D testset here
    end;

    # Then, 3D
    @testset "Plaquettes in 3D (hyperoctagon)" begin

        # get a test unitcell
        unitcell = getUnitcellHyperoctagon()
        # get a test lattice
        lattice = getLatticeByBondDistance(unitcell, 7)

        @testset "Finding plaquettes of one site" for testsite in [1,2,3]
            # does finding give the correct type
            @test typeof(getPlaquettesOfSite(lattice, testsite, 10)) == Array{Array{Int64, 1}, 1}
            # does finding without length give the correct type
            @test typeof(getPlaquettesOfSite(lattice, testsite)) == Array{Array{Int64, 1}, 1}
            # correct number of plaquettes found? 15 plaquettes of length 10 per site in hyperoctagon but nothing smaller or slightly bigger
            @test length(getPlaquettesOfSite(lattice, testsite)) == 15
            @test length(getPlaquettesOfSite(lattice, testsite, 10)) == 15
            @test sum([length(getPlaquettesOfSite(lattice, testsite, l)) for l in [3,4,5,6,7,8,9,11,12]]) == 0
        end;

        @testset "Finding plaquettes of complete lattice" begin
            # does finding give the correct type
            @test typeof(getPlaquettesOfLattice(lattice, 10)) == Array{Array{Int64, 1}, 1}
            # does finding without length give the correct type
            @test typeof(getPlaquettesOfLattice(lattice)) == Array{Array{Int64, 1}, 1}
            # correct number of plaquettes found? 3 hexagons per site in honeycomb but nothing smaller or slightly bigger
            @test length(getPlaquettesOfLattice(lattice)) > 0
            @test length(getPlaquettesOfLattice(lattice, 10)) > 0
            @test sum([length(getPlaquettesOfLattice(lattice, l)) for l in [3,4,5,6,7,8,9,11,12]]) == 0
        end;

    # end the 2D testset here
    end;

# end the testset here
end;

################################################################################
#
#    PRINTING INFORMATION ABOUT PLAQUETTES
#
################################################################################

# begin the testset
@testset "Printing Plaquette information" begin

    # get a test unitcell
    unitcell_2d = getUnitcellHoneycomb()
    unitcell_3d = getUnitcellHyperoctagon()
    # get a test lattice
    lattice_2d = getLatticeByBondDistance(unitcell_2d, 7)
    lattice_3d = getLatticeByBondDistance(unitcell_3d, 7)

    # ONE SITE
    @testset "Printing plaquette information for one site" for testsite in [1,2,3]
        @test printPlaquetteStatisticsOfSite(lattice_2d, testsite)==nothing
        @test printPlaquetteStatisticsOfSite(lattice_2d, testsite, detailed=true)==nothing
        @test printPlaquetteStatisticsOfSite(lattice_2d, testsite, max_length=10)==nothing
    end;

    # COMPLETE LATTICE
    @testset "Printing plaquette information for complete lattice" begin
        @test printPlaquetteStatisticsOfLattice(lattice_2d)==nothing
        @test printPlaquetteStatisticsOfLattice(lattice_2d, detailed=true)==nothing
        @test printPlaquetteStatisticsOfLattice(lattice_2d, max_length=10)==nothing
    end;

# end the test block here
end;

# end the plaquette test block here
end;