text:
type Parcel
  ID
  Index
  sDate
  Duration
  function Parcel()
    return new()
  end
end

"""
   Find gaps in time series (tsParm) structure and return duration of good data
     as well as an index into parcels of data which have contigously good data.

   Good data have values other than NaN.  Where there is a gap in the data
     larger than in.days4gap, a new parcel of data will be indicated.

   Any data which is NaN will be considered to form part of a gap in the time
     series such that where there are contiguous NaN's for a duration larger
     than days4gap will be marked as two separate parcels

 INPUT
  tsParm structure    -  from MetoceanTools package
  		      -  requires tsParm.sDate[1:nsDate] to be filled
		      -  requires tsParm.Data[1:nsDate,:] to be filled
  iDataColumn         -  column of tsParm.Data to use for data to be checked for gaps
  days4gap[1]         -  duration in days which is to be regarded as a gap in
                           the data

 OUTPUT:
   Parcel definitions as Parcel type
     ID 		- either 0 for a gap parcel or integer sequence number of parcel
     sDate		- start and end date of each parcel
     Index		- start and end index of each parcel where index is into tsParm.sDate
     Duration           - duration [dyas] of each parcel
   iParcel[1:nsDate]   -  either 0 for a gap larger than days4gap or a sequence number
                           giving the parcel number
   goodDurationDays    - the effective duration in days of the complete time series
                           (sum of the durations of each parcels)

 Created by Jason McConochie
 Revision 1 26Dec2016 by Jason McConochie
"""
function findGaps(in::tsParm,iDataColumn::Int64,days4gap::Float64)

   # A. Find the gaps greater than days4gap
   iParcel=zeros(Int64,size(in.sDate,1))
   ms2days=1/(1000 * 3600 * 24)
   inGap=false
   tParcel=1
   iSt=1
   iEd=1
   for i=1:length(in.sDate)
     isBad=isnan(in.Data[i,iDataColumn])
     if isBad & !inGap
       iSt=i
       inGap=true
     end
     if !isBad
       iEd=i-1
       dt=(Dates.value(in.sDate[iEd+1])-Dates.value(in.sDate[iSt]))*ms2days
       if inGap
         if 1.*dt>days4gap
           iParcel[iSt:iEd]=0
	   if iSt!=1
             tParcel+=1
	   end
	 else
           iParcel[iSt:iEd]=tParcel
	 end
       end
       iParcel[i]=tParcel
       inGap=false
     end
   end

   # B. Set duration of all parcels
   # B1. count number of parcels
   k=0
   lastParcel=-1
   for i=1:length(iParcel)-1
     if iParcel[i]!=lastParcel
        k=k+1
	lastParcel=iParcel[i]
     end
   end
   nParcels=k
   # B2. set start and end dates of parcels and parcel number associated
   p=Parcel()
   p.ID=Array{Int64,1}(nParcels)
   p.sDate=Array{DateTime,2}(nParcels,2)
   p.Index=Array{Int64,2}(nParcels,2)
   p.Duration=Array{Float64,1}(nParcels)
   lastParcel=-1
   k=0
   for i=1:length(iParcel)-1
     if iParcel[i]!=lastParcel
       k=k+1
       p.sDate[k,1]=in.sDate[i]
       p.Index[k,1]=i
       p.ID[k]=iParcel[i]
       lastParcel=iParcel[i]
     end
   end
   for i=1:nParcels-1
     p.sDate[i,2]=p.sDate[i+1,1]
     p.Index[i,2]=p.Index[i+1,1]
   end
   p.sDate[nParcels,2]=in.sDate[end]
   p.Index[nParcels,2]=length(in.sDate)

   goodDurationDays=0.
   for i=1:nParcels
     p.Duration[i]=ms2days*(Dates.value(p.sDate[i,2])-Dates.value(p.sDate[i,1]))
     if p.ID!=0
       goodDurationDays+=p.Duration[i]
     end
     #println("$(i) : $(p.ID[i]) : $(p.Index[i,1]),$(p.Index[i,2]) : $(p.sDate[i,1]),$(p.sDate[i,2]) : $(p.Duration[i])")
   end

   return p,goodDurationDays,iParcel

end

"""
   Function test_findGaps() to test findGaps function
"""
function test_findGaps()
   in=tsParm()
   in.sDate=[
      DateTime(1900,1,1),
      DateTime(1900,1,2),
      DateTime(1900,1,3),
      DateTime(1900,1,4),
      DateTime(1900,1,5),
      DateTime(1900,1,6),
      DateTime(1900,1,7),
      DateTime(1900,1,8),
      DateTime(1900,1,9),
      DateTime(1900,1,10),
      ]
   in.Data=zeros(10,1)
   in.Data[1:10,1]=[1,1,NaN,NaN,1,1,1,1,NaN,NaN]
   days4gap=1.5
   p,goodDurationDays,iParcel=findGaps(in,1,days4gap)

end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
module Unregistered
greet() = print("Hello World!")
end # module

ext: jl


text:
"""
    Boundary(body_pts,t_hats,n_hats)

Struct containing boundary information. The user must supply the location of the
body points (body_pts), and the tangent and normal vectors of the body at each
body point.

# ARGUMENTS
* `body_pts`   : Array of the form [[x1,y1,z1], [x2,y2,z2], ...] containing the
                 locations of the boundary nodes. Assumes a CCW orientation.
* `t_hats`     : Array of the same form containing the tangent unit vectors at
                 each ponit of the boundary. Same size as body_pts.
* `n_hats`     : Array of the same form containing the normal unit vectors at
                 each ponit of the boundary. Same size as body_pts.

# PROPERTIES
* `bodyPTS`    : Array of the form [[x1,y1,z1], [x2,y2,z2], ...] containing the
                 locations of the boundary nodes.
* `tHats`      : Array of the same form containing the tangent unit vectors at
                 each ponit of the boundary. Same size as body_pts.
* `nHats`      : Array of the same form containing the normal unit vectors at
                 each ponit of the boundary. Same size as body_pts.
* `NPTS_BODY`  : Number of points in the body.
"""
mutable struct Boundary

    bodyPTS::NUMB_MATRIX
    tHats::NUMB_MATRIX
    nHats::NUMB_MATRIX
    NPTS_BODY::Int

    function Boundary(body_pts::NUMB_MATRIX,
                      t_hats::NUMB_MATRIX,
                      n_hats::NUMB_MATRIX)
        new(body_pts,t_hats,n_hats,length(body_pts))
    end
end

"""
    `CalcRhoCoefs(self::Boundary)`

Takes a Boundary object and calculates the RBF coefficients for the eigenfunction
of the kernel decomposition.

# ARGUMENTS
* `self::Boundary`   : Boundary object
* `sigma = 0.2`      : Gaussian width

# RETURNS
* `etas::Array{Float64}` : Array containing the RBF coefficient values for rho
"""
function CalcRhoCoefs(self::Boundary; sigma = 0.2)

    # === Extract geometry from boundary ===
    X_body = [[point[1], point[2], 0.0] for point in self.bodyPTS]
    n_hats = [[n_hat[1], n_hat[2], 0.0] for n_hat in self.nHats]
    t_hats = [[t_hat[1], t_hat[2], 0.0] for t_hat in self.tHats]
    NPTS = self.NPTS_BODY

    # === Function definitions ===

    # i -> RBF index
    # j -> Integration index
    # k -> Evaluation index

    X_i(i) = X_body[i]                                  # Set_i of points on S
    X_j(j) = X_body[j]                                  # Set_j of points on S
    X_jp1(j) = (j != NPTS) ? X_body[j+1] : X_body[1]    # Point j + 1
    X_k(k) = X_body[k]                                  # Set_k of points on S

    R_kj(k,j) = X_k(k) - X_j(j)                         # X_k - X_j
    R_ji(j,i) = X_j(j) - X_i(i)                         # X_j - X_i
    R_ki(k,i) = X_k(k) - X_i(i)                         # X_k - X_i
    r_kj(k,j) = norm(R_kj(k,j))                         # | X_k - X_j |
    r_ji(j,i) = norm(R_ji(j,i))                         # | X_j - X_i |
    r_ki(k,i) = norm(R_ki(k,i))                         # | X_k - X_i |

    del_S(j) = norm(X_j(j) - X_jp1(j))                  # ΔS_j
    L = sum( [del_S(j) for j in 1:NPTS] )               # Length of surface

    # === Coefficient functions ===
    function Theta_ki(k,i)
        res = 0
        for j in 1:NPTS
            if j != k
                res = res + (1/(2*pi)) * (dot(R_kj(k,j), n_hats[j]) * RBF_gauss(r_ji(j,i), sigma=sigma) * del_S(j))/(r_kj(k,j)^2)
            end
        end
        return res
    end

    function phi_ki(k,i)
        return RBF_gauss(r_ki(k,i), sigma=sigma)
    end

    # === Build coefficient matrix ===
    A = zeros(NPTS,NPTS)
    for k in 1:NPTS
        for i in 1:NPTS
            A[k,i] = Theta_ki(k,i) - phi_ki(k,i)
        end
    end

    #println(A)

    # === Build RHS vector ===
    b = zeros(NPTS)

    # === Solve linear system ===
    etas = A\b
    return etas

end

"""
    `CalcRho(self::Boundary, etas, X_eval)`

Evaluates the RBF approxiamation function: rho(X) = sum(eta_i * phi(X - X_)).

# ARGUMENTS
* `self::Boundary`          : Boundary object
* `etas::Array{Float64}`    : Array containing RBF coefficients
* `X_eval::Array{Float64}`  : Point to evaluate function at
* `sigma = 0.2`             : Gaussian width

# RETURNS
* `rho::Float64`            : Evaluated function value
"""
function CalcRho(self::Boundary, etas::Array{T}, X_eval::Array{T}; sigma = 0.2) where {T<:Real}

    # === Unpack geometry ===
    XP = [[point[1], point[2], 0.0] for point in self.bodyPTS]
    NPTS = self.NPTS_BODY

    # === Summation over all points for gamma ===
    rho = sum( [etas[i] * RBF_gauss(norm(X_eval - XP[i]), sigma=sigma) for i in 1:NPTS ] )
    return rho

end

"""
    `Circle(R, O, N)`

Creates a Boundary object in the shape of a circle of radius R at origin O.

# ARGUMENTS
* `R::Float64`        : Radius of cirlce
* `O::Array{Float64}` : Origin of cirlce of form [X, Y, Z]
* `N::Int64`          : Number of points

# RETURNS
* `circle::Boundary`  : Boundary object of cirlce
"""
function Circle(R::Float64,
                O::Array{Float64},
                NPTS::Int64)

    del_theta = (2*pi)/(NPTS)              # Step in theta around body
    k_hat = [0.0, 0.0, 1.0]                 # Z direction unit vector
    X_coor = zeros(NPTS)                    # X-coordinates
    Y_coor = zeros(NPTS)                    # Y-coordinates
    X_tHat = zeros(NPTS)                    # X Tangent Vector
    Y_tHat = zeros(NPTS)                    # X Tangent Vector
    X_nHat = zeros(NPTS)                    # Y Normal Vector
    Y_nHat = zeros(NPTS)                    # Y Normal Vector
    for i=1:NPTS
        X_coor[i] = R*cos((i-1)*del_theta) + O[1]
        Y_coor[i] = R*sin((i-1)*del_theta) + O[2]

        X_tHat[i] = sin((i-1)*del_theta)
        Y_tHat[i] = -cos((i-1)*del_theta)

        X_nHat[i] = cos((i-1)*del_theta)
        Y_nHat[i] = sin((i-1)*del_theta)
    end
    body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
    t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
    n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]

    circle = VSB.Boundary(body_pts,t_hats,n_hats)
    return circle

end

"""
    `NACA4(numb, N, c)`

NACA Four-Digit Series Airfoil. Returns two lists of ordered pairs representing
a NACA Four-Digit Airfoil.

# ARGUMENTS
* `numb::String`     : Four digit series. Symmetric airfoil given by "00xx"
* `N::Int64`         : Number of points
* `c::Float64=1.0`   : Length of airfoil. Defaults to 1.0 for appropiate scaling
"""
function NACA4(numb::String,
               N::Int64,
               c::Float64=1.0)

    # Unpackage the digits
    m = parse(Int,numb[1])*.01
    p = parse(Int,numb[2])*.1
    tau = (parse(Int,numb[3])*10 + parse(Int,numb[4]))*.01

    Yt(x) = 5*tau*(0.2969*sqrt(x) - 0.1260*x - 0.3516*x^2 + 0.2843*x^3 - 0.1015*x^5)

    function Yc(x)
        if 0 <= x && x > p*c
            return (m/p^2) * (2*p*(x/c) - (x/c)^2)
        else
            return (m/(1 - p)^2) * ((1 - 2*p) + 2*p*(x/c) - (x/c)^2)
        end
    end

    function dYc_dX(x)
        if 0 <= x && x > p*c
            return (2*m/p^2) * (p - x/c)
        else
            return (2*m/(1 - p)^2) * (p - x/c)
        end
    end

    theta(x) = atan(dYc_dX(x))

    del_X = c/(N - 1)
    X = 0:del_X:c
    X_upper = [x - Yt(x)*sin(theta(x)) for x in X]
    X_lower = [x + Yt(x)*sin(theta(x)) for x in X]
    Y_upper = [Yc(x) + Yt(x)*cos(theta(x)) for x in X]
    Y_lower = [Yc(x) - Yt(x)*cos(theta(x)) for x in X]

    upper = [[X_upper[i], Y_upper[i]] for i in 1:N]
    lower = [[X_lower[i], Y_lower[i]] for i in 1:N]

    body_pts = [[0.0,0.0,0.0] for i in 1:2*N]
    for i in 0:N-1
        body_pts[i+1] = [X_upper[N-i], Y_upper[N-i], 0.0]
        body_pts[N+i+1] = [X_lower[i+1], Y_lower[i+1], 0.0]
    end

    t_hats = [[0.0,0.0,0.0] for i in 1:2*N]
    n_hats = t_hats
    for i in 1:2*N

        if i != 2*N
            R = body_pts[i] - body_pts[i+1]
        else
            R = body_pts[N] - body_pts[1]
        end

        r = norm(R)
        t_hats[i] = R/r
        n_hats[i] = cross(t_hats[i], [0,0,1])
    end

    NACA = Boundary(body_pts, t_hats, n_hats)

    return NACA

end


text:
struct Config
    apt::Vector{String}
    exclude_packages_from_sysimage::Vector{String}
    julia_cpu_target::String
    julia_version::Union{String, VersionNumber}
    make_sysimage::Bool
    no_test::Vector{String}
    packagecompiler_installation_command::String
    parent_image::String
    pkgs::Vector{Dict{Symbol,String}}
    precompile_execution_env_vars::Dict{String, String}
    simplecontainergenerator_installation_command::String
    tests_must_pass::Vector{String}
    wrapper_script_env_vars::Dict{String, String}
end

function Config(pkgs::AbstractVector{<:AbstractDict{<:Symbol,<:AbstractString}} =
                            Vector{Dict{Symbol,String}}(undef, 0);
                        #
                        additional_apt::AbstractVector{<:AbstractString} =
                            String[],
                        override_default_apt::AbstractVector{<:AbstractString} =
                            _default_apt(),
                        #
                        exclude_packages_from_sysimage =
                            String[],
                        julia_cpu_target =
                            _default_julia_cpu_target(),
                        julia_version::Union{AbstractString, VersionNumber} =
                            _default_julia_version(),
                        make_sysimage::Bool =
                            true,
                        no_test =
                            String[],
                        packagecompiler_installation_command::String =
                            _default_packagecompiler_installation_command(),
                        parent_image::String = _default_docker_parent_image(),
                        precompile_execution_env_vars =
                            _default_precompile_execution_env_vars(),
                        simplecontainergenerator_installation_command::String =
                            _default_simplecontainergenerator_installation_command(),
                        tests_must_pass =
                            String[],
                        wrapper_script_env_vars =
                            _default_wrapper_script_env_vars())
    apt = Vector{String}(undef, 0)
    append!(apt, override_default_apt)
    append!(apt, additional_apt)
    unique!(apt)
    config = Config(apt,
                    exclude_packages_from_sysimage,
                    julia_cpu_target,
                    julia_version,
                    make_sysimage,
                    no_test,
                    packagecompiler_installation_command,
                    parent_image,
                    pkgs,
                    precompile_execution_env_vars,
                    simplecontainergenerator_installation_command,
                    tests_must_pass,
                    wrapper_script_env_vars)
    return config
end

function Config(pkg::Union{AbstractDict{<:Symbol, <:AbstractString}, AbstractString, NamedTuple};
                        kwargs...)
    return Config([pkg]; kwargs...)
end

function Config(pkg_names::AbstractVector{<:AbstractString};
                        kwargs...)
    num_pkgs = length(pkg_names)
    pkgs = Vector{Dict{Symbol, String}}(undef, num_pkgs)
    for i = 1:num_pkgs
        pkgs[i] = Dict(:name => pkg_names[i])
    end
    return Config(pkgs; kwargs...)
end

function Config(pkg_tuples::AbstractVector{<:NamedTuple};
                        kwargs...)
    pkgs = Dict.(pairs.(pkg_tuples))
    return Config(pkgs; kwargs...)
end


text:
using PyCall

# PyCall情報表示
println(PyCall.pyversion)
println(PyCall.pyprogramname)
println(PyCall.libpython)

# python/selenium パッケージ import
const webdriver = pyimport("selenium.webdriver")
const WebDriverWait = pyimport("selenium.webdriver.support.ui").WebDriverWait
const EC = pyimport("selenium.webdriver.support.expected_conditions")
# python/chromedriver_binary パッケージ import
pyimport("chromedriver_binary")

# ChromeWebDriverで Google WEBサイト取得
driver = webdriver.Chrome()
driver.get("https://www.google.co.jp")
# ロード待ち
WebDriverWait(driver, 10).until(EC.presence_of_all_elements_located)
# 検索窓に "ChromeDriver" を入力し submit
search_box = driver.find_element_by_name("q")
search_box.send_keys("ChromeDriver")
search_box.submit()
# ロード待ち
WebDriverWait(driver, 10).until(EC.presence_of_all_elements_located)
# スクリーンショット保存
driver.save_screenshot("screenshot.png")
# ChromeWebDriver終了
driver.quit()

max_stars_repo_name: amenoyoya/wsgi-api-server

max_issues_repo_name: amenoyoya/wsgi-api-server

max_forks_repo_name: amenoyoya/wsgi-api-server


text:
"""
    solidBackground(
        raster::Array{Float64,3},
        color::Array{Float64,1}=[1.0,1.0,1.0]
        )::Array{Float64,3}

Take a 4-channel `raster` with alpha channel and convert it to 3-channel one by
adding a background `color`. Useful for saving non-transparent PNGs and JPEGs.
"""
function solidBackground(raster::Raster, color::Array{Float64,1} = [1.0, 1.0, 1.0])::Raster
    r = copy(raster)
    r[1, :, :] = r[1, :, :] .* r[4, :, :] + color[1] .* (1 .- r[4, :, :])
    r[2, :, :] = r[2, :, :] .* r[4, :, :] + color[2] .* (1 .- r[4, :, :])
    r[3, :, :] = r[3, :, :] .* r[4, :, :] + color[3] .* (1 .- r[4, :, :])
    r[1:3, :, :]
end

"""
    savePNG(filename::String, raster::Array{Float64,3})

Save the `raster` to a PNG `filename`. Transparency is choosen automatically
based on the number of raster channels.
"""
function savePNG(filename::String, raster::Raster)
    format = Images.RGBA
    if size(raster, 1) == 3
        format = Images.RGB
    elseif size(raster, 1) == 4
        format = Images.RGBA
    else
        throw(DomainError(size(raster, 1), "unsupported number of channels in raster"))
    end

    Images.save(FileIO.File(FileIO.format"PNG", filename), Images.colorview(format, raster))
end

"""
    saveJPEG(filename::String, raster::Array{Float64,3})

Save the `raster` to a JPEG `filename`. Only supports 3-channel rasters (see `solidBackground`).
"""
function saveJPEG(filename::String, raster::Raster)
    format = Images.RGB
    if size(raster, 1) == 3
        format = Images.RGB
    else
        throw(DomainError(size(raster, 1), "unsupported number of channels in raster"))
    end

    Images.save(
        FileIO.File(FileIO.format"JPEG", filename),
        Images.colorview(format, raster),
    )
end


text:
Lar = LinearAlgebraicRepresentation

function frag_face_channel(in_chan, out_chan, V, EV, FE, sp_idx)
    run_loop = true
    while run_loop
        sigma = take!(in_chan)
        if sigma != -1
            put!(out_chan, frag_face(V, EV, FE, sp_idx, sigma))
        else
            run_loop = false
        end
    end
end
function frag_face(V, EV, FE, sp_idx, sigma)
    vs_num = size(V, 1)

    sigmavs = (abs.(FE[sigma:sigma,:])*abs.(EV))[1,:].nzind
    sV = V[sigmavs, :]
    sEV = EV[FE[sigma, :].nzind, sigmavs]

    M = submanifold_mapping(sV)
    tV = ([V ones(vs_num)]*M)[:, 1:3]

    sV = tV[sigmavs, :]

    for i in sp_idx[sigma]
        tmpV, tmpEV = face_int(tV, EV, FE[i, :])

        sV, sEV = Lar.skel_merge(sV, sEV, tmpV, tmpEV)
    end

    sV = sV[:, 1:2]

    nV, nEV, nFE = planar_arrangement(sV, sEV, sparsevec(ones(Int8, length(sigmavs))))

    if nV == nothing
        return [], spzeros(Int8, 0,0), spzeros(Int8, 0,0)
    end

    nvsize = size(nV, 1)
    nV = [nV zeros(nvsize) ones(nvsize)]*inv(M)[:, 1:3]
    return nV, nEV, nFE
end

function merge_vertices(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp, err=1e-4)
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    facenum = size(FE, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    W = convert(Lar.Points, LinearAlgebra.transpose(V))
    kdtree = KDTree(W)

    global todelete = []

    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = Lar.inrange(kdtree, V[vi, :], err)

            newverts[nearvs] .= i

            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)

            i = i + 1
        end
    end

    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)

    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind

        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))

    end
    nedges = union(edges)
    nedges = filter(t->t[1]!=t[2], nedges)

    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()

    for ei in 1:nedgenum
    	begin
        	nEV[ei, collect(nedges[ei])] .= 1
        	nEV
        end
        etuple2idx[nedges[ei]] = ei
    end

    faces = [[
        map(x->newverts[x], FE[fi, ei] > 0 ? oedges[ei] : reverse(oedges[ei]))
        for ei in FE[fi, :].nzind
    ] for fi in 1:facenum]

    visited = []
    function filter_fn(face)

        verts = []
        map(e->verts = union(verts, collect(e)), face)
        verts = Set(verts)

        if !(verts in visited)
            push!(visited, verts)
            return true
        end
        return false
    end

    nfaces = filter(filter_fn, faces)

    nfacenum = length(nfaces)
    nFE = spzeros(Int8, nfacenum, size(nEV, 1))

    for fi in 1:nfacenum
        for edge in nfaces[fi]
            ei = etuple2idx[Tuple{Int, Int}(sort(collect(edge)))]
            nFE[fi, ei] = sign(edge[2] - edge[1])
        end
    end

    return Lar.Points(nV), nEV, nFE
end

"""
    spatial_arrangement(V::Points, EV::ChainOp, FE::ChainOp; [multiproc::Bool])

Compute the arrangement on the given cellular complex 2-skeleton in 3D.

A cellular complex is arranged when the intersection of every possible pair of cell
of the complex is empty and the union of all the cells is the whole Euclidean space.
The function returns the full arranged complex as a list of vertices V and a chain of borders EV, FE, CF.

## Additional arguments:
- `multiproc::Bool`: Runs the computation in parallel mode. Defaults to `false`.
"""
function spatial_arrangement(
V::Lar.Points,
EV::Lar.ChainOp,
FE::Lar.ChainOp, multiproc::Bool=false)

    fs_num = size(FE, 1)
    sp_idx = Lar.Arrangement.spatial_index(V, EV, FE)

    global rV = Lar.Points(undef, 0,3)
    global rEV = SparseArrays.spzeros(Int8,0,0)
    global rFE = SparseArrays.spzeros(Int8,0,0)

    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))

        @async begin
            for sigma in 1:fs_num
                put!(in_chan, sigma)
            end
            for p in Distributed.workers()
                put!(in_chan, -1)
            end
        end

        for p in Distributed.workers()
            @async Base.remote_do(
                frag_face_channel, p, in_chan, out_chan, V, EV, FE, sp_idx)
        end

        for sigma in 1:fs_num
            rV, rEV, rFE = Lar.skel_merge(rV, rEV, rFE, take!(out_chan)...)
        end

    else
        for sigma in 1:fs_num
            # print(sigma, "/", fs_num, "\r")
            nV, nEV, nFE = Lar.Arrangement.frag_face(
            	V, EV, FE, sp_idx, sigma)
            a,b,c = Lar.skel_merge(
            	rV, rEV, rFE, nV, nEV, nFE)
            global rV=a; global rEV=b; global rFE=c
        end

    end

    rV, rEV, rFE = merge_vertices(rV, rEV, rFE)

    rCF = minimal_3cycles(rV, rEV, rFE)

    return rV, rEV, rFE, rCF
end


text:
@testset "6.5 Hyperbolic secant" begin
    include("6.5.1 (c+d x)^m (a+b sech)^n.jl")
    include("6.5.2 (e x)^m (a+b sech(c+d x^n))^p.jl")
    include("6.5.3 Hyperbolic secant functions.jl")
    include("6.5.7 (d hyper)^m (a+b (c sech)^n)^p.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/.jl


text:
using Test
using DataWrangler

@testset "impute" begin

    x = 1:100;
    y = rand(100);

    @test isnothing(impute!(x,y))

    y = Array{Union{Missing,Float64}}(undef,100);
    y[:] = rand(100);

    yInt = Array{Union{Missing,Int64}}(undef,100);
    yInt[:] = rand(-1000:1000,100);

    # impute!
    ## loess
    x = sort(rand(100));
    y[rand(1:100,10)] .= missing;
    impute!(x,y)
    @test !any(ismissing,y)

    y[rand(1:100,10)] .= missing;
    impute!(y)
    @test !any(ismissing,y)

    yInt[rand(1:100,10)] .= missing;
    impute!(yInt)
    @test !any(ismissing,yInt)

    ## normal
    x = sort(rand(100));
    y[rand(1:100,10)] .= missing;
    impute!(x,y; method = "normal")
    @test !any(ismissing,y)

    y[rand(1:100,10)] .= missing;
    impute!(y; method = "normal")
    @test !any(ismissing,y)

    yInt[rand(1:100,10)] .= missing;
    impute!(yInt; method = "normal")
    @test !any(ismissing,yInt)

    ## uniform
    x = sort(rand(100));
    y[rand(1:100,10)] .= missing;
    impute!(x,y; method = "uniform")
    @test !any(ismissing,y)

    y[rand(1:100,10)] .= missing;
    impute!(y; method = "uniform")
    @test !any(ismissing,y)

    yInt[rand(1:100,10)] .= missing;
    impute!(yInt; method = "uniform")
    @test !any(ismissing,yInt)

    # impute
    ## loess
    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(x,y))

    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(y))

    yInt[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(yInt))

    ## normal
    x = sort(rand(100));
    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(x,y; method = "normal"))

    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(y; method = "normal"))

    yInt[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(yInt; method = "normal"))

    ## uniform
    x = sort(rand(100));
    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(x,y; method = "uniform"))

    y[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(y; method = "uniform"))

    yInt[rand(1:100,10)] .= missing;
    @test !any(ismissing,impute(yInt; method = "uniform"))

end


text:
@testset "filter_topk" begin

  @test_throws ArgumentError filter_topk(spzeros(5,5), -1)

  A = sparse([1.0 2.0 3.0;
       4.0 5.0 6.0;
       7.0 8.0 9.0])

  B = Array(filter_topk(A,1))
  @test B == [0.0 0.0 0.0;
              0.0 0.0 0.0;
              7.0 8.0 9.0]

  B = Array(filter_topk(A,0))
  @test B == [0.0 0.0 0.0;
              0.0 0.0 0.0;
              0.0 0.0 0.0]

  B = Array(filter_topk(sparse(
             [0.0 0.0 0.0;
              0.0 0.0 0.0]
  ), 1))
  @test B == [0.0 0.0 0.0;
              0.0 0.0 0.0]

  @test filter_topk!(spzeros(0,0),0) == spzeros(0,0)
  @test filter_topk!(spzeros(0,0),5) == spzeros(0,0)

  @test filter_topk!(spzeros(0,5),5) == spzeros(0,5)
  @test filter_topk!(spzeros(5,0),5) == spzeros(5,0)

  @test filter_topk!(spzeros(4,5),5) == spzeros(4,5)

  @test filter_topk!(spzeros(4,5),5) == spzeros(4,5)

  B = Array(filter_topk(sparse(
             [2.0 3.0 5.0;
              0.0 0.0 6.0;
              0.0 0.5 0.5]
  ), 3))
  @test B == [2.0 3.0 5.0;
   0.0 0.0 6.0;
   0.0 0.5 0.5]

end


text:
include("examples/doublependulum_disconnection.jl")

using ConstrainedDynamics: currentasknot!, getbody, geteqconstraint, getcomponent

dq = minimalCoordinates(mech)
@test true
setPosition!(mech,dq)
@test true

dF = ConstrainedDynamics.UnitDict(dq.keys,[[0.0] for i=1:2])
setForce!(mech,dF)
@test true

currentasknot!(mech)
@test true

body1 = getbody(mech, "body1")
body2 = getbody(mech, "body2")
eqc1 = geteqconstraint(mech, "jointb1")
eqc2 = geteqconstraint(mech, "joint12")
@test body1 === mech.bodies[1]
@test body2 === mech.bodies[2]
@test eqc1 === mech.eqconstraints[3]
@test eqc2 === mech.eqconstraints[4]

body1 = getcomponent(mech, "body1")
body2 = getcomponent(mech, "body2")
eqc1 = getcomponent(mech, "jointb1")
eqc2 = getcomponent(mech, "joint12")
@test body1 === mech.bodies[1]
@test body2 === mech.bodies[2]
@test eqc1 === mech.eqconstraints[3]
@test eqc2 === mech.eqconstraints[4]


text:
# Test external function calls / semantic attachments
@testset "external functions" begin

path = joinpath(dirname(pathof(PDDL)), "..", "test", "functions")

# Load domain and define external functions
domain = load_domain(joinpath(path, "domain.pddl"))
throw_range(v, θ) = v^2*sin(2*θ*π/180) / 9.81
domain.funcdefs[:range] = throw_range
throw_height(v, θ, x) = tan(θ*π/180)*x - 9.81*x^2 / (2*v^2 * cos(θ*π/180)^2)
domain.funcdefs[:height] = throw_height

# Load problem
problem = load_problem(joinpath(path, "problem.pddl"))
state = initstate(domain, problem)

# Check that evaluation with external functions works correctly
@test domain[state => pddl"(range 20 45)"] == throw_range(20, 45)
@test domain[state => pddl"(height 20 45 10)"] == throw_height(20, 45, 10)

# Execute plan
state = initstate(domain, problem)
state = execute(domain, state, pddl"(pick ball1)")
state = execute(domain, state, pddl"(throw ball1 85)")
state = execute(domain, state, pddl"(pick ball2)")
state = execute(domain, state, pddl"(throw ball2 75)")

# Check if goal is satisfied
@test domain[state => pddl"(< (loc ball1) 10)"]
@test domain[state => pddl"(> (loc ball2) 15)"]
@test satisfy(domain, state, problem.goal) == true

end # external functions


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

@testset "Meshes" begin

  # Load 2D .msh file from disk and create from vertices, elements
  mesh1 = Mesh(joinpath(meshdir_2D, "circle25_32.msh"));
  vtx, ele, elt = data(mesh1)
  mesh2 = Mesh(vtx,ele,elt)

  for m in [mesh1 mesh2]
    @test nodecount(m) == 3511
    @test elemcount(m) == 6840
    @test dimensions(m) == 2
    @test boundingbox(m) == [-25.0 -25.0; 25.0 25.0]
    @test fullsize(m) == 1963.0959349692218
    @test maxnodes(m) == 3
    vtx, ele = surface(m)
    @test vtx[10,:] == [19.7003,15.3915]
    @test ele[1:5] == [90,135,89,164,88]
  end

  # Load 3D .msh file from disk and create from vertices, elements
  mesh1 = Mesh(joinpath(meshdir_3D, "cyl3.msh"));
  vtx, ele, elt = data(mesh1)
  mesh2 = Mesh(vtx,ele,elt)

  for m in [mesh1 mesh2]
    @test nodecount(m) == 27084
    @test elemcount(m) == 141702
    @test dimensions(m) == 3
    @test boundingbox(m) == [-25.0 -25.0 -25.0; 25.0 25.0 25.0]
    @test fullsize(m) == 98056.71142861326
    @test maxnodes(m) == 4
    vtx, ele = surface(m)
    @test vtx[10,:] == [-1.29504,6.49008,25.0]
    @test ele[1:5] == [ 2791,1429,3014,2552,124]
  end

end


text:
##############################################################################################################################
# Conversion and Promotion rules
# Small convention: Convert from type T to type S
##############################################################################################################################

"""
	mat(x::T where T<:AbstractVector, dim=LearnBase.ObsDim.Constant{1})

Function that transforms a vector into a row or column matrix.
"""
mat(x::T where T<:AbstractVector) = mat(getobs(x), LearnBase.ObsDim.Constant{1}())
mat(x::T where T<:AbstractVector, dim::S where S<:LearnBase.ObsDim.Constant{1}) = reshape(getobs(x), size(x,1), 1)
mat(x::T where T<:AbstractVector, dim::S where S<:LearnBase.ObsDim.Constant{2}) = reshape(getobs(x), 1, size(x,1))
mat(x::T where T<:AbstractMatrix, dim::S where S<:LearnBase.ObsDimension) = getobs(x)
mat(x::T where T<:Vector) = mat(x, LearnBase.ObsDim.Constant{1}())
mat(x::T where T<:Vector, dim::S where S<:LearnBase.ObsDim.Constant{1}) = reshape(x, size(x,1), 1)
mat(x::T where T<:Vector, dim::S where S<:LearnBase.ObsDim.Constant{2}) = reshape(x, 1, size(x,1))
mat(x::T where T<:Matrix, dim::S where S<:LearnBase.ObsDimension) = x

# Convert from Vectors to Matrices
convert(::Type{S}, x::T where T<:AbstractVector) where S<:AbstractMatrix = mat(x)

# Convert from CellData to AbstractArray
convert(::Type{S}, x::T where T<:DataCell) where S<:AbstractArray = convert(S, getx!(x))

# Convert to DataCell
convert(::Type{S}, x::T where T<:DataElement) where S<:DataCell = DataCell([x], nothing, nothing, 0, "DataCell (convert)", eval(oinfoglobal))
convert(::Type{S}, x::T where T<:AbstractArray) where S<:DataCell = DataCell(x, nothing, nothing, 0, "DataCell (convert)", eval(oinfoglobal))

# Convert to FunctionCell
convert(::Type{S}, x::T where T<:DataElement) where S<:FunctionCell = FunctionCell(nothing, nothing, (args...)->x, (), (), 0, "Constant output=$x", eval(oinfoglobal))
convert(::Type{S}, f::T where T<:Function) where S<:FunctionCell = FunctionCell(f, (), string(f))

# Convert to AbstractCell (default conversions to Cells-like objects)
convert(::Type{S}, x::T where T<:DataElement) where S<:AbstractCell = DataCell([x], nothing, nothing, 0, "DataCell (convert)", eval(oinfoglobal))
convert(::Type{S}, x::T where T<:AbstractArray) where S<:AbstractCell = DataCell(x, "DataCell (convert)")
convert(::Type{S}, f::T where T<:Function) where S<:AbstractCell = FunctionCell(f, (), string(f))
convert(::Type{S}, x::T where T<:Tuple) where S<:AbstractCell = S(x...)


text:
module MixedDistributions

import Base: maximum, minimum, rand
import Statistics
using Statistics: mean, var, std
using Random: AbstractRNG

import Distributions
const Dst = Distributions

import RecipesBase

export MixedDistribution

include("mixed_distribution.jl")
include("plot.jl")

end # module


text:
struct Trans2Covariate{S<:AbstractFloat, T<:AbstractFloat, U<:AbstractFloat} <: AbstractCovariate{U}
    name::String
    basecovariate1::AbstractCovariate{S}
    basecovariate2::AbstractCovariate{T}
    transform::Function
end

Base.length(var::Trans2Covariate) = length(var.basecovariate1)

function Trans2Covariate(name::String, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T},
                        transform::Function) where {S<:AbstractFloat} where {T<:AbstractFloat}
    U = promote_type(S, T)
    Trans2Covariate{S, T, U}(name, basecovariate1, basecovariate2, transform)
end

function Trans2Covariate(::Type{U}, name::String, basecovariate1::AbstractCovariate{S}, basecovariate2::AbstractCovariate{T},
                         transform::Function) where {U<:AbstractFloat} where {S<:AbstractFloat} where {T<:AbstractFloat}
    Trans2Covariate{S, T, U}(name, basecovariate1, basecovariate2, transform)
end

function slice(covariate::Trans2Covariate{S, T, U}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractFloat} where {T<:AbstractFloat} where {U<:AbstractFloat}
    basecov1 = covariate.basecovariate1
    basecov2 = covariate.basecovariate2
    f = covariate.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip2(slice(basecov1, fromobs, toobs, slicelength), slice(basecov2, fromobs, toobs, slicelength))
    mapslice2(f, slices, slicelength, U)
end

function Base.map(covariate::Trans2Covariate, dataframe::AbstractDataFrame)
    Trans2Covariate(covariate.name, map(covariate.basecovariate1, dataframe), map(covariate.basecovariate1, dataframe),
                    covariate.transform)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
#=

    Inference.jl

    A Julia module code to perform ABC SMC with CDF/correlation matching.

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#
module Inference

    using AbstractMCMC
    using AdvancedMH
    using DataFrames
    using DataFramesMeta
    using Distributions
    using Interpolations
    using LinearAlgebra
    using MCMCChains
    using .Threads
    using StatsBase
    using Statistics
    using StatsPlots
    using Printf

    export abc_smc, distances, locations, abc_mcmc, acceptance_rate, ess, logistic, logit
    export makediscrepancy, makediscrepancy_univariate
    export plot_posteriors, plot_chains
    export pooled_data, pooled_std, pooled_cov, pooled_cor
    export Particle

    include("abc.jl")
    include("particles.jl")
    include("discrepancies.jl")
    include("plotting.jl")
    include("pool_data.jl")

end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using Test, Pkg

ENV["JULIA_LOG_LEVEL"] = "WARN"

for submodule in ["Utilities/ParametersType",
                  "Utilities/PlanetParameters",
                  "Utilities/RootSolvers",
                  "Utilities/MoistThermodynamics",
                  "Atmos/Parameterizations/SurfaceFluxes",
                  "Atmos/Parameterizations/TurbulenceConvection",
                  "Mesh",
                  "DGmethods",
                  "Atmos/Dycore"
                  ]

  println("Testing $submodule")
  include(joinpath("../src",submodule,"test/runtests.jl"))
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
module TestPipeline

using Random
using Test
using AutoMLPipeline
using AutoMLPipeline.Pipelines
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.DecisionTreeLearners
using AutoMLPipeline.Utils
using AutoMLPipeline.FeatureSelectors

global const data = getiris()
global const features = data[:,1:4]
global const X=data[:,1:5]
const Y=data[:,5] |> Vector
X[!,5]= X[!,5] .|> string

global const ohe = OneHotEncoder()
global const pca = SKPreprocessor("PCA")
global const ica = SKPreprocessor("FastICA")
global const fa = SKPreprocessor("FactorAnalysis")
global const disc = CatNumDiscriminator()
global const catf = CatFeatureSelector()
global const numf = NumFeatureSelector()
global const rf = RandomForest()
global const ada = Adaboost()
global const pt = PrunedTree()

function test_pipeline()
  # test initialization of types
  ohe = OneHotEncoder()
  linear1 = Pipeline(Dict(:name=>"lp",:machines => [ohe]))
  linear2 = Pipeline(Dict(:name=>"lp",:machines => [ohe]))
  combo1 = ComboPipeline(Dict(:machines=>[ohe,ohe]))
  combo2 = ComboPipeline(Dict(:machines=>[linear1,linear2]))
  linear1 = Pipeline([ohe])
  linear2 = Pipeline([ohe])
  combo1 = ComboPipeline([ohe,ohe])
  combo2 = ComboPipeline([linear1,linear2])
  # test fit/transform workflow
  fit!(combo1,X)
  res1=transform!(combo1,X)
  res2=fit_transform!(combo1,X)
  @test (res1 .== res2) |> Matrix |> sum == 2100
  fit!(combo2,X)
  res3=transform!(combo2,X)
  res4=fit_transform!(combo2,X)
  @test (res3 .== res4) |> Matrix |> sum == 2100
  pcombo1 = @pipeline ohe + ohe
  pres1 = fit_transform!(pcombo1,X)
  @test (pres1 .== res1) |> Matrix |> sum == 2100
end
@testset "Pipelines" begin
  Random.seed!(123)
  test_pipeline()
end

function test_sympipeline()
  pcombo5 = @pipeline :((numf |> pca) + (numf |> ica) |> (ada * rf * pt))
  @test crossvalidate(pcombo5,X,Y,"accuracy_score",10,false).mean >= 0.90
  expr = :((numf |> pca) + (numf |> ica) |> (ada * rf * pt))
  processexpr!(expr.args)
  @test crossvalidate(eval(expr),X,Y,"accuracy_score",10,false).mean >= 0.90
  expr = :((numf |> pca) + (numf |> ica) |> (ada * rf * pt))
  pcombo6 = sympipeline(expr) |> eval
  @test crossvalidate(pcombo6,X,Y,"accuracy_score",10,false).mean >= 0.90
  pcombo7 = (@pipelinez expr) |> eval
  @test crossvalidate(pcombo7,X,Y,"accuracy_score",10,false).mean >= 0.90
end
@testset "Symbolic Pipeline: Global Scope" begin
  Random.seed!(123)
  test_sympipeline()
end

function test_pipeline()
  ohe = OneHotEncoder()
  pca = SKPreprocessor("PCA")
  ica = SKPreprocessor("FastICA")
  fa = SKPreprocessor("FactorAnalysis")
  disc = CatNumDiscriminator()
  catf = CatFeatureSelector()
  numf = NumFeatureSelector()
  rf = RandomForest()
  ada = Adaboost()
  pt = PrunedTree()
  # test symbolic pipeline expression
  pcombo2 = @pipeline (pca |> ica) + ica + pca
  @test fit_transform!(pcombo2,features) |> Matrix |> size |> collect |> sum == 162
  pcombo2 = @pipeline pca |> ica |> fa
  @test fit_transform!(pcombo2,features) |> Matrix |> size |> collect |> sum == 154
  pcombo3 = @pipeline disc |> ((catf + numf) + (numf |> pca) + (numf |> ica) + (catf|>ohe)) |> rf
  (fit_transform!(pcombo3,X,Y)  .== Y) |> sum == 150
  pcombo4 = @pipeline (numf |> pca) + (numf |> ica) |> (ada * rf * pt)
  @test crossvalidate(pcombo4,X,Y,"accuracy_score",10,false).mean >= 0.90
end
@testset "Symbolic Pipeline: Local Scope" begin
  Random.seed!(123)
  test_pipeline()
end

end


text:
import Base: ==, convert, Dict

const zklmTuple = NamedTuple{(:z, :k, :l, :m), Tuple{Int16, IntS, IntS, IntS}}

"""
`AList` : datastructure to help compute the A_zklm density projections

* `i2zklm` : list of all admissible `(z,k,l,m)` tuples
* `zklm2i` : dictionary with (z,k,l,m) keys to compute  the map `(z,k,l,m) -> i`
* `firstz` : `firstz[iz]` stores the first index in the A_zklm array for with
             z = zi. This can be used to iterate over all A entries for which
             z = zi. (they are sorted by z first)
"""
struct AList
   i2zklm::Vector{zklmTuple}
   zklm2i::Dict{zklmTuple, IntS}
   firstz::Vector{IntS}
end

# --------------(de-)serialisation----------------------------------------
Dict(alist::AList) = Dict("__id__" => "SHIPs_AList",
                           "i2zklm" => zklm2vec.(alist.i2zklm))
AList(D::Dict) = AList(vec2zklm.(D["i2zklm"]))
zklm2vec(t) = [t.z, t.k, t.l, t.m]
vec2zklm(v) = (z=Int16(v[1]), k=IntS(v[2]), l=IntS(v[3]), m=IntS(v[4]))
==(al1::AList, al2::AList) = (al1.i2zklm == al2.i2zklm)
# ------------------------------------------------------------------------

Base.length(alist::AList) = length(alist.i2zklm)
Base.getindex(alist::AList, i::Integer) = alist.i2zklm[i]
Base.getindex(alist::AList, zklm::zklmTuple) = alist.zklm2i[zklm]

alloc_A(alist::AList, T=Float64) = zeros(Complex{T}, length(alist))

function AList(zklmlist::AbstractVector{zklmTuple})
   # sort the tuples - by z, then k, then l, then m
   i2zklm = sort(zklmlist)
   # create the inverse mapping
   zklm2i = Dict{zklmTuple, IntS}()
   for i = 1:length(i2zklm)
      zklm2i[i2zklm[i]] = IntS(i)
   end
   # find the first index for each z
   zmax = maximum( a.z for a in i2zklm )
   firstz = [ findfirst([a.z == iz for a in i2zklm])
              for iz = 1:zmax ]
   return AList( i2zklm, zklm2i, [firstz; length(i2zklm)+1] )
end

"""
This fills the A-array stored in tmp with the A_zklm density projections in
the order specified by AList. It also evaluates the radial and angular basis
functions along the way.
"""
function precompute_A!(A, tmp, alist, Rs, Zs, ship)
   fill!(A, 0)
   for (R, Z) in zip(Rs, Zs)
      # evaluate the r-basis and the R̂-basis for the current neighbour at R
      eval_basis!(tmp.J, tmp.tmpJ, ship.J, norm(R))
      eval_basis!(tmp.Y, tmp.tmpY, ship.SH, R)
      # add the contributions to the A_zklm
      iz = z2i(ship, Z)
      for i = alist.firstz[iz]:(alist.firstz[iz+1]-1)
         zklm = alist[i]
         A[i] += tmp.J[zklm.k+1] * tmp.Y[index_y(zklm.l, zklm.m)]
      end
   end
   return A
end

function precompute_dA!(A, dA, tmp, alist, Rs, Zs, ship)
   fill!(A, 0)
   fill!(dA, zero(eltype(dA)))
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      # precompute the derivatives of the Jacobi polynomials and Ylms
      eval_basis_d!(tmp.J, tmp.dJ, tmp.tmpJ, ship.J, norm(R))
      eval_basis_d!(tmp.Y, tmp.dY, tmp.tmpY, ship.SH, R)
      # deduce the A and dA values
      iz = z2i(ship, Z)
      R̂ = R / norm(R)
      for i = alist.firstz[iz]:(alist.firstz[iz+1]-1)
         zklm = alist[i]
         ik = zklm.k+1; iy = index_y(zklm.l, zklm.m)
         A[i] += tmp.J[ik] * tmp.Y[iy]
         # and into dA # grad_phi_Rj(R, iR, zklm, tmp)
         ∇ϕ_zklm = tmp.dJ[ik] * tmp.Y[iy] * R̂ + tmp.J[ik] * tmp.dY[iy]
         dA[iR, i] = ∇ϕ_zklm
      end
   end
   return dA
end

# ---------

"""
`AAList` : datastructure to help compute the A_𝐳𝐤𝐥𝐦 = ∏ A_zklm

* `i2Aidx` : indices in AList of the zklms to avoid the Dict lookup
* `len`    : len[i] is the number of relevant entries of i2zklm[i,:]
             i.e. the body-order of this entry
* `zklm2i` : dictionary of all (z,k,l,m) tuples to compute  the
             map `(zz,kk,ll,mm) -> i`
* `firstz` : `firstz[iz]` stores the first index in the A_zklm array for with
             z = zi. This can be used to iterate over all A entries for which
             z = zi. (they are sorted by z first)
"""
struct AAList
   i2Aidx::Matrix{IntS}    # where in A can we find these
   len::Vector{IntS}       # body-order
   zklm2i::Dict{Any, IntS} # inverse mapping
end

# --------------(de-)serialisation----------------------------------------

==
function Dict(aalist::AAList)
   ZKLM_list = Vector{Any}(undef, length(aalist))
   for (zzkkllmm, i) in aalist.zklm2i
      ZKLM_list[i] = zzkkllmm
   end
   return Dict("__id__" => "SHIPs_AAList", "ZKLM_list" => ZKLM_list)
end
zzkkllmm2vec(zzkkllmm) = Vector.([zzkkllmm...])
vec2zzkkllmm(v) = (SVector(Int16.(v[1])...),
                   SVector(IntS.(v[2])...),
                   SVector(IntS.(v[3])...),
                   SVector(IntS.(v[4])...))
AAList(D::Dict, alist) = AAList(vec2zzkkllmm.(D["ZKLM_list"]), alist)
==(aal1::AAList, aal2::AAList) = (aal1.i2Aidx == aal2.i2Aidx)
# ------------------------------------------------------------------------

bodyorder(aalist::AAList) = maximum(aalist.len) + 1

Base.length(aalist::AAList) = length(aalist.len)

Base.getindex(aalist::AAList, t::Tuple) = aalist.zklm2i[t]

alloc_AA(aalist::AAList, T = Float64) = zeros(Complex{T}, length(aalist))

"""
`AAList(ZKLM_list, alist)` : standard AAList constructor,
* `ZKLM_list` : a collection of (izz=?, kk=?, ll=?, mm=?) tuples
* `alist` a compatible `AList`
"""
function AAList(ZKLM_list, alist)
   BO = maximum(ν -> length(ν[1]), ZKLM_list)  # body-order -> size of iAidx

   # create arrays to construct AAList
   iAidx = IntS[]
   len = IntS[]
   zklm2i = Dict{Any, IntS}()

   idx = 0
   for (izz, kk, ll, mm) in ZKLM_list
      # store in the index of the current row in the reverse map
      idx += 1
      zklm2i[(izz, kk, ll, mm)] = idx
      # store the body-order of the current ∏A function
      push!(len, length(ll))

      # fill the row of the i2Aidx matrix
      for α = 1:length(ll)
         zklm = (z=izz[α], k=kk[α], l=ll[α], m=IntS(mm[α]))
         iA = alist[zklm]
         push!(iAidx, iA)
      end
      # fill up the iAidx vector with zeros up to the body-order
      # this will create 0 entries in the matrix after reshaping
      for α = (length(ll)+1):BO
         push!(iAidx, 0)
      end
   end
   return AAList( reshape(iAidx, (BO, idx))', len, zklm2i )
end

function precompute_AA!(AA, A, aalist) # tmp, ship::SHIPBasis{T}, iz0) where {T}
   fill!(AA, 1)
   for i = 1:length(aalist)
      for α = 1:aalist.len[i]
         iA = aalist.i2Aidx[i, α]
         AA[i] *= A[iA]
      end
   end
   return nothing
end

# --------------------------------------------------------
# this section of functions is for computing
#  ∂∏A / ∂R_j

# function grad_phi_Rj(Rj, j, zklm, tmp)
#    ik = zklm.k + 1
#    iy = index_y(zklm.l, zklm.m)
#    return ( tmp.dJJ[j, ik] *  tmp.YY[j, iy] * (Rj/norm(Rj))
#            + tmp.JJ[j, ik] * tmp.dYY[j, iy] )
# end

function grad_AA_Ab(iAA, b, alist, aalist, A)
   g = one(eltype(A))
   for a = 1:aalist.len[iAA]
      if a != b
         iA = aalist.i2Aidx[iAA, a]
         g *= A[iA]
      end
   end
   return g
end

function grad_AAi_Rj(iAA, j, Rj::JVec{T}, izj::Integer,
                     alist, aalist, A, AA, dA, tmp) where {T}
   g = zero(JVec{Complex{T}})
   for b = 1:aalist.len[iAA] # body-order
      # A_{n_b} = A[iA]
      iA = aalist.i2Aidx[iAA, b]
      # daa_dab = ∂(∏A_{n_a}) / ∂A_{n_b}
      daa_dab = grad_AA_Ab(iAA, b, alist, aalist, A)
      zklm = alist[iA] # (zklm corresponding to A_{n_b})
      ∇ϕ_zklm = dA[j, iA] # grad_phi_Rj(Rj, j, zklm, tmp)
      g += daa_dab * ∇ϕ_zklm
   end
   return g
end

function grad_AA_Rj!(tmp, ship, j, Rs, Zs, iz0) where {T}
   for iAA = 1:length(ship.aalists[iz0])
      # g = ∂(∏_a A_a) / ∂Rj
      tmp.dAAj[iz0][iAA] = grad_AAi_Rj(iAA, j, Rs[j], z2i(ship, Zs[j]),
                                       ship.alists[iz0], ship.aalists[iz0],
                                       tmp.A[iz0], tmp.AA[iz0], tmp.dA[iz0],
                                       tmp)
   end
   return tmp.dAAj[iz0]
end

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

"""
convert the "old" `(NuZ, ZKL)` format into the simpler (zz, kk, ll, mm)
format, and at the same time extract the one-particle basis (z, k, l, m)
"""
function alists_from_bgrps(bgrps::Tuple)
   NZ = length(bgrps)
   zzkkllmm_list = [ Tuple[] for _=1:NZ ]
   zklm_set = [ Set() for _=1:NZ ]
   for iz0 = 1:NZ
      for (izz, kk, ll) in bgrps[iz0], mm in _mrange(ll)
         push!(zzkkllmm_list[iz0], (izz, kk, ll, IntS.(mm)))
         for α = 1:length(ll)
            zklm = (z=izz[α], k=kk[α], l=ll[α], m=IntS(mm[α]))
            push!(zklm_set[iz0], zklm)
         end
      end
   end

   alist =  ntuple(iz0 -> AList([ zklm for zklm in collect(zklm_set[iz0]) ]), NZ)
   aalist = ntuple(iz0 -> AAList(zzkkllmm_list[iz0], alist[iz0]), NZ)
   return alist, aalist
end

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

using SparseArrays: SparseMatrixCSC

function _my_mul!(C::AbstractVector, A::SparseMatrixCSC, B::AbstractVector)
   A.n == length(B) || throw(DimensionMismatch())
   A.m == length(C) || throw(DimensionMismatch())
   nzv = A.nzval
   rv = A.rowval
   fill!(C, zero(eltype(C)))
   @inbounds for col = 1:A.n
      b = B[col]
      for j = A.colptr[col]:(A.colptr[col + 1] - 1)
         C[rv[j]] += nzv[j] * b
      end
   end
   return C
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
module ScottishTaxBenefitModel

#
# A microsimulation tax benefit model of Scotland.
#
# This is the `parent module` that just imports (and re-exports) all its child modules.
#

include("Utils.jl" )
export Utils

include( "Randoms.jl")
export Randoms

include( "Monitor.jl")
export Monitor

include( "TimeSeriesUtils.jl" )
export TimeSeriesUtils

include("Definitions.jl" )
export Definitions

include( "RunSettings.jl" )
export RunSettings

include("STBIncomes.jl")
export STBIncomes

include( "EquivalenceScales.jl")
export EquivalenceScales

include("GeneralTaxComponents.jl" )
export GeneralTaxComponents

include("Uprating.jl" )
export Uprating

include("ModelHousehold.jl" )
export ModelHousehold

include("Weighting.jl" )
export Weighting

include( "HistoricBenefits.jl")
export HistoricBenefits

include( "Pensions.jl")
export Pensions

include("HouseholdFromFrame.jl" )
export HouseholdFromFrame

include("ExampleHouseholdGetter.jl" )
export ExampleHouseholdGetter

include("FRSHouseholdGetter.jl" )
export FRSHouseholdGetter

include("STBParameters.jl" )
export STBParameters

include( "BenefitGenerosity.jl")
export BenefitGenerosity

include( "Intermediate.jl")
export Intermediate

include( "Results.jl")
export Results

include( "NonMeansTestedBenefits.jl")
export NonMeansTestedBenefits

include("IncomeTaxCalculations.jl" )
export IncomeTaxCalculations

include("NationalInsuranceCalculations.jl" )
export NationalInsuranceCalculations

include("LocalLevelCalculations.jl" )
export LocalLevelCalculations

include("LegacyMeansTestedBenefits.jl" )
export LegacyMeansTestedBenefits

include( "UniversalCredit.jl")
export UniversalCredit

include( "ScottishBenefits.jl")
export ScottishBenefits

include( "BenefitCap.jl" )
export BenefitCap

include( "UCTransition.jl")
export UCTransition

include( "UBI.jl")
export UBI

include("SingleHouseholdCalculations.jl" )
export SingleHouseholdCalculations

include( "ExampleHelpers.jl")
export ExampleHelpers

include( "BCCalcs.jl")
export BCCalcs

include( "SimplePovertyCounts.jl")
export SimplePovertyCounts

# pro. tem include("WebModelLibs.jl" )
# export WebModelLibs

include( "STBOutput.jl")
export STBOutput

include( "Runner.jl" )
export Runner

include( "TheEqualiser.jl" )
export TheEqualiser

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using KentFinancialEconomics
using Test

@testset "KentFinancialEconomics.jl" begin
    # Write your tests here.
end

max_stars_repo_name: alfredjmduncan/KentFinancialEconomics

max_issues_repo_name: alfredjmduncan/KentFinancialEconomics

max_forks_repo_name: alfredjmduncan/KentFinancialEconomics


text:
const examples_dir = joinpath(dirname(dirname(@__FILE__)), "examples")

const Examples = Dict(
    "HydroValleys" => [
            "case3.jl",
            # "case3_soc.jl",
            "case3deterministic.jl",
            "case3deterministic_nowater.jl",
            "case3deterministic_overwater.jl"
        ]

)

@testset "Examples" begin
    for (key, examples) in Examples
        @testset "$(key)" begin
            for example in examples
                @testset "$example" begin
                    println("Running $(key):$(example)")
                    include(joinpath(examples_dir, key, example))
                end
            end
        end
    end
end


text:
function rebalance_every_n_days(t::Weights, interval::Integer)::Weights
    counter = 0
    last_value = 0.0

    function mapping(x)
        counter += 1
        if counter <= interval
            return NaN
        else
            counter = 0
            last_value = x
            return x
        end
    end
    return TimeArray(timestamp(t), map(mapping, values(t))) |> dropnan
end

function rebalance_once(t::Weights)::Weights
    return t[1]
end


text:
import LinearAlgebra: ldiv!

"""
    tridiagonalsolve!(x, b, u, d, l)

Solve the tridiagonal system `A*x = b`, where `A` is tridiagonal and defined
by the upper diagonal `u`, diagonal `d`, and lower diagonal `l`

Following https://julialang.org/blog/2016/02/iteration for the multidimensional implementation.
"""
function ldiv!(x, A::Tridiagonal, b)
    N = length(x)
    N == length(b) || throw("x and b must have the same length.")
    _tridiagonalsolve!(x, b, A.dl, A.d, A.du, N)
    return nothing
end

function tridiagonalsolve!(x, d, a, b, c)
    length(x) == length(d) == (length(a)+1) == length(b) == (length(c)+1) || (
        throw(ArgumentError("The size of x, d, a, b, and c do not match.")))
    _tridiagonalsolve!(x, d, a, b, c, length(x))
end

function _tridiagonalsolve!(x, d, a, b, c, N)
    for i = 2:N
        @inbounds begin
            b[i] -= a[i-1] * c[i-1] / b[i-1]
            d[i] -= a[i-1] * d[i-1] / b[i-1]
        end
    end

    @inbounds x[N] = d[N] / b[N]

    for i = N-1:-1:1
        # Solve for x by bulk substitution
        @inbounds x[i] = (d[i] - c[i]*x[i+1] ) / b[i]
    end

    return nothing
end


text:
# From JuliaLang/Combinatorics

immutable Combinations{T}
    a::T
    t::Int
end

start(c::Combinations) = [1:c.t;]
function next(c::Combinations, s)
    comb = [c.a[si] for si in s]
    if c.t == 0
        # special case to generate 1 result for t==0
        return (comb,[length(c.a)+2])
    end
    s = copy(s)
    for i = length(s):-1:1
        s[i] += 1
        if s[i] > (length(c.a) - (length(s)-i))
            continue
        end
        for j = i+1:endof(s)
            s[j] = s[j-1]+1
        end
        break
    end
    (comb,s)
end
done(c::Combinations, s) = !isempty(s) && s[1] > length(c.a)-c.t+1

length(c::Combinations) = binomial(length(c.a),c.t)

eltype{T}(::Type{Combinations{T}}) = Vector{eltype(T)}

"Generate all combinations of `n` elements from an indexable object. Because the number of combinations can be very large, this function returns an iterator object. Use `collect(combinations(array,n))` to get an array of all combinations.
"
function combinations(a, t::Integer)
    if t < 0
        # generate 0 combinations for negative argument
        t = length(a)+1
    end
    Combinations(a, t)
end

"""
generate combinations of all orders, chaining of order iterators is eager,
but sequence at each order is lazy
"""
combinations(a) = chain([combinations(a,k) for k=1:length(a)]...)


text:
abstract type SetConstraint  <: JuMP.AbstractConstraint end

## Loading  ##

function need_variablify(set::Union{Polyhedra.Rep, Polyhedra.RepElement,
                                    Sets.AbstractSet})
    return false
end
function variablify(set::Union{Polyhedra.Rep, Polyhedra.RepElement,
                               Sets.AbstractSet})
    return set
end
need_variablify(v::SetVariableRef) = true
variablify(v::SetVariableRef) = v.variable

function load(model::JuMP.Model, constraint::SetConstraint)
    @assert need_variablify(constraint)
    return JuMP.add_constraint(model, variablify(constraint))
end

## Adding ##

struct SetShape <: JuMP.AbstractShape end
struct ConstraintIndex
    value::Int
end

function JuMP.add_constraint(model::JuMP.Model, constraint::SetConstraint,
                             name::String="")
    d = data(model)
    if d.state != Modeling
        error("Constraint of type $(typeof(constraint)) not supported yet")
    end
    @assert need_variablify(constraint)
    index = ConstraintIndex(d.last_index += 1)
    d.constraints[index] = constraint
    d.names[index] = name
    return JuMP.ConstraintRef(model, index, SetShape())
end

### SetConstraintRef ###

const SetConstraintRef{M} = JuMP.ConstraintRef{M, ConstraintIndex, SetShape}
function JuMP.name(cref::SetConstraintRef)
    return data(cref.model).names[cref.index]
end
function JuMP.constraint_object(cref::SetConstraintRef)
    return data(cref.model).constraints[cref.index]
end
function transformed_constraint(cref::SetConstraintRef)
    return data(cref.model).transformed_constraints[cref.index]
end
function JuMP.dual(cref::SetConstraintRef)
    return JuMP.dual(transformed_constraint(cref))
end
function SumOfSquares.moment_matrix(cref::SetConstraintRef)
    return SumOfSquares.moment_matrix(transformed_constraint(cref))
end

include("membership.jl")
include("inclusion.jl")
include("tangent.jl")


text:
using LPIP, Test

include("lpip_interface_test.jl")
# include("moi_test.jl")


text:
# Get paths.
markdown_path = joinpath(@__DIR__, "markdown")

function polish_latex(path::String)
    txt = open(f -> read(f, String), path)
    open(path, "w+") do f
        txt = replace(txt, raw"$$" => raw"\$\$")

        beg_align = r"\\\$\\\$\n?\\begin{align}"
        end_align = r"\\end{align}\n?\\\$\\\$?"

        txt = replace(txt, beg_align => "\$\$\n\\begin{align}")
        txt = replace(txt, end_align => "\\end{align}\n\$\$")

        # Make inline math work.
        inline = r"(?<=[^\\])\$(?!\$)"
        txt = replace(txt, inline => raw"$$")

        # If there's more than three dollar signs together, use only two.
        extra_dollars = r"\${3,}"
        txt = replace(txt, extra_dollars => raw"$$")

        write(f, txt)
    end
end

function add_yaml(path)
    lines = readlines(path, keep=true)

    title = ""
    permalink = "/:collection/:name/"

    # Find the ID sections.
    for i = 1:length(lines)
        line = lines[i]

        if startswith(line, "# ")
            title = replace(line, "# " => "")
            title = replace(title, "\n" => "")

            lines[i] = ""
            break
        end
    end

    yaml_lines = ["---\n", "title: $title\n", "permalink: $permalink\n", "---\n"]

    # Write lines back.
    open(path, "w+") do f
        # Prepend YAML.
        for line in yaml_lines
            write(f, line)
        end

        # Add remaining content.
        for line in lines
            write(f, replace(line, "![](figure" => "![](/tutorials/figure"))
        end
    end
end

function handle_file(fn)
    polish_latex(fn)
    add_yaml(fn)
end

max_stars_repo_name: SebastianCallh/TuringTutorials

max_issues_repo_name: SebastianCallh/TuringTutorials

max_forks_repo_name: SebastianCallh/TuringTutorials


text:
struct NNKolmogorov{C,O,S} <: NeuralNetDiffEqAlgorithm
    chain::C
    opt::O
    sdealg::S
end
NNKolmogorov(chain  ; opt=Flux.ADAM(0.1) , sdealg = EM()) = NNKolmogorov(chain , opt , sdealg)

function DiffEqBase.solve(
    prob::KolmogorovPDEProblem,
    alg::NeuralNetDiffEqAlgorithm;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    save_everystep = false,
    dt,
    kwargs...
    )

    tspan = prob.tspan
    phi = prob.phi
    xspan = prob.xspan
    sigma = prob.sigma
    μ = prob.μ
    d = prob.d
    ts = tspan[1]:dt:tspan[2]
    xs = xspan[1]:0.001:xspan[2]
    N = size(ts)
    T = tspan[2]
    #hidden layer
    chain  = alg.chain
    opt    = alg.opt
    sdealg = alg.sdealg
    ps     = Flux.params(chain)
    xi     = rand(xs , d , N[1])
    #Finding Solution to the SDE having initial condition xi. Y = Phi(S(X , T))
    sdeproblem = SDEProblem(μ,sigma,xi,tspan)
    sol = solve(sdeproblem, sdealg ,dt=0.01 , save_everystep=false , kwargs...)
    x_sde = sol[end]
    y = phi(x_sde)
    data   = Iterators.repeated((xi , y), maxiters)

    #MSE Loss Function
    loss(x , y) =Flux.mse(chain(x), y)

    cb = function ()
        l = loss(xi, y)
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end

    Flux.train!(loss, ps, data, opt; cb = cb)
    xi , chain(xi)
 end #solve


text:
K80(κ::Float64, safe::Bool=true) = K80rel(κ, safe)

K80(α::Float64, β::Float64, safe::Bool=true) = K80abs(α, β, safe)

function K80(θ_vec::A,
              safe::Bool=true) where A <: AbstractArray
  if length(θ_vec) == 1
    return K80rel(θ_vec[1], safe)
  elseif length(θ_vec) == 2
    return K80abs(θ_vec[1], θ_vec[2], safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function K80rel(θ_vec::A,
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 1
      error("Incorrect parameter vector length")
    end
  end
  return K80rel(θ_vec[1], safe)
end

function K80abs(θ_vec::A,
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 2
      error("Incorrect parameter vector length")
    end
  end
  return K80abs(θ_vec[1], θ_vec[2], safe)
end


text:
module InstallCmd

using JSON
using ProgressMeter
using Crayons.Box
using Pkg.BinaryPlatforms
using Comonicon.Tools: prompt, cmd_error
using SHA: sha256
using Downloads: download, Downloader
using ..IonBase: dot_ion, ca_roots
using ..Options

# NOTE:
# when the user tries to install a certain julia version, it is very possible
# that the user wants to use it, so we create a symlink to julia by default
function install(target::String, version::String, activate::Bool=true, yes::Bool=false, cache::Bool=false)
    if target == "julia"
        install_julia(version, activate, yes, cache)
    else
        cmd_error("unkown target: $target")
    end
end

function install_julia(version_string::String="stable", activate::Bool=true, yes::Bool=false, cache::Bool=false)
    toml = Options.read()

    if Options.find_julia_bin(version_string, toml) !== nothing
        # TODO: update Comonicon to change the default answer to no
        if !prompt("$(CYAN_FG("julia-", version_string)) exists, do you want to install again?", yes)
            return
        end
    end

    version, file = download_julia(version_string, cache)

    # make sure the path is available
    ispath(toml.julia.install_dir) || mkpath(toml.julia.install_dir)

    julia_bin = if Sys.islinux()
        install_julia_linux(version, file, toml)
    elseif Sys.isapple()
        install_julia_mac(version, file, toml)
    elseif Sys.iswindows()
        install_julia_win(version, file, toml)
    else
        cmd_error("unsupported system")
    end

    # successful installed update meta
    if version_string == "stable"
        toml.julia.stable = version
    end

    bin = dot_ion("bin")

    if activate
        toml.julia.active = julia_bin
        create_symlink(julia_bin, "julia", bin)
    end

    if version isa VersionNumber
        julia_minor = "julia-$(version.major).$(version.minor)"
        for each in readdir(bin)
            if startswith(each, julia_minor)
                @info "clean up old patch version: $(CYAN_FG(each))"
                rm(each; force=true)
            end
        end

        # allocate a list so we can delete the keys
        @info "updating $(CYAN_FG("ion.toml"))"
        for key in collect(keys(toml.julia.versions))
            if (key.major == version.major) && (key.minor == version.minor)
                delete!(toml.julia.versions, key)
            end
        end

        create_symlink(julia_bin, "julia-$(version.major).$(version.minor)", bin)
        toml.julia.versions[version] = julia_bin
    else # nightly
        # we create two link: julia-nightly and julia-latest
        create_symlink(julia_bin, "julia-latest", bin)
        toml.julia.nightly = julia_bin
    end

    create_symlink(julia_bin, "julia-$version", bin)
    Options.dump(toml)

    if !cache
        # NOTE: installation won't be very frequently
        # so we just rm tarballs intermediately by default
        rm(file; force=true)
    end
    return julia_bin
end

function _triplet()
    if Sys.islinux()
        BinaryPlatforms.triplet(Linux(Sys.ARCH))
    elseif Sys.isapple()
        BinaryPlatforms.triplet(MacOS(Sys.ARCH))
    elseif Sys.iswindows()
        BinaryPlatforms.triplet(Windows(Sys.ARCH))
    else
        cmd_error("unsupported system")
    end
end

function download_julia_version_json()
    io = IOBuffer()
    downloader = Downloader()
    @static if Sys.isapple()
        downloader.ca_roots = ca_roots()
    end
    download("https://julialang-s3.julialang.org/bin/versions.json", io; downloader=downloader)
    raw = JSON.parse(String(take!(io)))
    version_info = Dict{VersionNumber, Any}()
    for (k, v) in raw
        version_info[VersionNumber(k)] = v
    end
    return version_info
end

function find_julia_stable(version_info::Dict{VersionNumber, Any}=download_julia_version_json())
    stable = v"0.0.0"
    for (v, info) in version_info
        if stable < v && info["stable"]
            stable = v
        end
    end

    return query_julia_downloads(stable, version_info)
end

function query_julia_downloads(version_string::String, version_info = download_julia_version_json())
    version = Options.find_version(version_string, keys(version_info))
    version === nothing && cmd_error("cannot find version: $version_string")
    return query_julia_downloads(version::VersionNumber, version_info)
end

function query_julia_downloads(version::VersionNumber, version_info = download_julia_version_json())
    for file in version_info[version]["files"]
        if _triplet() == file["triplet"]
            return version, file["sha256"], file["url"]
        end
    end
    cmd_error("cannot find julia binary for this platform")
end

function find_julia_installer_info(version_string::String)
    if version_string == "stable"
        return find_julia_stable()
    elseif version_string == "latest" || version_string == "nightly"
        return find_nightly()
    else
        return query_julia_downloads(version_string)
    end
end

function find_nightly()
    if Sys.ARCH == :x86_64
        arch = 64
        platform = "x64"
    elseif Sys.ARCH == :i686
        arch = 32
        platform = "x86"
    elseif Sys.ARCH == :aarch64
        arch = 64
        platform = "aarch64"
    end

    if Sys.isapple()
        sys = "mac"
        installer = string("mac", arch, ".dmg")
    elseif Sys.islinux()
        sys = "linux"
        if platform == "aarch64"
            installer = "linuxaarch64.tar.gz"
        else
            installer = string("linux", arch, ".tar.gz")
        end
    elseif Sys.iswindows()
        sys = "winnt"
        installer = string("win", arch, ".exe")
    end

    url = "https://julialangnightlies-s3.julialang.org/bin/$sys/$platform/julia-latest-$installer"
    return "nightly", nothing, url
end

"""
returns `VersionNumber` or "nightly"
"""
function download_julia(version::String="stable", cache::Bool=false)
    version, sha, url = find_julia_installer_info(version)

    if cache
        ispath(dot_ion("cache")) || mkpath(dot_ion("cache"))
        file = joinpath(dot_ion("cache"), basename(url))
    else
        file = joinpath(tempdir(), basename(url))
    end

    @info "downloading from $url"
    p = Progress(0, 0.1, "julia-$version  ")

    function progress(total, now)
        p.n = total
        if now <= total
            update!(p, now)
        end
        return
    end

    if cache && isfile(file) && validate_julia_download(sha, file)
        @info "julia installer is cached, using existing installer"
        return version, file
    end

    downloader = Downloader()
    @static if Sys.isapple()
        downloader.ca_roots = ca_roots()
    end

    file = download(url, file; progress=progress, downloader=downloader)
    if validate_julia_download(sha, file)
        @info "checksum: $(CYAN_FG("true"))"
    else
        cmd_error("download is not complete, please try again")
    end
    return version, file
end

# do not validate nightly downloads
validate_julia_download(::Nothing, ::String) = true

function validate_julia_download(sha::String, file::String)
    SHA = open(file) do io
        bytes2hex(sha256(io))
    end

    return SHA == sha
end

function withdmg(f, file::String, mountpoint::String = splitext(basename(file))[1]; maxtry::Int = 5)
    @info "mounting $mountpoint"
    cmd = ignorestatus(`hdiutil attach $file -quiet -mount required -mountpoint $mountpoint`)
    code = 1
    for _ in 1:maxtry
        code = run(cmd).exitcode
        if code == 0
            break
        end
    end
    code == 0 || cmd_error("unable to mount $file")

    ret = f(mountpoint)
    @info "umounting $mountpoint"
    run(`umount $mountpoint`)
    return ret
end

function withtar(f, file::String, mountpoint::String = mktempdir())
    @info "extracting tar $file"
    run(`tar -zxf $file -C $mountpoint`)
    ret = f(mountpoint)
    @info "removing $mountpoint"
    rm(mountpoint; force=true, recursive=true)
    return ret
end

function install_julia_linux(version, file::String, ion::Options.Ion = Options.read())
    postfix = version == "nightly" ? version : "$(version.major).$(version.minor)"
    install_path = joinpath(ion.julia.install_dir, "julia-$postfix")

    if ispath(install_path)
        @info "removing old path: $install_path"
        rm(install_path; force=true, recursive=true)
    end

    withtar(file) do mountpoint
        src = mountpoint
        @info "copying new julia files from $src"
        cp(src, install_path; force=true, follow_symlinks=true)
    end

    return joinpath(install_path, first(readdir(install_path)), "bin", "julia")
end

function install_julia_mac(version, file::String, ion::Options.Ion = Options.read())
    return withdmg(file) do mountpoint
        dirs = readdir(mountpoint)
        idx = findfirst(startswith("Julia"), dirs)
        idx !== nothing || cmd_error("julia installer does not contain julia binary")
        julia_dir = dirs[idx]
        install_path = joinpath(ion.julia.install_dir, julia_dir)
        julia_bin = "$install_path/Contents/Resources/julia/bin/julia"

        if ispath(install_path)
            @info "removing old path: $(CYAN_FG(install_path))"
            rm(install_path; force=true, recursive=true)
        end

        src = joinpath(mountpoint, julia_dir)
        @info "copying new julia files from $(CYAN_FG(src))"
        cp(src, install_path; force=true, follow_symlinks=true)
        return julia_bin
    end
end

function install_julia_win(version, file::String, ion::Options.Ion = Options.read())
    cmd_error("not implemented")
end

# NOTE: we don't use joinpath(homedir(), .local), since we prefer to allow users to
# remove all effects that ion creates by simply `rm -rf .ion`
# this also makes testing easier, since we can define where .ion is
function create_symlink(src::String, name::String, bin=dot_ion("bin"))
    if !ispath(bin)
        mkpath(bin)
    end

    link = joinpath(bin, name)
    if islink(link) || ispath(link)
        @info "$(CYAN_FG(link)) exists, removing..."
        rm(link; force=true)
    end
    @info "creating symlink: $(CYAN_FG(link))"
    symlink(src, link)
    return
end

end # InstallCmd

# TODO: use traits to implement this so we can make it extensible

"""
install a given target. Currently only supports "julia",

# Arguments

- `target`: install given target, currently only supports `julia`.
- `version`: version you want to install, default is `stable`.

# Flags

- `-y, --yes`: always choose yes.
- `--cache`: cache downloads in `.ion/cache`.
- `--no-activate`: do not activate the `julia` installation.
"""
@cast function install(target::String, version::String="stable"; no_activate::Bool=false, yes::Bool=false, cache::Bool=false)
    try
        InstallCmd.install(target, version, !no_activate, yes, cache)
    catch e
        if (e isa InterruptException) || (e isa TaskFailedException && e.task.result isa InterruptException)
            print("canceled by user")
            return
        else
            rethrow(e)
        end
    end
    return
end


text:
abstract type AbstractBinomialModel end

struct BinomialModel{T1,T2} <: AbstractBinomialModel
    N::T1
    p::T2
    q::T2
    σ::T2
    τ::T2
end

struct BinomialGridModel{T1, T2, T3, T4, T5} <: AbstractBinomialModel
    Nind::T1
    pind::T1
    qind::T1
    σind::T1
    τind::T1
    Nrng::T2
    prng::T3
    qrng::T3
    σrng::T3
    τrng::T3
    N::T4
    p::T5
    q::T5
    σ::T5
    τ::T5
end

function BinomialGridModel(
    Nind, pind, qind, σind, τind,
    Nrng, prng, qrng, σrng, τrng
)
    N = Nrng[Nind]
    p = prng[pind]
    q = qrng[qind]
    σ = σrng[σind]
    τ = τrng[τind]

    return BinomialGridModel(
        Nind, pind, qind, σind, τind,
        Nrng, prng, qrng, σrng, τrng,
        N,    p,    q,    σ,    τ
    )
end

function BinomialGridModel(m_out::Int, my_Nrng, my_prng, my_qrng, my_σrng, my_τrng)
    Nrng = CuArray(Int.(my_Nrng))
    prng = CuArray(Float32.(my_prng))
    qrng = CuArray(Float32.(my_qrng))
    σrng = CuArray(Float32.(my_σrng))
    τrng = CuArray(Float32.(my_τrng))

    Nind = cu(rand(1:length(Nrng), m_out))
    pind = cu(rand(1:length(prng), m_out))
    qind = cu(rand(1:length(qrng), m_out))
    σind = cu(rand(1:length(σrng), m_out))
    τind = cu(rand(1:length(τrng), m_out))

    return BinomialGridModel(
                Nind, pind, qind, σind, τind,
                Nrng, prng, qrng, σrng, τrng
            )
end

function refresh!(model::BinomialGridModel)
    model.N .= model.Nrng[model.Nind]
    model.p .= model.prng[model.pind]
    model.q .= model.qrng[model.qind]
    model.σ .= model.σrng[model.σind]
    model.τ .= model.τrng[model.τind]
    return model
end

# special outer constructor to convert a BinomialGridModel into a BinomialModel
function BinomialModel(model::BinomialGridModel)
    refresh!(model)
    return BinomialModel(model.N, model.p, model.q, model.σ, model.τ)
end

function BinomialModel(nmax::Int, m_out::Int, device = :gpu)
    if device == :gpu
        N = CuArray(rand(1:nmax, m_out))
        p = CUDA.rand(m_out)
        q = CUDA.rand(m_out)
        σ = CUDA.rand(m_out)
        τ = CUDA.rand(m_out)
    elseif device == :cpu
        N = rand(1:nmax, m_out)
        p = rand(m_out)
        q = rand(m_out)
        σ = rand(m_out)
        τ = rand(m_out)
    end
    return BinomialModel(N, p, q, σ, τ)
end

function BinomialModel(m_out::Int, my_Nrng, my_prng, my_qrng, my_σrng, my_τrng)
    gridmodel = BinomialGridModel(m_out, my_Nrng, my_prng, my_qrng, my_σrng, my_τrng)
    return BinomialModel(gridmodel)
end

function ScalarBinomialModel(nmax, device = :cpu)
    return BinomialModel(nmax, 1, device)
end

function ScalarBinomialModel(N::Int, p, q, σ, τ, device = :cpu)
    if device == :cpu
        Ns = N .* ones(Int, 1)
        ps = p .* ones(1)
        qs = q .* ones(1)
        σs = σ .* ones(1)
        τs = τ .* ones(1)
    elseif device == :gpu
        Ns = N .* CUDA.ones(Int, 1)
        ps = Float32(p) .* CUDA.ones(1)
        qs = Float32(q) .* CUDA.ones(1)
        σs = Float32(σ) .* CUDA.ones(1)
        τs = Float32(τ) .* CUDA.ones(1)
    end
    return BinomialModel(Ns, ps, qs, σs, τs)
end

struct BinomialState{T}
    n::T
    k::T
end

function BinomialState(nmax::Int, m_out::Int, m_in::Int, device = :gpu)
    if device == :gpu
        n = CuArray(rand(1:nmax, m_out, m_in))
        k = CUDA.zeros(Int, m_out, m_in)
    elseif device == :cpu
        n = rand(1:nmax, m_out, m_in)
        k = zeros(Int, m_out, m_in)
    end
    return BinomialState(n, k)
end

function ScalarBinomialState(nmax::Int, device = :cpu)
    return BinomialState(nmax, 1, 1, device)
end

function ScalarBinomialState(n::Int, k::Int, device = :cpu)
    if device == :cpu
        ns = n .* ones(Int, 1, 1)
        ks = k .* ones(Int, 1, 1)
    elseif device == :gpu
        ns = n .* CUDA.ones(Int, 1, 1)
        ks = k .* CUDA.ones(Int, 1, 1)
    end
    return BinomialState(ns, ks)
end

struct BinomialObservation{T1, T2}
    EPSP::T1
    dt::T2
end


text:
# packages required for test
using Test
using LinearAlgebra
using Random
using ForwardDiff
using Zygote

# packages tested
using KFEstimate

@testset "KFEstimate.jl" begin
    # Write your tests here.
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "Test KF" begin
        include(joinpath(testdir, "test_kf.jl"))
    end
end


text:
"""
`get_hierarchy`
Return a vector indicating how many consecutive falses or trues
preceded an event
"""

function nextcount(count::T, rewarded) where {T <: Number}
    if rewarded
        count > 0 ? count + 1 : one(T)
    else
        count < 0 ? count - 1 : -one(T)
    end
end
signedcount(v) = accumulate(nextcount, v;init=0.0)
get_hierarchy(v) = lag(signedcount(v), default = NaN)

Base.findall(v::AbstractVector{Union{Bool,Missing}}) = findall(collect(skipmissing(v)))
Base.findlast(v::AbstractVector{Union{Bool,Missing}}) = findlast(collect(skipmissing(v)))
Base.findprev(v::AbstractVector{Union{Bool,Missing}},idx) = findprev(collect(skipmissing(v)),idx)

function exp_calendar(df::AbstractDataFrame)
    ExpCalendar = Dict(d => n for (n,d) in enumerate(sort(union(df[:,:Day]))))
    [get(ExpCalendar,x,Date(2000,12,31)) for x in df[:,:Day]]
end

function exp_calendar(df::AbstractDataFrame,Phase::Symbol)
    PhaseCalendar = Dict()
    combine(groupby(df,Phase)) do dd
        for (n,d) in enumerate(sort(union(dd[:,:Day])))
            PhaseCalendar[d] = n
        end
    end
    [get(PhaseCalendar,x,Date(2000,12,31)) for x in df[:,:Day]]
end


text:
matrix = [1 2 3 4
          3 4 5 6
          4 4 4 6
          3 3 3 3] #Define the matrix
f(x,y) = 2x+y #Create an inline function
@time for i=1:4,j=1:4
  matrix[i,j] += f(i,j) #update matrix[i,j] = matrix[i,j] + f(i,j)
end #end statement is required for control statements like for
Pkg.add("ASCIIPlots") #Only use the first time!
using ASCIIPlots #Load the package
scatterplot(1:4,matrix[end,:],sym='*') #Plot the last row
scatterplot(1:4,matrix[2,:],sym='^') #Plot the second row
package = "ASCIIPlots" #Define a string
println("$package is a super cool package") #Print a string

ran = 1:5 #Define a range
vec = collect(ran) #Turn the range into a vector
vec2=Vector{Float64}(2) #Define a size 2 vector of 64-bit floats
vec2[1]=2; vec2[2]=3   # vec2 = [2,3], ; allows multiple commands
push!(vec2,5) # vec2 = [2,5]
mat = Array{Int64,2}(5,5) # Create an empty 5x5 matrix
mat2 = ones(5,5) #Create a 5x5 matrix of ones
Q,R = qr(mat2) #QR decompose mat2 (returns a tuple) into Q and R
tup = ran,vec,vec2 #Define a tuple of many different objects
tup[2] #This returns vec
dict = Dict("a" => 1, "b" => 2, "c" => 3) #Define a dictionary
dict["a"] #Returns 1

function testFunction(x,y;z=0) #z is an optional argument
  if z==0
    return x+y,x*y #Return a tuple
  else
  return x*y*z,x+y+z #Return a different tuple
  #whitespace is optional
  end #End if statement
end #End function definition
a,b = testFunction(2,2) #Returns 4,4
a2,b2 = testFunction(2,3,z=3) #Returns 18,8


text:
using OAuth, Test

@test typeof(oauth_timestamp()) <: AbstractString

@test parse(Int, oauth_timestamp()) > 1422235471

@test typeof(oauth_nonce(15)) <: AbstractString

@test length(oauth_nonce(15)) == 15
@test length(oauth_nonce(20)) == 20
@test length(oauth_nonce(25)) == 25
@test length(oauth_nonce(30)) == 30
@test length(oauth_nonce(32)) == 32

test_oauth1a = OAuth1a("consumer_key", "9djdj82h48djs9d2", "token", "kkk9d7dh3k39sjv7")

@test oauth_sign_hmac_sha1("randy", test_oauth1a) == "TrJTbr3xSI6bhJUggXo7tJicLOg="

@test oauth_signing_key(test_oauth1a) == "9djdj82h48djs9d2&kkk9d7dh3k39sjv7"

@test oauth_signature_base_string("POST", "http://example.com", "?julia=fast&lang=elegant") ==
        "POST&http%3A%2F%2Fexample.com&%3Fjulia%3Dfast%26lang%3Delegant"

@test oauth_percent_encode_keys!(Dict("badkey!"   => "value", "goodkey" => "value")) == Dict("badkey%21" => "value", "goodkey" => "value")

@test oauth_serialize_url_parameters(Dict("language" => "julia", "result" => "awesome")) == "language=julia&result=awesome"

@test encodeURI!(Dict("iv" => 10, "s" => "value!")) == Dict("iv" => 10, "s" => "value%21")

@test oauth_body_hash_encode("Hello, World!") == "CgqfKmdylCVXq1NV12r0Qvj2XgE="

@test oauth_body_hash_data("Hello, World!") == "oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="

test_file = joinpath(dirname(@__FILE__), "auth_body_hash_file.txt")
@test oauth_body_hash_file(test_file) == "oauth_body_hash=CgqfKmdylCVXq1NV12r0Qvj2XgE="

#TODO
#oauth_header()

#TODO
#oauth_request_resource()


text:
#=
jim.jl
jiffy image display
2019-02-23 Jeff Fessler, University of Michigan
=#

export jim

using Plots: heatmap, plot, plot!, annotate!
using MosaicViews: mosaicview
using FFTViews: FFTView

# global default key/values
jim_def = Dict([
	:aspect_ratio => :equal,
	:clim => nothing,
	:color => :grays,
	:line3plot => true, # lines around sub image for 3d mosaic?
	:line3type => (:yellow),
	:ncol => 0,
	:padval => nothing,
	:mosaic_npad => 1,
	:tickdigit => 1,
	:title => "",
	:xlabel => "",
	:ylabel => "",
	:fft0 => false,
	:yflip => nothing, # defer to minimum value of y - see default below
	:yreverse => nothing, # defer to whether y is non-ascending
	:abswarn => isinteractive(), # warn when taking abs of complex images?
	])

minfloor = x -> floor(minimum(x), digits=jim_def[:tickdigit])
maxceil = x -> ceil(maximum(x), digits=jim_def[:tickdigit])

"""
    nothing_else(x, y)
return `y` if `x` is nothing, else return `x`
"""
function nothing_else(x, y)
	return x == nothing ? y : x
end

"""
    jim(z, ...)

A jiffy image display of `x` using `heatmap`

in
- `z` image, can be 2D or higher, if higher then it uses `mosaicviews`

option
- `aspect_ratio`; default `:equal`
- `clim`; default `(minimum(z),maximum(z))`
- `color` (colormap, e.g. `:hsv`); default `:grays`
- `ncol` for mosaicview for 3D and higher arrays; default `0` does auto select
- `padval` padding value for mosaic view; default `minimum(z)`
- `line3plot` lines around sub image for 3d mosaic; default `false`
- `line3type` line type around sub image for 3d mosaic; default `(:yellow)`
- `mosaic_npad` # of pixel padding for mosaic view; default `1`
- `fft0` if true use FFTView to display; default `false`
- `title`; default `""`
- `xlabel`; default `""`
- `ylabel`; default `""`
- `yflip`; default `true` if `minimum(y) >= 0`
- `yreverse`; default `true` if `y[1] > y[end]`
- `x` values for x axis; default `1:size(z,1)`
- `y` values for y axis; default `1:size(z,2)`
- `xtick`; default `[minimum(x),maximum(x)]`
- `ytick`; default `[minimum(y),maximum(y)]`

out
- returns plot handle, type `Plots.Plot`

2019-02-23 Jeff Fessler, University of Michigan
"""
function jim(z::AbstractArray{<:Real} ;
		aspect_ratio = jim_def[:aspect_ratio],
		clim = nothing_else(jim_def[:clim], (minimum(z), maximum(z))),
		color = jim_def[:color],
		line3plot = jim_def[:line3plot],
		line3type = jim_def[:line3type],
		ncol::Int = jim_def[:ncol],
		padval = nothing_else(jim_def[:padval], minimum(z)),
		mosaic_npad::Int = jim_def[:mosaic_npad],
		title::AbstractString = jim_def[:title],
		xlabel::AbstractString = jim_def[:xlabel],
		ylabel::AbstractString = jim_def[:ylabel],
		fft0::Bool = jim_def[:fft0],
		x = fft0 ? Int.(-size(z,1)/2:size(z,1)/2-1) : (1:size(z,1)),
		y = fft0 ? Int.(-size(z,2)/2:size(z,2)/2-1) : (1:size(z,2)),
		xtick = (minimum(x) < 0 && maximum(x) > 0) ?
			[minfloor(x),0,maxceil(x)] : [minfloor(x),maxceil(x)],
		ytick = (minimum(y) < 0 && maximum(y) > 0) ?
			[minfloor(y),0,maxceil(y)] : [minfloor(y),maxceil(y)],
		yflip::Bool = nothing_else(jim_def[:yflip], minimum(y) >= 0),
		yreverse::Bool = nothing_else(jim_def[:yreverse], y[1] > y[end]),
		abswarn::Bool = jim_def[:abswarn], # ignored here
	)

	# because some backends require y to be in ascending order
	if yreverse
		y = reverse(y)
		z = reverse(z, dims=2)
	end

	n1,n2,n3 = size(z,1), size(z,2), size(z,3)
	if n3 == 1
		z = z[:,:,1] # revert to 2D
	end
	xy = (x, y)
	if ndims(z) > 2
		n1 += mosaic_npad
		n2 += mosaic_npad
		n3 = size(z,3)
		if ncol == 0
			ncol = Int(floor(sqrt(prod(size(z)[3:end]))))
		end
		z = mosaicview(z ; fillvalue=padval, ncol=ncol, npad=mosaic_npad)
		xy = () # no x,y for mosaic
	elseif fft0
		z = FFTView(z)[x,y]
	end

	if minimum(z) ≈ maximum(z) # uniform or nearly uniform image
		x = 1:size(z,1)
		y = 1:size(z,2)
		plot(aspect_ratio=aspect_ratio,
			xlim = [x[1], x[end]],
			ylim = [y[1], y[end]],
			title=title,
			yflip=yflip,
			xlabel=xlabel,
			ylabel=ylabel,
			xtick=xtick,
			ytick=ytick,
		)
		tmp = (minimum(z) == maximum(z)) ? "Uniform $(z[1])" :
			"Nearly uniform $((minimum(z),maximum(z)))"
		annotate!((sum(x)/length(x), sum(y)/length(y), tmp, :red))
	else
		heatmap(xy..., z', transpose=false,
			aspect_ratio=aspect_ratio,
			clim=clim,
			color=color,
			title=title,
			yflip=yflip,
			xlabel=xlabel,
			ylabel=ylabel,
			xtick=xtick,
			ytick=ytick,
		)
	end

	if n3 > 1 && line3plot # lines around each subimage
		m1 = (1+size(z,1)) / n1 # add one because of mosaicview non-edge
		m2 = (1+size(z,2)) / n2
		plot_box! = (ox, oy) -> plot!(
			ox .+ [0,1,1,0,0]*n1 .+ 0.0,
			oy .+ [0,0,1,1,0]*n2 .+ 0.0,
			line=jim_def[:line3type], label="")
		for ii=0:n3-1
			i1 = mod(ii, m1)
			i2 = floor(ii / m1)
			plot_box!(i1*n1, i2*n2)
		end
	end
	plot!()

end # jim

# handle case of complex image
function jim(z::AbstractArray{<:Number} ;
			abswarn::Bool = jim_def[:abswarn],
			kwargs...
		)

	if !(eltype(z) <: Real)
		abswarn && (@warn "magnitude at $(caller_name())")
		z = abs.(z)
	end
	jim(z ; kwargs...)
end

"""
    jim(z, title ; kwargs...)
"""
jim(z::AbstractArray{<:Number}, title::AbstractString ; kwargs...) =
	jim(z ; title=title, kwargs...)

"""
    jim(x, y, z ; kwargs...)
"""
jim(x, y, z ; kwargs...) = jim(z ; x=x, y=y, kwargs...)

"""
    jim(x, y, z, title ; kwargs...)
"""
jim(x, y, z, title::AbstractString ; kwargs...) =
	jim(z ; x=x, y=y, title=title, kwargs...)

"""
    jim()
return docstring if user calls `jim()` with no arguments
"""
function jim()
	@doc jim
end

#"""
#`jim(abswarn=false)`
#"""
#function jim( ; abswarn::Bool=jim_state_abswarn)
#	global jim_state_abswarn = abswarn
#end

"""
    jim(key, value)
set default value for one of the keys
"""
function jim(key::Symbol, value)
	global jim_def
	!haskey(jim_def, key) && throw(ArgumentError("no key $key"))
	jim_def[key] = value
end

"""
    jim(:test)

`jim(:keys)` return default keys

`jim(:defs)` return `Dict` of default keys / vals

`jim(:key)` return `Dict[key]` if possible
"""
function jim(test::Symbol)
	global jim_def
	if test === :keys
		return keys(jim_def)
	end
	if test === :defs
		return jim_def
	end
	if haskey(jim_def, test)
		return jim_def[test]
	end

	throw("symbol $test")
end