text:
#################### Post order traversal ####################
"""
    post_order(root::T, traversal::Vector{T})::Vector{T} where T<:AbstractNode

This function performs a post order traversal through the tree. It is assumed that `root` is the
root of the tree. Thus, if `root` is not the root, the subtree defined by the root `root` is
used for the post order traversal.

Returns vector of Nodes.

* `root` : root Node of tree.

* `traversal` : vector of Nodes; nodes are pushed to this vector as the tree is traversed.
"""
function post_order(root::T, traversal::Vector{T}) where T<:AbstractNode
   if root.nchild != 0
        for child in root.children
            post_order(child, traversal)
        end
   end # if
   push!(traversal, root)
   return traversal
end # function post_order_trav

"""
    post_order(root::T)::Vector{T} where T<:AbstractNode

This function does post order traversal. Only the root node needs to be supplied.

Returns vector of Nodes.

* `root` : root Node of tree.
"""
function post_order(root::T)::Vector{T} where T<:AbstractNode
    t::Vector{T} = []
    post_order(root, t)
    return t
end # function post_order

"""
    get_leaves(root::T, traversal::Vector{T})::Vector{T} where T<:AbstractNode

This function returns leaf nodes of a tree. It is assumed that `root` is the
root of the tree. Thus, if `root` is not the root, the subtree defined by the root `root` is
used.

Returns a vector of leaf Nodes.

* `root` : root Node of tree to traverse.

* `traversal` : vector of Nodes; leaf Nodes, once found, are pushed to this vector.

"""
function get_leaves(root::T, traversal::Vector{T})::Vector{T} where T<:AbstractNode
   if root.nchild != 0
        for child in root.children
            get_leaves(child, traversal)
        end
   else
       push!(traversal, root)
   end # if

   return traversal
end # function post_order_trav

"""
    get_leaves(root::T)::Vector{T} where T<:AbstractNode

This function returns the leaves of a tree. Only the root node needs to be supplied.

Returns vector of leaf Nodes.

* `root` : root Node of tree.
"""
function get_leaves(root::T)::Vector{T} where T<:AbstractNode
    t::Vector{T} = []
    get_leaves(root, t)
    return t
end # function post_order

#################### Pre order traversal ####################

"""
    pre_order(root::T, traversal::Vector{T})::Vector{T} where T<:AbstractNode

This function performs a pre order traversal through the tree. It is assumed that `root` is the
root of the tree. Thus, if `root` is not the root, the subtree defined by the root `root` is
used for the pre order traversal.

Returns vector of Nodes.

* `root` : root Node of tree.

* `traversal` : vector of Nodes; nodes are pushed to this vector as the tree is traversed.
"""
function pre_order(root::T, traversal::Vector{T})::Vector{T} where T<:AbstractNode
    push!(traversal, root)
    if root.nchild != 0
        for child in root.children
            pre_order(child, traversal)
        end
    end # if
    return traversal
end # function pre_order!

"""
    pre_order(root::T)::Vector{T} where T<:AbstractNode

This function does pre order traversal. Only the root node needs to be supplied.

Returns vector of Nodes.

* `root` : root Node of tree.
"""
function pre_order(root::T)::Vector{T} where T<:AbstractNode
    t::Vector{T} = []
    pre_order(root, t)
    return t
end # function pre_order

#################### Level order traversal ####################

"""
    level_order(node::T)::Array{T} where T<:AbstractNode

This function does level order traversal. Only the root node needs to be supplied.

Returns Array of Nodes.

* `node` : root Node of tree.
"""
function level_order(node::T)::Array{T} where T<:AbstractNode
    level = 1
    stack::Array{T} = []
    while level_traverse(node, level, stack)
        level += 1
    end # while
    stack
end # function level_order

"""
    level_traverse(node::T, level::Int64, stack::Array{T})::Bool where T <:AbstractNode

This function traverses a level of the tree specified through `node`. The level
is specified via the `level` argument and the nodes visited are stored in the
`stack`.
This function is intended as the internal worker for the level_order function.
"""
function level_traverse(node::T, level::Int64, stack::Array{T})::Bool where T <:AbstractNode

    if level == 1
        # level which needs to be traversed right now
        push!(stack, node)
        return true
    else
        # move down the tree to the correct level
        boolqueue = [false] # this is used to look for the correct level
        for child in node.children
            push!(boolqueue, level_traverse(child, level-1, stack))
        end # for
        return reduce(|, boolqueue)
    end # if
end # function level_traverse


text:
#
# @description Implementation of PKCS#7 padding functions.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export Pkcs7Padder
export pad, unpad

# @description Defines PKCS#7 the padding algorithm data struture.
struct Pkcs7Padder <: PaddingMethod end

# @description Computes the padded block for the given plaintext chunk.
#
# @param {Pkcs7Padder}   self       The padding algorithm.
# @param {Vector{UInt8}} chunk      The plaintext chunk to pad. This must be shorter than the block size.
# @param {Integer}       block_size The cipher block size.
#
# @return {Vector{UInt8}} The padded block.
function pad(self::Pkcs7Padder, chunk::Vector{UInt8}, block_size::Integer)
  len = block_size - (length(chunk) % block_size)
  return vcat(chunk, fill(convert(UInt8, len), len))
end

# @description Removes the pad bytes from the given plaintext block.
#
# @param {Pkcs7Padder}   self  The padding algorithm.
# @param {Vector{UInt8}} block The padded block.
#
# @return {Vector{UInt8}} The block without pad.
function unpad(self::Pkcs7Padder, block::Vector{UInt8})
  if block[end] == 0x00 || block[end] > length(block)
    throw(PaddingError())
  end

  if findfirst(x -> x != block[end], block[(end - block[end] + 1) : (end - 1)]) != 0
    throw(PaddingError())
  end

  return block[1 : (end - block[end])]
end

ext: jl


text:
function show_args(args)
  @show args
end

# Example of running this routine:

include("run_hmc_background_student_rescale.jl")

# File containing the transit times:
#fname = "../../data/T1_timings_20191203.txt"
fname = "timing_simon.txt"

felements = "../../data/elements_noprior_students.txt"

# Name an output file:
#foutput = string("T1_run_hmc_student_ecc_lgndof_V1exp2nuinv_nstep2000_eps0.1_nleap20_",show_args(ARGS)[1],".jld2")
foutput = string("T1_run_hmc_student_ecc_lgndof_V1exp2nuinv_nstep10_eps0.1_nleap20_",show_args(ARGS)[1],".jld2")

# Call the likelihood profile function:
#run_hmc(fname,felements,foutput,2000,0.1,20,false)
run_hmc(fname,felements,foutput,10,0.1,20,false)


text:
module A
	export foo, B

	function foo()
		"hello"
	end

	function bar()
		"goodbye"
	end

	module B
		baz = 5
	end
end

module C
	export quux

	function quux()
		"kaboom"
	end
end


text:
#=
f is per t in planning graph
x0 goes with y1, xt-1 goes with yt
no-op are part of y

1. init and goal constraints
=#

function symb_get_fluent_terminal_indices(graph, max_nf; loc=:init)
    indices = []
    not_ind = []
    T = graph.num_levels - 1
    if loc == :init
        for i=1:length(graph.initprops)
            ind = findall(x->x == graph.initprops[i], graph.props[0])[1]
            push!(indices, ind)
        end
    else
        for i=1:length(graph.goalprops)
            ind = findall(x->x == graph.goalprops[i], graph.props[T])[1]
            push!(indices, ind)
        end
    end
    not_ind = [i for i=1:max_nf if !(i in indices)]
    return indices, not_ind
end

function symb_get_actions_precondition_sat(fluent, level, graph)
    actions = graph.acts[level]
    # actions = [act for act in graph.acts[level] if act.name != :NoOp]
    act_inds = []
    for i=1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_prec
            push!(act_inds, i)
        end
    end
    return act_inds
end

function symb_get_actions_effect_sat(fluent, level, graph)
    actions = graph.acts[level]
    # actions = [act for act in graph.acts[level] if act.name != :NoOp]
    act_inds = []
    for i=1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_eff
            push!(act_inds, i)
        end
    end
    return act_inds
end

function symb_get_fluent_index(fluent, level, graph)
    for i=1:length(graph.props[level])
        if graph.props[level][i] == fluent
            return i
        end
    end
    return nothing
end

function symb_get_mutex_actions(fluent, level, graph)
    pref = symb_get_actions_precondition_sat(fluent, level, graph)
    preneff = union(pref, symb_get_actions_effect_sat(fluent, level, graph))
    act_mutexes = []
    [push!(act_mutexes, [i,j]) for i in pref for j in preneff if i!=j && !([j,i] in act_mutexes)]

    return act_mutexes
end

function symb_get_mutex_indices(level, graph)
    act_mutexes = []
    for pair in graph.μacts[level]
        i = findall(x->x == pair[1], graph.acts[level])[1]
        j = findall(x->x == pair[2], graph.acts[level])[1]

        if !([j,i] in act_mutexes) push!(act_mutexes, [i,j]) end
    end
    return act_mutexes
end

function symb_test_mutexes(domain_name, problem_name; max_levels=10)
    domain = load_domain(domain_name)
    problem = load_problem(problem_name)
    graph = create_symbolic_graph(domain, problem; max_levels=max_levels)
    T = graph.num_levels
    for t = 1:T-1 #xt, yt, bk: xt+1
        println("\n\nLevel: ",t)
        println("------------------")
        for f=1:length(graph.props[t-1])
            fluent = graph.props[t-1][f]
            # println("fluent: ",fluent)
            m1 = symb_get_mutex_actions(fluent, t, graph)
            m2 = symb_get_mutex_actions2(fluent, t, graph)
            # println("m1: ",m1)
            println("m2: ",m2)
            println(" ")
        end
        println("gp mutexes:")
        println(symb_get_mutex_indices(t, graph))
    end

end

function pure_symbolic_solve(domain_name, problem_name; max_levels=10)
    domain = load_domain(domain_name)
    problem = load_problem(problem_name)
    graph = create_graph(domain, problem; max_levels=max_levels)

    model = Model(Gurobi.Optimizer)
    set_silent(model)

    T = graph.num_levels
    max_nf = max([length(graph.props[i-1]) for i=1:T]...)
    max_na = max([length(graph.acts[i]) for i=1:T-1]...)
    # max_na = max([length([act for act in graph.acts[i] if act.name != :NoOp]) for i=1:T-1]...)

    @variable(model, x[i=1:max_nf, j=1:T], Bin)
    @variable(model, y[i=1:max_na, j=1:T-1], Bin)

    inits, notinits = symb_get_fluent_terminal_indices(graph, max_nf; loc=:init)
    [@constraint(model, x[f, 1] == 1) for f in inits]
    [@constraint(model, x[f, 1] == 0) for f in notinits]

    goals, _ = symb_get_fluent_terminal_indices(graph, max_nf; loc=:goal)
    [@constraint(model, x[f, T] == 1) for f in goals]
    [@constraint(model, x[f, T] == 0) for f = 1:max_nf if !(f in goals)]

    for t = 1:T-1 #xt, yt, bk: xt+1
        for f=1:length(graph.props[t-1])
            fluent = graph.props[t-1][f]
            act_inds = symb_get_actions_precondition_sat(fluent, t, graph)
            [@constraint(model, y[a,t] <= x[f, t]) for a in act_inds]
        end
    end

    for t = 1:T-1 #xt, yt, bk: xt+1
        for f=1:length(graph.props[t])
            fluent = graph.props[t][f]
            act_inds = symb_get_actions_effect_sat(fluent, t, graph)
            fnext = symb_get_fluent_index(fluent, t, graph)
            if !isnothing(fnext) && !isempty(act_inds)
                @constraint(model, x[fnext, t+1] <= sum([y[a, t] for a in act_inds]))
            end
        end
        # na = length([act for act in graph.acts[t] if act.name != :NoOp])
        na = length(graph.acts[t])
        @constraint(model, y[na+1:max_na, t] .== 0)
    end

    for t = 1:T-1
        for f=1:length(graph.props[t])
            fluent = graph.props[t][f]
            act_pairs = symb_get_mutex_actions(fluent, t, graph)
            [@constraint(model, y[a[1], t] + y[a[2], t] <= 1) for a in act_pairs]
        end
    end

    @objective(model, Min, sum(y))
    @time optimize!(model)
    sol = value.(y)
    render_action(sol, graph)
    # graph

end

function render_action(y, graph)
    T = graph.num_levels-1
    actions = []
    for t=1:T
        # try
        ind = findall(x->x == 1.0, y[1:end,t])[1]
        act  = graph.acts[t][ind]
        # act = all_actions[ind]
        push!(actions, (act.name, act.args))
        # catch
        # end
    end
    # println(actions)
    return actions

end


text:
#=
This file reads Dynare .mod files, and extracts the following information:
  variables
  exo_variables
  parameters
  calibration
  model
  initval
  endval
  shocks
=#

using DataStructures

function modfile_parser(mod_file_name)
  f = open(mod_file_name)
  lines = readlines(f)  # reads in line by line
  close(f)

  # Remove lines beginning with comments (%, //, /*, \*)
  filter!(x -> ~ismatch(r"^\s*%", x), lines)
  filter!(x -> ~ismatch(r"^\s*\/\/", x), lines)
  filter!(x -> ~ismatch(r"^\s*\/\*", x), lines)

  # Remove end-of-line comments
  for ln = 1:length(lines)
    lines[ln] = split(lines[ln], r"%.*$")[1]    # %  comments
    lines[ln] = split(lines[ln], r"\/\/.*$")[1]  # \\ comments
    lines[ln] = split(lines[ln], r"/\*.*$")[1]   # \* comments
  end

  # Smoosh text back together
  text = join(lines)

  # Remove new lines, spaces, carriage returns,
  text = replace(text, "\t", " ")
  text = replace(text, "\r", " ")
  text = replace(text, "\n", " ")

  # Get variable names
  tmp = match(r"var\s(.*?);", text)
  endogenous = split(tmp[1], " ")
  filter!(x -> ismatch(r"\S", x), endogenous) # Remove lines if no chars present

  # Get exogenous variable names
  tmp = match(r"varexo\s(.*?);", text)
  shocks = split(tmp[1], " ")
  filter!(x -> ismatch(r"\S", x), shocks) # Remove lines if no chars present

  # Get parameter names
  tmp = match(r"parameters\s(.*?);", text)
  parameters = split(tmp[1], " ")
  filter!(x -> ismatch(r"\S", x), parameters) # Remove lines if no chars present

  # get equations, fill a list
  tmp  = match(r"model;\s(.*?)end;", text)
  equations = split(tmp[1], ";")
  filter!(x -> ismatch(r"\S", x), equations) # Remove lines if no chars present
  for ln = 1:length(equations)
    equations[ln] = strip(equations[ln])
  end

  # Get calibration values, fill a dictionary
  parameter_values = OrderedDict()
  tmp  = match(r"parameters\s(.*?);(.*)model", text)
  tmp = split(tmp[2], ";")
  filter!(x -> ismatch(r"\S", x), tmp) # Remove lines if no chars present
  for ln = 1:length(tmp)
    key = strip(match(r"(.*)=", tmp[ln])[1])
    entry = strip(match(r"=(.*)", tmp[ln])[1])
    parameter_values[key] =  entry
  end

  # Get initial values, fill a dictionary
  initval = OrderedDict()
  tmp  = match(r"initval;(.*?)end;", text)
  if tmp != nothing
    tmp = split(tmp[1], ";")
    filter!(x -> ismatch(r"\S", x), tmp) # Remove lines if no chars present
    for ln = 1:length(tmp)
      key = strip(match(r"(.*)=", tmp[ln])[1])
      entry = strip(match(r"=(.*)", tmp[ln])[1])
      initval[key] =  entry
    end
  end

  # Get end values, fill a dictionary
  endval = OrderedDict()
  tmp  = match(r"endval;(.*?)end;", text)
  if tmp != nothing
    tmp = split(tmp[1], ";")
    filter!(x -> ismatch(r"\S", x), tmp) # Remove lines if no chars present
    for ln = 1:length(tmp)
      key = strip(match(r"(.*)=", tmp[ln])[1])
      entry = strip(match(r"=(.*)", tmp[ln])[1])
      endval[key] =  entry
    end
  end

  # Get the calibrated shock values, fill matrix
  shockvaldict = Dict()
  shock_matrix = fill("0", length(shocks), length(shocks))
  tmp = match(r"shocks;(.*?)(.*?)end;", text)
  if tmp != nothing
    if contains(tmp[2], "stderr")
      tmpkey = matchall(r"var(.*?);", tmp[2])
      tmpentry = matchall(r"stderr(.*?);", tmp[2])
      for ln = 1:length(tmpkey)
       key = strip(match(r"var\s(.*)", tmpkey[ln])[1])
       key = match(r"(.*);", key)[1]
       entry = strip(match(r"stderr\s(.*)", tmpentry[ln])[1])
       entry = match(r"(.*);", entry)[1]
       shockvaldict[key] =  entry
      end
    else
      tmp = split(tmp[2], ";")
      filter!(x -> ismatch(r"\S", x), tmp) # Remove lines if no chars present
      for ln = 1:length(tmp)
        key = match(r"var\s(.*)\s=", tmp[ln])[1]
        entry = match(r"=\s(.*)", tmp[ln])[1]
        shockvaldict[key] =  entry
      end
    end
    # Fill matrix in order of shocks in "shocks" dictionary
    for ln = 1:length(shocks)
      if haskey(shockvaldict, shocks[ln])
        shock_matrix[ln, ln] = shockvaldict[shocks[ln]]
      else
        shock_matrix[ln, ln] = "0.0"
      end
    end
  end

  data = Dict()
  data["symbols"] = Dict("endogenous"=>endogenous, "exogenous"=>shocks, "parameters"=>parameters)
  data["equations"] = equations
  data["calibration"] = merge(parameter_values, initval)
  exogenous = Dict("tag"=>"Normal", "sigma"=>shock_matrix)
  data["exogenous"] = exogenous
  data["optional"] = Dict("endval"=>endval)

  return data

end


text:
include("./params.jl")

function fitgamma(isi)

    epsilon = 1e-6;

    n = length(isi);
    T = sum(isi);

    lambda = n / T;

    alpha = 0.2
    beta = 1.0

    eta = 0.001

    dll_dalpha = epsilon * 10
    dll_dbeta = epsilon * 10

    while max(abs(dll_dalpha), abs(dll_dbeta)) > epsilon

        dll_dalpha = - n / alpha + n * log(beta) + sum(log(isi)) - n * digamma(alpha);
        dll_dbeta = n * (1+alpha) / beta - T;

        alpha = alpha + eta*dll_dalpha;
        beta = beta + eta*dll_dbeta;

    end

    k = alpha + 1
    lambda = beta / k

    lambda, k

end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
using StatsBase

type Node
    id::Int16
    item_ids::Array{Int16,1}
    transactions::BitArray{1}
    children::Array{Node,1}
    mother::Node
    supp::Int

    function Node(id::Int16, item_ids::Array{Int16,1}, transactions::BitArray{1})
        nd = new()
        nd.id = id
        nd.item_ids = item_ids
        nd.transactions = transactions
        nd.children = Array{Node,1}(0)
        return nd
    end

    function Node(id::Int16, item_ids::Array{Int16,1}, transactions::BitArray{1}, mother::Node, supp::Int)
        nd = new()
        nd.id = id
        nd.item_ids = item_ids
        nd.transactions = transactions
        nd.children = Array{Node,1}(0)
        nd.mother = mother
        nd.supp = supp
        return nd
    end
end

type Rule
    p::Array{Int16,1}
    q::Int16
    supp::Float64
    conf::Float64
    lift::Float64

    function Rule(node::Node, mask::BitArray{1}, supp_dict::Dict{Array{Int16,1}, Int}, num_transacts::Int)

        rule = new()
        rule.p = node.item_ids[mask]
        rule.supp = node.supp/num_transacts
        rule.conf = rule.supp/supp_dict[node.item_ids[mask]]
        unmask = !mask
        q_idx = findfirst(unmask)
        rule.q = node.item_ids[q_idx]
        rule.lift = rule.conf/supp_dict[node.item_ids[unmask]]
        return rule
    end
end

# @code_warntype Node(Int16(1), Int16[1], trues(3))
n1 = Node(Int16(1), Int16[1], trues(3))

# @code_warntype Node(Int16(1), Int16[1, 2], trues(3), n1)
n2 = Node(Int16(1), Int16[1, 2], trues(3), n1, 1)
n3 = Node(Int16(1), Int16[1, 3], trues(3), n1, 1)
n4 = Node(Int16(1), Int16[1, 4], trues(3), n1, 1)
n5 = Node(Int16(1), Int16[1, 5], trues(3), n1, 1)
n6 = Node(Int16(1), Int16[1, 6], trues(3), n1, 1)
n7 = Node(Int16(1), Int16[2, 3], trues(3), n1, 1)
n8 = Node(Int16(1), Int16[2, 4], trues(3), n1, 1)
n9 = Node(Int16(1), Int16[2, 5], trues(3), n1, 1)
n10 = Node(Int16(1), Int16[2, 6], trues(3), n1, 1)

push!(n1.children, n2)
push!(n1.children, n3)
push!(n1.children, n4)
push!(n1.children, n5)
push!(n1.children, n6)
push!(n1.children, n7)
push!(n1.children, n8)
push!(n1.children, n9)
push!(n1.children, n10)

function has_children(nd::Node)
    res = length(nd.children) > 0
    res
end

# @code_warntype has_children(n1)

function younger_siblings(nd::Node)
    n_sibs = length(nd.mother.children)
    return view(nd.mother.children, (nd.id + 1):n_sibs)
end

# @code_warntype younger_siblings(n1.children[1])
younger_siblings(n1.children[1])

function update_support_cnt!(supp_dict::Dict, nd::Node)
    supp_dict[nd.item_ids] = nd.supp
end

# This function is used internally and is the workhorse of the frequent()
# function, which generates a frequent itemset tree. The growtree!() function
# builds up the frequent itemset tree recursively.
function growtree!(nd::Node, minsupp, k, maxdepth)
    sibs = younger_siblings(nd)

    for j = 1:length(sibs)
        transacts = nd.transactions & sibs[j].transactions
        supp = sum(transacts)

        if supp ≥ minsupp
            items = zeros(Int16, k)
            items[1:k-1] = nd.item_ids[1:k-1]
            items[end] = sibs[j].item_ids[end]

            child = Node(Int16(j), items, transacts, nd, supp)
            push!(nd.children, child)
        end
    end
    # Recurse on newly created children
    maxdepth -= 1
    if maxdepth > 1
        for child in nd.children
            growtree!(child, minsupp, k+1, maxdepth)
        end
    end
end

# @code_warntype growtree!(n2, 1, 3, 3)
growtree!(n2, 1, 3, 3)

function get_unique_items{M}(all_transacts::Array{Array{M, 1}, 1})
    dict = Dict{M, Bool}()

    for t in all_transacts
        for i in t
            dict[i] = true
        end
    end
    return collect(keys(dict))
end

t = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e"]]

@code_warntype get_unique_items(t);
@time get_unique_items(t);

# This function is used internally by the frequent() function to create the
# initial bitarrays used to represent the first "children" in the itemset tree.
function occurrence(all_transacts::Array{Array{String, 1}, 1}, uniq_items::Array{String, 1})
    n = length(all_transacts)
    p = length(uniq_items)

    itm_pos = Dict(zip(uniq_items, 1:p))
    res = falses(n, p)
    for i = 1:n
        for itm in all_transacts[i]
            j = itm_pos[itm]
            res[i, j] = true
        end
    end
    res
end

unq = get_unique_items(t)
@code_warntype occurrence(t, unq)
@time occurrence(t, unq)

"""
    frequent(all_transacts, minsupp, maxdepth)

This function creates a frequent itemset tree from an array of transactions.
The tree is built recursively using calls to the growtree!() function. The
`minsupp` and `maxdepth` parameters control the minimum support needed for an
itemset to be called "frequent", and the max depth of the tree, respectively
"""
function frequent(all_transacts::Array{Array{String, 1}, 1}, minsupp, maxdepth)
    uniq_items = get_unique_items(all_transacts)
    sort!(uniq_items)

    occ = occurrence(all_transacts, uniq_items)
    root = Node(Int16(1), Int16[-1], BitArray(0))
    n_items = length(uniq_items)

    # This loop creates 1-item nodes (i.e., first children)
    for j = 1:n_items
        supp = sum(occ[:, j])
        if supp ≥ minsupp
            nd = Node(Int16(j), Int16[j], occ[:, j], root, supp)
            push!(root.children, nd)
        end
    end
    n_kids = length(root.children)

    # Grow nodes in breadth-first manner
    for j = 1:n_kids
        growtree!(root.children[j], minsupp, 2, maxdepth)
    end
    root
end

t = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e", "f"]]

@code_warntype frequent(t, 0.01, 3)
xtree1 = frequent(t, 0.01, 4)

function prettyprint(node::Node, k::Int = 0)
    if has_children(node)
        for nd in node.children
            print("k = $(k + 1): ")
            println(nd.item_ids)
        end
        for nd in node.children
            prettyprint(nd, k+1)
        end
    end
end

prettyprint(xtree1)

function randstring(n::Int, len::Int = 16)
    vals = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m",
            "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"]
    upper = map(uppercase, vals)
    append!(vals, upper)
    append!(vals, map(string, 0:9))
    res = Array{String,1}(n)
    for i = 1:n
        res[i] = join(rand(vals, len))
    end
    res
end

itemlist = randstring(100);

n = 10_000
m = 25              # number of items in transactions
t = [sample(itemlist, m, replace = false) for _ in 1:n];

# @code_warntype frequent(t, 1)
@time unq2 = get_unique_items(t);
@time occ2 = occurrence(t, unq2);
@time f = frequent(t, round(Int, 0.01*n), m);

function grow_support_dict!(supp_cnt::Dict{Array{Int16,1}, Int}, node::Node)
    if has_children(node)
        for nd in node.children
            update_support_cnt!(supp_cnt, nd)
            grow_support_dict!(supp_cnt, nd)
        end
    end
end

# This function generates a dictionary whose keys are the frequent
# itemsets (their integer represenations, actually), and whose values
# are the support count for the given itemset. This function is used
# for computing support, confidence, and lift of association rules.
function gen_support_dict(root::Node)
    supp_cnt = Dict{Array{Int16, 1}, Int}()
    grow_support_dict!(supp_cnt, root)
    return supp_cnt
end

t1 = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e", "f"]]

@code_warntype frequent(t1, 1, 3)
xtree1 = frequent(t1, 1, 4);
@code_warntype gen_support_dict(xtree1)
xsup = gen_support_dict(xtree1)

function gen_rules(node::Node, supp_dict::Dict{Array{Int16,1}, Int}, k, num_transacts)
    mask = trues(k)
    rules = Array{Rule, 1}(k)
    for i = 1:k
        mask[i] = false
        if i > 1
            mask[i-1] = true
        end
        rules[i] = Rule(node, mask, supp_dict, num_transacts)
    end
    return rules
end

@code_warntype gen_rules(xtree1.children[1].children[1].children[1], xsup, 3, 8)

# function compute_metrics(root::Node)
#     # supp_dict = gen_support_dict(root)
# end

# function get_cousins(node::Node)
#     cousins = Array{Node,1}(0)
#     if isdefined(node, :mother) && isdefined(node.mother, :mother)
#         for aunt in younger_siblings(node.mother)
#             for nd in aunt.children
#                 push!(cousins, nd)
#             end
#         end
#     end
#     cousins
# end

# get_cousins(xtree1.children[1])

# function prettyprint2(node::Node, k::Int)
#     if has_children(node)
#         for nd in node.children
#             print("k = $(k + 1): ")
#             println(nd.item_ids)
#         end
#         for nd in node.children
#             if (k + 1) ≥ 1 && has_children(nd)
#                 for nd in vcat(nd.children, get_cousins(nd.children[1]))
#                     print("k = $(k + 2): ")
#                     println(nd.item_ids)
#                 end

#             end
#         end
#         prettyprint2(node.children[1], k + 1)
#     end
# end

# prettyprint2(xtree1, 0)


text:
@testset "Download" begin

    pfam_code = "PF11591"
    @test_throws ErrorException downloadpfam("2vqc")
    filename = downloadpfam(pfam_code, filename=tempname()*".gz")
    try
        aln = read(filename, Stockholm)
        if size(aln) == (6,34)
            @test getannotfile(aln, "ID") == "2Fe-2S_Ferredox"
        end
    finally
        rm(filename)
    end
end

@testset "PDB code" begin

    msa = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm)
    @test getseq2pdb(msa)["F112_SSV1/3-112"] == [("2VQC","A")]
end

@testset "Mapping PDB/Pfam" begin

    msa_file   = joinpath(pwd(), "data", "PF09645_full.stockholm")
    sifts_file = joinpath(pwd(), "data", "2vqc.xml.gz")
    pdb_file   = joinpath(pwd(), "data", "2VQC.xml")
    msa = read(msa_file, Stockholm, generatemapping=true, useidcoordinates=true)
    cmap = msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645", sifts_file)
    res = residuesdict(read(pdb_file, PDBML), "1", "A", "ATOM", All)

    #     -45              20 pdb
    #.....QTLNSYKMAEIMYKILEK  msa seq
    #     123456789012345678  msa col
    #     345678901234567890  uniprot 3-20
    #    ****              *
    #12345678901234567890123  ColMap

    @test_throws KeyError cmap[5]  # insert
    @test cmap[6]  == ""           # missing
    @test cmap[7]  == "4"
    @test cmap[8]  == "5"
    @test cmap[23] == "20"

    #.....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ..... msa seq
    #.....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..... pdb ss
    #                                                                                111111111111111111111 ColMap hundreds
    #         111111111122222222223333333333444444444455555555556666666666777777777788888888889999999999000000000011111111112 ColMap tens
    #123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 ColMap ones
    #                                                                           **                                     **

    @test_throws KeyError cmap[116]   # insert
    @test cmap[115] == ""             # missing
    @test cmap[77]  == ""             # missing
    @test cmap[76]  == "73"

    @testset "Residues" begin

        msares = msaresidues(msa, res, cmap)

        @test_throws KeyError msares[5]   # insert
        @test_throws KeyError msares[6]   # missing
        @test msares[7].id.name  == "THR" # T
        @test msares[8].id.name  == "LEU" # L
        @test msares[23].id.name == "LYS" # K

        @test_throws KeyError msares[116]   # insert
        @test_throws KeyError msares[117]   # missing
        @test_throws KeyError msares[77]    # missing
        @test msares[76].id.name == "LYS"   # K
    end

    @testset "Contacts" begin

        @test msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645",
                                   sifts_file, strict=true, checkpdbname=true) == cmap

        @test length(res) == 70
        @test length(unique(values(cmap))) == 71

        #     -45              20 pdb
        #.....QTLNSYKMAEIMYKILEK  msa seq
        #     123456789012345678  msa col
        #12345678901234567890123  ColMap
        #     345678901234567890  uniprot 3-20
        #     ***              *

        @test_throws KeyError res["3"]
        @test res[cmap[7]].id.number == "4"

        contacts = msacontacts(msa, res, cmap, 6.05)
        missings = sum(Int.(isnan.(contacts)), dims=1)

        @test size(contacts) == (110,110)

        @test missings[1]  == 110              # missing
        @test missings[2]  == (110 - (70 - 1)) # 70 residues, 1 diagonal NaN
        @test missings[3]  == (110 - (70 - 1))
        @test missings[18] == (110 - (70 - 1))

        @test missings[110] == 110                 # missing
        @test missings[72]  == 110                 # missing
        @test missings[71]  == (110 - (70 - 1))

        ncontacts = sum(Int.(contacts .== 1.0), dims=1)

        @test ncontacts[1]  == 0
        @test ncontacts[2]  == 2
        @test ncontacts[3]  == 6

        @testset "using msaresidues" begin
        # Test MSA contact map using PDBResidues from the MSA

            msares = msaresidues(msa, res, cmap)
            ncol = ncolumns(msa)
            colmap = getcolumnmapping(msa)

            for i in 1:(ncol-1), j in (i+1):ncol
                if (colmap[i] in keys(msares)) && (colmap[j] in keys(msares))
                    @test (contacts[i,j] .== 1.0) ==
                           contact(msares[colmap[i]], msares[colmap[j]], 6.05)
                end
            end

            @testset "hasresidues" begin

                mask = hasresidues(msa, cmap)

                @test mask[1] == false
                @test mask[2] == true
                @test sum(mask) == length(msares)
            end
        end

        @testset "AUC" begin

            @test round(AUC(buslje09(msa,lambda=0.05,threshold=62.0,samples=0)[2],contacts), digits=4) == 0.5291
        end
    end
end

@testset "AUC" begin

    ntru = 90
    nfal = 100
    score_tru = Float16[ 2 + 2x for x in randn(ntru)]
    score_fal = Float16[-2 + 2x for x in randn(nfal)]
    msacontacts = NamedArray(PairwiseListMatrix{Float16,false,Vector{Float16}}(
        vcat(ones(Float16, ntru), zeros(Float16, nfal)), Float16[1.0 for x in 1:20], 20))
    score = NamedArray(PairwiseListMatrix{Float16,false,Vector{Float16}}(
        vcat(score_tru, score_fal), Float16[NaN for x in 1:20], 20))
    correct = 1.0 - auc(roc(score_tru, score_fal))

    @test AUC(score, msacontacts) == correct

    score[1:end,2] .= NaN
    msacontacts[1:end,3] .= NaN
    list_score = getlist(getarray(score))
    list_contact = getlist(getarray(msacontacts))
    n_values = length(list_score)
    @test n_values == length(list_contact)
    tar = [ list_score[i] for i in 1:n_values if
            !isnan(list_score[i]) & !isnan(list_score[i]) & (list_contact[i] == 1.0) ]
    non = [ list_score[i] for i in 1:n_values if
            !isnan(list_score[i]) & !isnan(list_score[i]) & (list_contact[i] == 0.0) ]
    @test AUC(score, msacontacts) == AUC(roc(tar, non))
end


text:
using SCS
const SCS_PRIMAL_OPT = SCS.Optimizer(verbose = false)
const SCS_DUAL_OPT = DualOptimizer(SCS.Optimizer(verbose = false))
const SCS_PRIMAL_FACTORY = with_optimizer(SCS.Optimizer, verbose = false)
const SCS_DUAL_FACTORY = with_optimizer(DualOptimizer, SCS.Optimizer(verbose = false))

push!(primal_power_cone_factory, SCS_PRIMAL_FACTORY)
push!(dual_power_cone_factory, SCS_DUAL_FACTORY)
push!(dual_power_cone_optimizer, SCS_DUAL_OPT)
push!(primal_power_cone_optimizer, SCS_PRIMAL_OPT)

@testset "SCS Exponential Cone Problems" begin
    list_of_exp_problems = [
        exp1_test,
        exp2_test
    ]
    test_strong_duality(list_of_exp_problems, SCS_PRIMAL_FACTORY)
end

@testset "SCS Power Cone Problems" begin
    list_of_pow_problems = [
        pow1_test,
        pow2_test
    ]
    test_strong_duality(list_of_pow_problems, SCS_PRIMAL_FACTORY; atol = 1e-3, rtol = 1e-3)
end


text:
type PivotMetadata
  rowNames::Vector{String}
  colNames::Vector{String}
  rowLevels::Vector{Vector{String}}
  colLevels::Vector{Vector{String}}
  statisticNames::Vector
  statisticInCols::Bool

  function PivotMetadata(rowNames::Vector{String}, colNames::Vector{String},
                         rowLevels, colLevels, statisticNames, statisticInCols)
    nRowLevel, nColLevel = length(rowLevels), length(colLevels)
    nRows,     nCols     = length(rowNames),  length(colNames)

    nRowLevel == nRows ? nothing :
      error("dim mismatch. nrowLevels != nrows, $mRowLevel != $nRows")
    nColLevel == nCols ? nothing :
      error("dim mismatch. ncolLevels != ncols, $nColLevel != $mCols")

    return new(rowNames, colNames, rowLevels, colLevels, statisticNames, statisticInCols)
  end
end

function _sanitizelabels(labels::Vector)
  if length(labels) == 0
    return String[""]
  else
    return [string(lbl) for lbl in labels]
  end
end

function PivotMetadata(rowNames, colNames,
                       rowLevels, colLevels,
                       statisticNames, statisticInCols)
  rowNames = _sanitizelabels(rowNames)
  colNames = _sanitizelabels(colNames)
  return PivotMetadata(rowNames, colNames,
                       rowLevels, colLevels,
                       statisticNames, statisticInCols)
end

function ncols(pm::PivotMetadata)
  levels = copy(pm.colLevels)
  if pm.statisticInCols
    push!(levels, statisticNames)
  end
  return prod([length(level) for level in levels])
end

function nrows(pm::PivotMetadata)
  levels = copy(pm.rowLevels)
  if !pm.statisticInCols
    push!(levels, statisticNames)
  end
  return prod([length(level) for level in levels])
end

max_stars_repo_name: cpritcha/Pivots

max_issues_repo_name: cpritcha/Pivots

max_forks_repo_name: cpritcha/Pivots


text:
"""
    Gibbs(n_iters, algs...)

Compositional MCMC interface. Gibbs sampling combines one or more
sampling algorithms, each of which samples from a different set of
variables in a model.

Example:
```julia
@model gibbs_example(x) = begin
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
        ...
end

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(1000, HMC(1, 0.2, 3, :v1), PG(20, 1, :v2))
```

Tips:
- `HMC` and `NUTS` are fast samplers, and can throw off particle-based
methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including
more particles in the particle sampler.
"""
mutable struct Gibbs{A} <: InferenceAlgorithm
    n_iters   ::  Int     # number of Gibbs iterations
    algs      ::  A   # component sampling algorithms
    thin      ::  Bool    # if thinning to output only after a whole Gibbs sweep
    gid       ::  Int
end
Gibbs(n_iters::Int, algs...; thin=true) = Gibbs(n_iters, algs, thin, 0)
Gibbs(alg::Gibbs, new_gid) = Gibbs(alg.n_iters, alg.algs, alg.thin, new_gid)

const GibbsComponent = Union{Hamiltonian,MH,PG}

function Sampler(alg::Gibbs, model::Model)
    n_samplers = length(alg.algs)
    samplers = Array{Sampler}(undef, n_samplers)

    space = Set{Symbol}()

    for i in 1:n_samplers
        sub_alg = alg.algs[i]
        if isa(sub_alg, GibbsComponent)
            samplers[i] = Sampler(typeof(sub_alg)(sub_alg, i), model)
        else
            @error("[Gibbs] unsupport base sampling algorithm $alg")
        end
        space = union(space, sub_alg.space)
    end

    # Sanity check for space
    @assert issubset(Set(pvars(model)), space) "[Gibbs] symbols specified to samplers ($space) doesn't cover the model parameters ($(Set(pvars(model))))"

    if Set(pvars(model)) != space
        @warn("[Gibbs] extra parameters specified by samplers don't exist in model: $(setdiff(space, Set(pvars(model))))")
    end

    info = Dict{Symbol, Any}()
    info[:samplers] = samplers

    Sampler(alg, info)
end

function sample(
                model::Model,
                alg::Gibbs;
                save_state=false,         # flag for state saving
                resume_from=nothing,      # chain to continue
                reuse_spl_n=0             # flag for spl re-using
                )

    # Init the (master) Gibbs sampler
    spl = reuse_spl_n > 0 ? resume_from.info[:spl] : Sampler(alg, model)

    @assert typeof(spl.alg) == typeof(alg) "[Turing] alg type mismatch; please use resume() to re-use spl"

    # Initialize samples
    sub_sample_n = []
    for sub_alg in alg.algs
        if isa(sub_alg, GibbsComponent)
            push!(sub_sample_n, sub_alg.n_iters)
        else
            @error("[Gibbs] unsupport base sampling algorithm $alg")
        end
    end

    # Compute the number of samples to store
    n = reuse_spl_n > 0 ? reuse_spl_n : alg.n_iters
    sample_n = n * (alg.thin ? 1 : sum(sub_sample_n))

    # Init samples
    time_total = zero(Float64)
    samples = Array{Sample}(undef, sample_n)
    weight = 1 / sample_n
    for i = 1:sample_n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    # Init parameters
    varInfo = if resume_from == nothing
        vi_ = VarInfo()
        model(vi_, HamiltonianRobustInit())
        vi_
    else
        resume_from.info[:vi]
    end

    n = spl.alg.n_iters; i_thin = 1

    # Gibbs steps
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[Gibbs] Sampling...", 0))
    for i = 1:n
        @debug "Gibbs stepping..."

        time_elapsed = zero(Float64)
        lp = nothing; epsilon = nothing; lf_num = nothing; eval_num = nothing

        for local_spl in spl.info[:samplers]
            last_spl = local_spl
      # PROGRESS[] && haskey(spl.info, :progress) && (local_spl.info[:progress] = spl.info[:progress])

            @debug "$(typeof(local_spl)) stepping..."

            if isa(local_spl.alg, GibbsComponent)
                for _ = 1:local_spl.alg.n_iters
                    @debug "recording old θ..."
                    time_elapsed_thin = @elapsed varInfo, is_accept = step(model, local_spl, varInfo, Val(i==1))

                    if ~spl.alg.thin
                        samples[i_thin].value = Sample(varInfo).value
                        samples[i_thin].value[:elapsed] = time_elapsed_thin
                        if ~isa(local_spl.alg, Hamiltonian)
                            # If statement below is true if there is a HMC component which provides lp and epsilon
                            if lp != nothing samples[i_thin].value[:lp] = lp end
                            if epsilon != nothing samples[i_thin].value[:epsilon] = epsilon end
                            if lf_num != nothing samples[i_thin].value[:lf_num] = lf_num end
                            if eval_num != nothing samples[i_thin].value[:eval_num] = eval_num end
                        end
                        i_thin += 1
                    end
                    time_elapsed += time_elapsed_thin
                end

                if isa(local_spl.alg, Hamiltonian)
                    lp = getlogp(varInfo)
                    epsilon = getss(local_spl.info[:wum])
                    lf_num = local_spl.info[:lf_num]
                    eval_num = local_spl.info[:eval_num]
                end
            else
                @error("[Gibbs] unsupport base sampler $local_spl")
            end
        end

        time_total += time_elapsed

        if spl.alg.thin
            samples[i].value = Sample(varInfo).value
            samples[i].value[:elapsed] = time_elapsed
            # If statement below is true if there is a HMC component which provides lp and epsilon
            if lp != nothing samples[i].value[:lp] = lp end
            if epsilon != nothing samples[i].value[:epsilon] = epsilon end
            if lf_num != nothing samples[i].value[:lf_num] = lf_num end
            if eval_num != nothing samples[i].value[:eval_num] = eval_num end
        end

        if PROGRESS[]
            if haskey(spl.info, :progress)
                ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
            end
        end
    end

    @info("[Gibbs] Finished with")
    @info("  Running time    = $time_total;")

    if resume_from != nothing   # concat samples
        pushfirst!(samples, resume_from.value2...)
    end
    c = Chain(0, samples)       # wrap the result by Chain

    if save_state               # save state
        save!(c, spl, model, varInfo)
    end

    return c
end


text:
include("header.jl")

@testset "getindex" begin

    for T in (Float32, Float64)
        @testset "Array" begin
            a = rand(T,3,3)
            for i in (:, 1, 1:2, 1:2:3, [1,2], [2,2], [])
                #@show i
                @test gradcheck(getindex, a, i; args=1)
                for j in (:, 1, 1:2, 1:2:3, [1,2], [2,2], [])
                    #@show i,j
                    @test gradcheck(getindex, a, i, j; args=1)
                end
            end
            j = [true,false,true]
            for i in (:, 1, 1:2, 1:2:3, [1,2], [2,2], [])
                #@show i,j
                @test gradcheck(getindex, a, i, j; args=1)
                #@show j,i
                @test gradcheck(getindex, a, j, i; args=1)
            end
            @test gradcheck(getindex, a, [1 2;1 2]; args=1)
            @test gradcheck(getindex, a, a.>0.5; args=1)
            @test gradcheck(getindex, a, CartesianIndex(1,2); args=1)
        end

        @testset "Tuple" begin
            a = Any[rand(T,3)...]
            t = tuple(a...)
            f(a,i)=(ai=a[i];s=0;for j=1:length(ai); s+=j*ai[j]; end; s)
            g = grad(f)
            for i in (1,1:2,1:2:3,[1,2],[2,2],[],[true,false,true])
                #@show i
                @test gradcheck(getindex, t, i; args=1)
                @test g(t,i)==g(a,i)==nothing || g(t,i) == tuple(g(a,i)...)
            end
        end

        @testset "Dict" begin
            g = grad(getindex)
            d = Dict(1=>rand(T), 2=>rand(T), 3=> rand(T))
            #@show d
            @test gradcheck(getindex, d, 1; args=1)
            @test collect(g(d,2)) == Any[(2=>1.0)]
        end

        @testset "size" begin
            #@show size
            f0(x)=(p=size(x); p[1]*sum(abs2,x))
            @test grad(f0)(ones(3)) == fill(6, 3)

            f1(x)=(p=size(x, 1); p*sum(abs2,x))
            @test grad(f1)(ones(3, 3)) == fill(6, 3, 3)

            # Issue #18: ambiguity error in size(rec, dims...)
            f2(x)=(p=(size(x, 1), size(x, 2)); p[1]*sum(abs2,x))
            @test grad(f2)(ones(3, 3)) == fill(6, 3, 3)
        end

        @testset "indexing" begin

            # info("Test indexing...")

            a1 = rand(T,3)
            t1 = (a1...,)
            d1 = Dict(); for i=1:length(a1); d1[i]=a1[i]; end

            s0(x)=x[1]^2+x[2]^2
            s1 = grad(s0)
            s1sum(x)=(y=s1(x);y[1]+y[2])
            s2 = grad(s1sum)

            @test gradcheck(s0,a1)
            @test gradcheck(s0,t1)
            @test gradcheck(s0,d1)
            @test gradcheck(s1sum,a1)
            @test gradcheck(s1sum,t1)
            @test gradcheck(s1sum,d1)

            f0(x)=(a=0;for i=1:length(x);a+=x[i]^2;end;a)
            f1=grad(f0)
            f1sum(x)=sum(values(f1(x)))
            f2=grad(f1sum)
            @test gradcheck(f0,a1)
            @test gradcheck(f0,t1)
            @test gradcheck(f0,d1)
            @test gradcheck(f1sum,a1)
            @test gradcheck(f1sum,t1)
            @test gradcheck(f1sum,d1)

            r0(x)=(s=0; for i=2:length(x); s+=(1-x[i-1])^2 + 2*(x[i]-x[i-1]^2)^2; end; s)
            r1 = grad(r0)
            r1sum(x)=sum(values(r1(x)))
            r2 = grad(r1sum)
            @test gradcheck(r0,a1)
            @test gradcheck(r0,t1)
            @test gradcheck(r0,d1)
            @test gradcheck(r1sum,a1)
            @test gradcheck(r1sum,t1)
            @test gradcheck(r1sum,d1)

        end

        @testset "view" begin
            a = rand(T,3,3)
            g1 = grad(x->sum(x[1:2,1:2].^2.5))
            g2 = grad(x->sum(view(x,1:2,1:2).^2.5))
            @test g1(a) == g2(a)

            g1 = grad(x->sum(x[:,1:2].^2.0))
            g2 = grad(x->sum(view(x,:,1:2).^2.0))
            g3 = grad(x->sum(selectdim(x,2,1:2).^2.0))
            @test g1(a) == g2(a) == g3(a)
        end
    end

    # Issue #73: incorrect gradient when indexing into a matrix of vectors
    a = Matrix{Array{Float32}}(undef, 2, 1)
    a[1] = [2,3,4]
    a[2] = [4,5,6]
    g1 = grad(x->x[1]' * x[2])
    g2 = grad(x->x[1, 1]' * x[2, 1])
    @test g1(a) == g2(a)

end

nothing


text:
#
# Factors Main
#
# Factors definition, constructors, and basic functions

"""
    Factor(dims::AbstractVector{Dimension}, potential::Array{Real})
    Factor(dim::Dimension, potential::Vector{Real})
    Factor(dims::AbstractVector{Symbol}, potential::Array{Real})
    Factor(dim::Symbol, potential::Vector{Real})

Create a Factor with scope `dims` and potential `potential` (or a zero
potential).
"""
type Factor{N}
    dimensions::Vector{Dimension}
    # the value mapped to each instance
    potential::Array{Float64, N}

    function Factor{D<:Dimension, K}(dimensions::Vector{D}, potential::Array{Float64, K})
        _check_dims_unique(dimensions)
        _check_dims_singleton(dimensions)
        _check_negatives(potential)

        (length(dimensions) == N) || not_enough_dims_error()
        ([size(potential)...] == [length(d) for d in dimensions]) || invalid_dim_sizes()

        (:potential in map(name, dimensions)) &&
                warn("having a dimension called `potential` will cause problems")

        return new(dimensions, potential)
    end
end

####################################################################################################
#                                   Constructors

# because of type parametrics, Vector{Dimension} ~= D<:Dimension, Vector{D}
#  the latter will work more often than not

#                                   Default Constructors
Factor{V<:Real}(dim::Dimension, potential::Vector{V}) = Factor{1}([dim], float(potential))

# convert from ints and other reals to floats
# also ,reshape
function Factor{D<:Dimension, V<:Real}(dims::Vector{D}, potential::Array{V})
    (length(dims) != 1 && ndims(potential) == 1) &&
        (potential = reshape(potential, _dims_sizes(dims)))

    return Factor{ndims(potential)}(dims, float(potential))
end

# for ranges, etc
Factor{D<:Dimension, V<:Real}(dims::Vector{D}, potential::AbstractArray{V}) =
        Factor(dims, collect(potential))

#                                   Default Value / Unitialized
"""
    Factor(dims::AbstractVector{Dims}, [value=0])

Create a factor with a potential initialized to `value`.
Use `value=nothing` for an unitialized potential
"""
function Factor{D<:Dimension}(dims::Vector{D}, value::Real=0)
    potential = fill(float(value), _dims_sizes(dims))

    return Factor(dims, potential)
end

function Factor{D<:Dimension}(dims::Vector{D}, ::Void)
    potential = Array{Float64}( _dims_sizes(dims) )

    return Factor(dims, potential)
end

Factor(dim::Dimension, value::Real=0) = Factor([dim], value)

#                                   Symbols and Potential
Factor{V<:Real}(dim::Symbol, potential::AbstractVector{V}) = Factor([dim], potential)

function Factor{V<:Real}(dims::Vector{Symbol}, potential::AbstractArray{V})
    (length(dims) == ndims(potential)) || not_enough_dims_error()
    new_dims = map(Dimension, dims, size(potential))

    return Factor(new_dims, potential)
end

#                                   Dictionary and Pairs
"""
    Factor(dims::Dict{Symbol}, [value=0])

    Factor([value=0], dims::Pair{Symbol}...)
    Factor(potential::Array{Real}, dims::Pair{Symbol}...)

Create factor from mappings from dimension names (symbols) to lengths,
ranges, or arrays.
"""
function Factor(dims::Dict{Symbol}, value=0)
    new_dims = [Dimension(name, state) for (name, state) in dims]

    return Factor(new_dims, value)
end

# this way, its order preserving, as opposed to a Dict
function Factor(potential, dims::Pair{Symbol}...)
    new_dims = [Dimension(name, state) for (name, state) in dims]

    return Factor(new_dims, potential)
end

function Factor(dims::Pair{Symbol}...)
    new_dims = [Dimension(name, state) for (name, state) in dims]

    return Factor(new_dims)
end

#                                   Zero Dimensional
"""
    Factor(potential::Real)

Create a zero-dimensional Factor.
"""
Factor(potential::Real) = Factor(Dimension[], reshape(Float64[potential], () ))

#                                   Similar
"""
    similar(ϕ::Factor)

Return an uninitialized Factor with the same dimensions as ϕ.
"""
Base.similar(ϕ::Factor) = Factor(copy(scope(ϕ)), nothing)

####################################################################################################
#                                   Methods

Base.eltype(::Factor) = Float64

"""
    names(ϕ::Factor)

Get an array of the names (symbols) of the dimensions of `ϕ`.
"""
Base.names(ϕ::Factor) = map(name, ϕ.dimensions)

"""
    potential(ϕ::Factor)

Get the potential of a Factor.
"""
potential(ϕ::Factor) = ϕ.potential

"""
    scope(ϕ)

Get the dimensions of `ϕ`.
"""
scope(ϕ::Factor) = ϕ.dimensions

Base.ndims{N}(::Factor{N}) = N

"""
    size(ϕ, [dims...])

Return a tuple of the dimensions of `ϕ`.
"""
Base.size(ϕ::Factor) = size(ϕ.potential)
Base.size(ϕ::Factor, dim::Symbol) = size(ϕ.potential, indexin(dim, ϕ))
Base.size{N}(ϕ::Factor, dims::Vararg{Symbol, N}) = ntuple(k -> size(ϕ, dims[k]), Val{N})

"""
    length(ϕ)

Return thelength of `ϕ`.
"""
Base.length(ϕ::Factor) = length(ϕ.potential)

"""
    values(ϕ)

Return the potential of `ϕ`.
"""
Base.values(ϕ::Factor) = ϕ.potential

"""
    in(d, ϕ::Factor)

Check if `Dimension` (or `Symbol`) `d` is in `ϕ`.
"""
Base.in(d::Dimension, ϕ::Factor) = d in scope(ϕ)
Base.in(d::Symbol, ϕ::Factor) = d in names(ϕ)

"""
    indexin(dims, ϕ::Factor)

Find the index of `dims` in `ϕ`. Return 0 if not in `ϕ`.
"""
Base.indexin(dim::Symbol, ϕ::Factor) = findnext(names(ϕ), dim, 1)
Base.indexin(dims::Vector{Symbol}, ϕ::Factor) = indexin(dims, names(ϕ))

Base.indexin(dim::Dimension, ϕ::Factor) = findnext(scope(ϕ), dim, 1)
Base.indexin{D<:Dimension}(dims::AbstractVector{D}, ϕ::Factor) = indexin(dims, scope(ϕ))

# Julia has a hashing issue (?), so if 2 dims have the same name & support, Julia
#  may not find find the first via indexin ...
# Base.indexin(dim::Dimension, ϕ::Factor) = findnext(scope(ϕ), dim, 1)
# Base.indexin{D<:Dimension}(dims::AbstractVector{D}, ϕ::Factor) = indexin(dims, scope(ϕ))

"""
    getdim(ϕ::Factor, dim::Symbol)
    getdim(ϕ::Factor, dims::Vector{Symbol})

Get the dimension with name `dim` in `ϕ`.
"""
@inline function getdim(ϕ::Factor, dim::Symbol)
    i = indexin(dim, ϕ)

    (i == 0) && not_in_factor_error(dim)

    return ϕ.dimensions[i]
end

@inline function getdim(ϕ::Factor, dims::Vector{Symbol})
    _check_dims_valid(dims, ϕ)

    inds = indexin(dims, ϕ)
    return ϕ.dimensions[inds]
end

@inline function Base.permutedims!(ϕ::Factor, perm)
    ϕ.potential = permutedims(ϕ.potential, perm)
    ϕ.dimensions = ϕ.dimensions[perm]

    return ϕ
end

Base.permutedims(ϕ::Factor, perm) = permutedims!(deepcopy(ϕ), perm)

"""
    push(ϕ::Factor, dim::Dimension)

Return a new Factor with `dim` appended to it.
"""
@inline function push(ϕ::Factor, dim::Dimension)
    (name(dim) in names(ϕ)) && error("Dimension $(name(dim)) already exists")

    new_dims = [scope(ϕ)..., dim]
    new_pot = duplicate(ϕ.potential, length(dim))

    return Factor(new_dims, new_pot)
end


text:
export PetscVec, PetscVecSetType, PetscVecSetValues, PetscVecAssemblyBegin, PetscVecAssemblyEnd, PetscVecSetSizes, PetscVecGetSize, PetscVecNorm, PetscVecGetValues, PetscVecGetOwnershipRange, PetscVecGetArray, PetscVecRestoreArray, PetscVecGetArrayRead, PetscVecRestoreArrayRead, PetscVecSet, PetscVecSqrtAbs, PetscVecLog, PetscVecExp, PetscVecAbs, PetscVecMax, PetscVecMin, PetscVecCopy, PetscVecDuplicate, PetscVecAXPY, PetscVecAXPBY, PetscVecAYPX, PetscVecWAXPY, PetscVecMAXPY, PetscVecAXPBYPCZ, PetscVecScale, PetscVecDot, PetscVecTDot, PetscVecSum, PetscVecSwap, PetscVecReciprocal, PetscVecShift, PetscVecPointwiseMult, PetscVecPointwiseDivide, set_values1!, get_values1!

type PetscVec
  pobj::Ptr{Void}
  function PetscVec(comm::MPI_Comm)
#    comm = PETSC_COMM_SELF();
    vec = Array(Ptr{Void},1)
    err = ccall(( :VecCreate, libpetsclocation ), PetscErrorCode,(comm_type,Ptr{Void}),comm,vec);
    vec = new(vec[1])
#    finalizer(vec,PetscDestroy)
    # does not seem to be called immediately when vec is no longer visible, is it called later during garbage collection? - yes
    return vec

  end

  function PetscVec(pobj::Ptr{Void})  # default constructor
    return new(pobj)
  end

end

  function PetscDestroy(vec::PetscVec)
    if (vec.pobj != 0)
      err = ccall(( :VecDestroy, libpetsclocation), PetscErrorCode, (Ptr{Ptr{Void}},), &vec.pobj);
    end
    vec.pobj = 0  # unnecessary? vec no longer has any references to it
#    println("VecDestroy called")
  end

  function PetscVecSetType(vec::PetscVec,name)
    err = ccall((:VecSetType,  libpetsclocation), PetscErrorCode, (Ptr{Void}, Cstring), vec.pobj,name);
  end

  function PetscVec(array::Array{PetscScalar})
    vec = PetscVec()
    err = ccall(( :VecSetType,  libpetsclocation), PetscErrorCode,(Ptr{Void},Cstring), vec.pobj,"seq");
    err = ccall( (:VecSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, PetscInt), vec.pobj,length(array),length(array));
    # want a PetscInt array so build it ourselves
    idx = Array(PetscInt, length(array));
    for i=1:length(array);  idx[i] = i-1;  end
    err = ccall( ( :VecSetValues,  libpetsclocation), PetscErrorCode,(Ptr{Void},PetscInt, Ptr{PetscInt},Ptr{PetscScalar},Int32), vec.pobj,length(idx),idx,array,PETSC_INSERT_VALUES);
    err = ccall( ( :VecAssemblyBegin,  libpetsclocation), PetscErrorCode,(Ptr{Void},), vec.pobj);
    err = ccall( ( :VecAssemblyEnd,  libpetsclocation), PetscErrorCode, (Ptr{Void},), vec.pobj);
    return vec
  end

  function PetscVecSetValues(vec::PetscVec,idx::Array{PetscInt},array::Array{PetscScalar},flag::Integer)

    err = ccall( ( :VecSetValues,  libpetsclocation), PetscErrorCode,
                 (Ptr{Void},PetscInt ,Ptr{PetscInt},Ptr{PetscScalar},Int32),
                 vec.pobj,length(idx), idx,array,flag);
    return err
  end

  function PetscVecSetValues(vec::PetscVec,idx::Array{PetscInt},array::Array{PetscScalar})
    PetscVecSetValues(vec,idx,array,PETSC_INSERT_VALUES)
  end

  function PetscVecSetValues(vec::PetscVec,array::Array{PetscScalar})
    idx = Array(PetscInt,length(array))
    for i=1:length(array);  idx[i] = i-1;  end
    PetscVecSetValues(vec,idx,array,PETSC_INSERT_VALUES)
  end

  # 1-based indexing that unifies regular matrices and Petsc matrices
  #----------------------------------------------------------------------------
  function set_values1!(vec::PetscVec, idx::Array{PetscInt}, vals::Array{PetscScalar}, flag::Integer=PETSC_INSERT_VALUES)
    for i=1:length(idx)
      idx[i] -= 1
    end

    err = PetscVecSetValues(vec, idx, vals, flag)

    for i=1:length(idx)
      idx[i] += 1
    end

    return err
  end

  function set_values1!{T}(vec::AbstractVector, idx::Array{PetscInt}, vals::Array{T}, flag::Integer=PETSC_INSERT_VALUES)
    if flag == PETSC_INSERT_VALUES
      for i in idx
        vec[i] = vals[i]
      end
    elseif flag == PETSC_ADD_VALUES
      for i in idx
        vec[i] += vals[i]
      end
    end

    return PetscErrorCode(0)
  end

  function get_values1!(vec::PetscVec, idx::Array{PetscInt}, vals::Array{PetscScalar})
    for i=1:length(idx)
      idx[i] -= 1
    end

    err = PetscVecGetValues(vec, idx, vals)

    for i=1:length(idx)
      idx[i] += 1
    end

    return err
  end

  function get_values1!(vec::AbstractVector, idx::Array{PetscInt}, vals::Array)
    for i=1:length(idx)
      vals[i] = vec[idx[i]]
    end

    return PetscErrorCode(0)
  end

  function PetscVecAssemblyBegin(obj::PetscVec)
    err = ccall( ( :VecAssemblyBegin,  libpetsclocation), PetscErrorCode, (Ptr{Void},), obj.pobj);
  end

  function PetscVecAssemblyEnd(obj::PetscVec)
    err = ccall( ( :VecAssemblyEnd,  libpetsclocation), PetscErrorCode,(Ptr{Void},), obj.pobj);
  end

  function PetscVecSetSizes(vec::PetscVec,n::PetscInt, N::PetscInt)
    err = ccall( ( :VecSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, PetscInt), vec.pobj,n,N);
  end

  function PetscView(obj::PetscVec,viewer)
    err = ccall( ( :VecView,  libpetsclocation), PetscErrorCode, (Ptr{Void},Int64),obj.pobj,0);
  end

  function PetscVecGetSize(obj::PetscVec)
    n = Array(PetscInt, 1)
    err = ccall( ( :VecGetSize,  libpetsclocation), PetscErrorCode, (Ptr{Void},Ptr{PetscInt}), obj.pobj,n);
    return n[1]
  end

  function PetscVecNorm(obj::PetscVec,normtype::Integer)
    n = Array(PetscReal,1)
    err = ccall( ( :VecNorm,  libpetsclocation), PetscScalar, (Ptr{Void},Int32,Ptr{PetscReal}), obj.pobj,normtype, n);
    return n[1]
  end

  function PetscVecNorm(obj::PetscVec)
    return PetscVecNorm(obj,PETSC_NORM_2)
  end

function PetscVecGetValues(vec::PetscVec, ni::Integer, ix::AbstractArray{PetscInt,1}, y::AbstractArray{PetscScalar,1})

     # need indices to be PetscInt
#     ix_local = Array(PetscInt, ni)
#     for i=1:ni
#       ix_local[i] = ix[i]
#     end

    err = ccall((:VecGetValues,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}),vec.pobj, ni, ix, y)

    return nothing
end

function PetscVecGetValues(vec::PetscVec,idx::AbstractArray{PetscInt,1}, y::AbstractArray{PetscScalar, 1})

  PetscVecGetValues(vec, length(idx), idx, y)
end

function PetscVecGetOwnershipRange(vec::PetscVec)
    low = Array(PetscInt, 1)
    high = Array(PetscInt, 1)

    ccall((:VecGetOwnershipRange,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),vec.pobj, low, high)

  return low[1], high[1]
end

# new functions

function PetscVecGetArray(vec::PetscVec)
# gets a pointer to the data underlying a Petsc vec, turns it into a Julia
# array
# ptr_arr must be passed into PetscVecRestoreArray

#TODO: use pointer(arr) to avoid returning ptr_arr

    ptr_arr = Array(Ptr{PetscScalar}, 1)
    ccall((:VecGetArray,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)

    first, last = PetscVecGetOwnershipRange(vec)
    len = last - first
    arr = pointer_to_array(ptr_arr[1], len)
    return arr, ptr_arr
end

function PetscVecRestoreArray(vec::PetscVec, ptr_arr::Array{Ptr{PetscScalar}, 1})
    ccall((:VecRestoreArray,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)
end

function PetscVecGetArrayRead(vec::PetscVec)
    ptr_arr = Array(Ptr{PetscScalar}, 1)
    ccall((:VecGetArrayRead,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)

    first, last = PetscVecGetOwnershipRange(vec)
    len = last - first
    arr = pointer_to_array(ptr_arr[1], len)
    return arr, ptr_arr
end

function PetscVecRestoreArrayRead(vec::PetscVec, ptr_arr::Array{Ptr{PetscScalar}, 1})

    ccall((:VecRestoreArrayRead,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)
#    ccall((:VecRestoreArrayRead,petsc),PetscErrorCode,(Vec,Ptr{Ptr{PetscScalar}}),arg1,arg2)
end

function PetscVecSet(vec::PetscVec, val::PetscScalar)
    ccall((:VecSet,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, val)
end

function PetscVecSqrtAbs(vec::PetscVec)
    ccall((:VecSqrtAbs,petsc),PetscErrorCode,(Ptr{Void},), vec.pobj)
end

function PetscVecLog(vec::PetscVec)
    ccall((:VecLog,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

function PetscVecExp(vec::PetscVec)
    ccall((:VecExp,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

function PetscVecAbs(vec::PetscVec)
    ccall((:VecAbs,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

function PetscVecMax(vec::PetscVec)
    r = Array(PetscReal, 1) # max value
    idx = Array(PetscInt, 1)  # index of max value
    ccall((:VecMax,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscInt},Ptr{PetscReal}), vec.pobj, idx, r)

    return r[1], idx[1]
end

function PetscVecMin(vec::PetscVec)
    r = Array(PetscReal, 1) # min value
    idx = Array(PetscInt, 1)  # index of min value

    ccall((:VecMin,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscReal}), vec.pobj, idx, r)

    return r[1], idx[1]
end

function PetscVecReciprocal(vec::PetscVec)
    ccall((:VecReciprocal,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

function PetscVecShift(vec::PetscVec, a::PetscScalar)
    ccall((:VecShift,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, a)
end

function PetscVecPointwiseMult(w::PetscVec, x::PetscVec,y::PetscVec)
    ccall((:VecPointwiseMult,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, Ptr{Void}), w.pobj, x.pobj, y.pobj)
end

function PetscVecPointwiseDivide(w::PetscVec, x::PetscVec, y::PetscVec)
    ccall((:VecPointwiseDivide,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, Ptr{Void}), w.pobj, x.pobj, y.pobj)
end

function PetscVecCopy(vec::PetscVec , vec2::PetscVec)
    ccall((:VecCopy,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}), vec.pobj, vec2.pobj)
end

function PetscVecDuplicate( vec::PetscVec)
    ptr_arr = Array(Ptr{Void}, 1)

    ccall((:VecDuplicate,petsc),PetscErrorCode,( Ptr{Void}, Ptr{Ptr{Void}}), vec.pobj, ptr_arr)

    return PetscVec(ptr_arr[1])
end

# Some vector linear algebra
function PetscVecAXPY( vec1::PetscVec, a::PetscScalar, vec2::PetscVec)
    ccall((:VecAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}), vec1.pobj, a, vec2.pobj)
end

function PetscVecAXPBY(vec1::PetscVec, a::PetscScalar, b::PetscScalar, vec2::PetscVec)
    ccall((:VecAXPBY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,PetscScalar,Ptr{Void}),vec1.pobj, a, b, vec2.pobj)
end

function PetscVecMAXPY(vec1::PetscVec, n::PetscInt, a::AbstractArray{PetscScalar, 1}, x::AbstractArray{Ptr{Void}, 1})
# the vector x must contains the pointers from the PetscVec objects, not the PetscVec objects themselves

    ccall((:VecMAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscScalar},Ptr{Ptr{Void}}),vec1.pobj, n, a, x)
end

function PetscVecAYPX(vec1::PetscVec, a::PetscScalar, vec2::PetscVec)
    ccall((:VecAYPX,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}),vec1.pobj, a, vec2.pobj)
end

function PetscVecWAXPY(w::PetscVec, a::PetscScalar, x::PetscVec, y::PetscVec)
    ccall((:VecWAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void},Ptr{Void}), w.pobj, a, x.pobj, y.pobj)
end

function PetscVecAXPBYPCZ(z::PetscVec, alpha::PetscScalar, beta::PetscScalar, gamma::PetscScalar, x::PetscVec, y::PetscVec)
    ccall((:VecAXPBYPCZ,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,PetscScalar,PetscScalar,Ptr{Void},Ptr{Void}), z.pobj, alpha,beta, gamma, x.pobj, y.pobj)
end

function PetscVecScale(vec::PetscVec, a::PetscScalar)
    ccall((:VecScale,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, a)
end

function PetscVecDot(x::PetscVec, y::PetscVec)
    r = Array(PetscScalar, 1)
    err = ccall((:VecDot,petsc),PetscErrorCode,( Ptr{Void}, Ptr{Void}, Ptr{PetscScalar}), x.pobj, y.pobj, r)

    println("PetscVecDot error code = ", err)
    return r[1]
end

function PetscVecTDot(x::PetscVec, y::PetscVec)
    r = Array(PetscScalar, 1)
    ccall((:VecTDot,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void},Ptr{PetscScalar}), x.pobj, y.pobj, r)

    return r[1]
end

function PetscVecSum(vec::PetscVec)
    r = Array(PetscScalar, 1)
    ccall((:VecSum,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscScalar}),vec.pobj, r)
    return r[1]
end

function PetscVecSwap(x::PetscVec, y::PetscVec)
    ccall((:VecSwap,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}), x.pobj, y.pobj)

end


text:
__precompile__(false)

module JAXTAM

using DataFrames
using DelimitedFiles
using Dates
using FileIO
using JSON
using JLD2
using FITSIO
using Arrow
using Feather
using DataStructures

using FFTW
using Statistics
using StatsBase
using OnlineStats
using LombScargle
using DSP
using LinearAlgebra

using Hyperscript

using Measures
using Plots
gr()

# Fix loading JLD2 files not finding types in workspace
export DataFrames
export Dates
export Arrow

abstract type JAXTAMData end

# @__DIR__ returns the location of this file
const __sourcedir__ = abspath(@__DIR__, "..")

include("missions/mission_dispatch.jl")
include("utils/logging.jl")
include("utils/errors.jl")
include("io/master_tables.jl")
include("io/master_query.jl")
include("io/misc.jl")
include("io/data_download.jl")
include("science/read_events.jl")
include("science/calibrate.jl")
include("science/lcurve.jl")
include("science/gtis.jl")
include("science/fspec.jl")
include("science/pgram.jl")
include("science/plots.jl")
include("science/sgram.jl")
include("science/espec.jl")
include("web/webgen.jl")
include("web/subgen.jl")
include("web/reports.jl")
include("web/autoreport.jl")

end


text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress

cd(@__DIR__)
#@load "../nn_scalar.jld2" nn
@load "/home2/vavrines/Coding/Flowmachine/src/1d/sampler/nn_rif.jld2" nn

D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "sod"
    D[:space] = "1d1f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "fix"
    D[:cfl] = 0.3
    D[:maxTime] = 0.1

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 100
    D[:vMeshType] = "rectangle"
    D[:nug] = 0

    D[:knudsen] = 1e-4
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 0.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

ks = SolverSet(D)
ctr, face = init_fvm(ks, ks.ps, :static_array; structarray = true)
#=for i in eachindex(ctr)
    prim = [2.0 * rand(), 0.0, 1 / rand()]

    ctr[i].prim .= prim
    ctr[i].w .= prim_conserve(prim, ks.gas.γ)
    ctr[i].f .= maxwellian(ks.vs.u, prim)
end=#
for i in eachindex(face)
    face[i].fw .= 0.0
    face[i].ff .= 0.0
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ctr[1].w)
for iter = 1:nt
    println("iteration: $iter")

    #reconstruct!(ks, ctr)

    #evolve!(ks, ctr, face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))
    for i in eachindex(face)
        w = (ctr[i-1].w .+ ctr[i].w) ./ 2
        prim = (ctr[i-1].prim .+ ctr[i].prim) ./ 2
        sw = (ctr[i].w .- ctr[i-1].w) / ks.ps.dx[i]
        tau = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)

        #regime = nn([w; sw; tau]) |> onecold
        regime = 1 - Int(round(nn([w; sw; tau])[1]))

        if regime == 1
            flux_gks!(
                face[i].fw,
                face[i].ff,
                ctr[i-1].w .+ ctr[i-1].sw .* ks.ps.dx[i-1] / 2,
                ctr[i].w .- ctr[i].sw .* ks.ps.dx[i] / 2,
                ks.vs.u,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1] / 2,
                ks.ps.dx[i] / 2,
                ctr[i-1].sw,
                ctr[i].sw,
            )
        else
            flux_kfvs!(
                face[i].fw,
                face[i].ff,
                ctr[i-1].f .+ ctr[i-1].sf .* ks.ps.dx[i-1] / 2,
                ctr[i].f .- ctr[i].sf .* ks.ps.dx[i] / 2,
                ks.vs.u,
                ks.vs.weights,
                dt,
                ctr[i-1].sf,
                ctr[i].sf,
            )
        end
    end

    KitBase.update!(ks, ctr, face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

    t += dt
    #if t > ks.set.maxTime || maximum(res) < 5.e-7
    #    break
    #end
end

plot(ks, ctr)

# test
regime = zeros(Int, ks.ps.nx)
for i = 1:ks.ps.nx
    sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
    tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
    x, y = regime_data(ks, ctr[i].w, sw, ctr[i].f)

    #regime[i] = nn(x) |> onecold
    regime[i] = round(nn(x)[1])
end
plot!(ks.ps.x[1:ks.ps.nx], regime)

regime1 = zeros(Int, ks.ps.nx)
KnGLL = zeros(ks.ps.nx)
for i = 1:ks.ps.nx
    sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
    L = abs(ctr[i].w[1] / sw[1])
    ℓ = (1/ctr[i].prim[end])^ks.gas.ω / ctr[i].prim[1] * sqrt(ctr[i].prim[end]) * ks.gas.Kn

    KnGLL[i] = ℓ / L
    regime1[i] = ifelse(KnGLL[i] >= 0.05, 1, 0)
end
plot!(ks.ps.x[1:ks.ps.nx], regime1)

plot(ks.ps.x[1:ks.ps.nx], KnGLL)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
if abspath(PROGRAM_FILE) == @__FILE__
    value = join((string(x) for x in 0:9 if x > 4 && x < 8), " ")
    @assert(value == "5 6 7")
    value2 = join((string(x) * " " * string(y) for x in 0:4 if x > 1 for y = 0:1), " ")
    @assert(value2 == "2 0 2 1 3 0 3 1 4 0 4 1")
    value3 = join((string(x + y) for x in 0:4 if x > 1 for y = 0:3), " ")
    @assert(value3 == "2 3 4 5 3 4 5 6 4 5 6 7")
end


text:
using DifferentialEquations
using Plots
plotlyjs()

#constants
m = 1.0                             # mass[Kg]
g = 9.81                            # gravitational acceleration [m/s²]

#function to compute the ODEs
function pendulum!(du,u,p,t)
    r,v,θ,ω=u
    du[1]=dr=v
    du[2]=dv=r*ω^2+g*cos(θ)
    du[3]=dθ=ω
    du[4] =dω=-(g\r)*sin(θ) - (2*v*ω)/r
end

#initial conditions
r₀=4                            #initial radial distance
v₀=-0.1                          #initial radial velocity
θ₀ = 1π/18                        # initial angular deflection [rad]
ω₀ = 1                       # initial angular velocity [rad/s]
u₀ = [r₀,v₀,θ₀, ω₀]                   # initial state vector
tspan = (0.0,5.0)             # time interval

#defining the problem and solving
prob = ODEProblem(pendulum!,u₀,tspan)
sol = solve(prob,saveat=0.1)

#plotting
p1=plot(sol,vars=(0,1),label="r",
    color="purple", xaxis="t (s)", yaxis="r (m)")
p2=plot(sol,vars=(0,2),label="v",
    color="blue", xaxis="t (s)", yaxis="v (m/s)")
p3=plot(sol,vars=(0,3),label="θ",
    color="green", xaxis="t (s)", yaxis="θ (rad)")
p4=plot(sol,vars=(0,4),label="ω",
    color="orange", xaxis="t (s)", yaxis="ω (rad/s)")
display(plot(p1,p2,p3,p4,layout=(2,2)))
savefig("./plots/simple length changing pendulum-plot.png")

#trajectory
t=sol.t
U=sol[1:end,:]
(r,v,θ,ω)=[U[x,:] for x in 1:size(U,1)]
x=r.*sin.(θ)
y=zeros(size(t))
z=r.*cos.(θ)
display(plot(x,y,z,linewidth=2,label="trajectory",
    xaxis="x (m)",yaxis="y (m)",zaxis="z (m)",grid=(:on,:black)))
savefig("./plots/simple length changing pendulum-trajectory.png")

max_stars_repo_path: Modeling the Dragon Module as a Pendulum/Scripts/length changing pendulum.jl


text:
include("src/analysis.jl")
# include("src/experiments.jl")

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
#TODO should this be <: AbstractArray{T,4}?
"""
    AbstractStorage{T,scalar}

Abstract supertype for storage of different Keldysh components of nonequilibrium Green's functions taking into account symmetries.
Represents an array with dimensions (norb,norb,N,M).
Supports basic arithmetic and indexing.
`scalar=true` statically expresses norb=1.
"""
abstract type AbstractStorage{T, scalar} end

Base.eltype(::Type{<:AbstractStorage{T,scalar}}) where {T, scalar} = T
Base.eltype(X::AbstractStorage) = eltype(typeof(X))

Base.getindex(X::AbstractStorage{T,true}, i, j) where {T} = X[1,1,i,j]
Base.getindex(X::AbstractStorage{T,false}, i, j) where {T} = X[:,:,i,j]
Base.setindex!(X::AbstractStorage{T,true}, v, i, j) where {T} = X[1,1,i,j] = v
Base.setindex!(X::AbstractStorage{T,false}, v, i, j) where {T} = X[:,:,i,j] = v

function Base.size(X::AbstractStorage)
  return (X.norb, X.norb, X.N, X.M)
end

function Base.similar(X::T) where T <: AbstractStorage
  T(similar(X.data))
end

function Base.zero(X::T) where T <: AbstractStorage
  T(zero(X.data))
end

function Base.:+(X::T, Y::T) where T <: AbstractStorage
  @assert size(X) == size(Y)
  return T(X.data .+ Y.data)
end

function Base.:-(X::T, Y::T) where T <: AbstractStorage
  @assert size(X) == size(Y)
  return T(X.data .- Y.data)
end

function Base.:*(X::T, α::Number) where T <: AbstractStorage
  return T(X.data * α)
end

function Base.:*(α::Number, X::T) where T <: AbstractStorage
  return X * α
end

include("storage/generic.jl")
include("storage/herm.jl")
include("storage/herm_toeplitz.jl")
include("storage/periodic.jl")


text:
export isfqdn

function checkparts(part::AbstractString, allowUnderscores::Bool=false)
    if length(part) > 63
        return false
    end

    if !contains(part, r"^[a-z_\u00a1-\uffff0-9-]+$"i)
        return false
    end

    # disallow full-width chars
    if contains(part, r"[\uff01-\uff5e]")
        return false
    end

    # disallow parts starting or ending with hyphen
    if contains(part, r"^-|-$")
        return false
    end

    if !allowUnderscores && contains(part, r"_")
        return false
    end

    return true
end

function isfqdn(str::AbstractString, requireTld::Bool=true, allowUnderscores::Bool=false, allowTrailingDot::Bool=false, allowNumericTld::Bool=false, allowWildcard::Bool=false)::Bool
    if allowTrailingDot && str[end] === '.'
        str = str[1:end-1]
    end

    if allowWildcard && findfirst("*.", str)[1] === 1
        str = str[3:end]
    end

    parts = split(str, ".")
    tld = parts[end]

    if requireTld
        if length(parts) < 2
            return false
        end

        tldReg = r"^([a-z\u00A1-\u00A8\u00AA-\uD7FF\uF900-\uFDCF\uFDF0-\uFFEF]{2,}|xn[a-z0-9-]{2,})$"i
        if !contains(tld, tldReg)
            return false
        end

        # disallow whitespaces
        if contains(tld, r"\s"i)
            return false
        end
    end

    # reject numeric TLDs
    if !allowNumericTld && contains(tld, r"^\d+$")
        return false
    end

    return all(map(part -> checkparts(part, allowUnderscores), parts))
end


text:
# ---------------------------------------------------------------
# take the content of a `file` and print it into a set of `ios`
function tee_file(ios::Vector, file::AbstractString;
        print_fun = print,
        finish_time::Ref{Float64} = Ref{Float64}(time() + 1e20),
        buffsize = 150,
        upfrec = 1.0,
        printlk::ReentrantLock = ReentrantLock(),
        append = false
    )

    # check
    isempty(ios) && return

    # empty
    !append && _empty!.(ios)

    buff = Vector{Char}(undef, buffsize)
    bi = 0
    open(file, "r") do fileio
        while true
            while eof(fileio) # will wait at the end of the file
                if (bi > 0) # will print any new content
                    str = join(buff[1:bi])
                    lock(printlk) do
                        for io in ios
                            print_fun(io, str)
                        end
                    end
                    bi = 0
                end
                (time() > finish_time[]) && return # time out
                sleep(upfrec)
            end

            ch = read(fileio, Char)
            bi += 1
            buff[bi] = ch
            if (bi == buffsize) # will print the whole buff
                str = join(buff)
                lock(printlk) do
                    for io in ios
                        print_fun(io, str)
                    end
                end
                bi = 0
            end
        end
    end
end


text:
@testset "memory" begin

let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# Mem.free()
    #=@test b == =# Mem.total()
    #=@test b-a == =# Mem.used()
end

# pointer-based
for async in [false, true]
    src = 42

    buf1 = Mem.alloc(sizeof(src))

    Mem.set!(buf1, UInt32(0), sizeof(Int)÷sizeof(UInt32); async=async)

    Mem.upload!(buf1, Ref(src), sizeof(src); async=async)

    dst1 = Ref(0)
    Mem.download!(dst1, buf1, sizeof(src); async=async)
    async && synchronize()
    @test src == dst1[]

    buf2 = Mem.alloc(sizeof(src))

    Mem.transfer!(buf2, buf1, sizeof(src); async=async)

    dst2 = Ref(0)
    Mem.download!(dst2, buf2, sizeof(src); async=async)
    async && synchronize()
    @test src == dst2[]

    Mem.free(buf2)
    Mem.free(buf1)
end

# array-based
for async in [false, true]
    src = [42]

    buf1 = Mem.alloc(src)

    Mem.upload!(buf1, src; async=async)

    dst1 = similar(src)
    Mem.download!(dst1, buf1; async=async)
    async && synchronize()
    @test src == dst1

    buf2 = Mem.upload(src)

    dst2 = similar(src)
    Mem.download!(dst2, buf2; async=async)
    async && synchronize()
    @test src == dst2

    Mem.free(buf1)
end

# type-based
for async in [false, true]
    buf = Mem.alloc(Int)

    # there's no type-based upload, duh
    src = [42]
    Mem.upload!(buf, src)

    dst = Mem.download(eltype(src), buf; async=async)
    async && synchronize()
    @test src == dst
end

let
    @test_throws ArgumentError Mem.alloc(Function, 1)   # abstract
    @test_throws ArgumentError Mem.alloc(Array{Int}, 1) # UnionAll
    @test_throws ArgumentError Mem.alloc(Integer, 1)    # abstract
    # TODO: can we test for the third case?
    #       !abstract && leaftype seems to imply UnionAll nowadays...

    # zero-width allocations should be permitted
    null = Mem.alloc(Int, 0)
    Mem.free(null)

    # double-free should throw
    x = Mem.alloc(1)
    Mem.free(x)
    @test_throws_cuerror CUDAdrv.ERROR_INVALID_VALUE Mem.free(x)
end

let
    @eval mutable struct MutablePtrFree
        foo::Int
        bar::Int
    end
    buf = Mem.alloc(MutablePtrFree)
    Mem.upload!(buf, [MutablePtrFree(0,0)])
    Mem.free(buf)
end

let
    @eval mutable struct MutableNonPtrFree
        foo::Int
        bar::String
    end
    @test_throws ArgumentError Mem.alloc(MutableNonPtrFree)
end

end


text:
using ElectricFields
using Test
using Unitful
using UnitfulAtomic
using IntervalSets
using FFTW

@testset "ElectricFields.jl" begin
    include("namespace_macro.jl")
    include("units.jl")
    include("rotations.jl")
    include("field_creation.jl")
    include("gaussian_pulses.jl")
    include("envelopes.jl")
    include("carriers.jl")
    include("field_types.jl")
    include("arithmetic.jl")

    @testset "Time axis" begin
        @field(F) do
            I₀ = 1.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
        end

        @test ElectricFields.default_sampling_frequency(F) == 100/2

        @test steps(F) == 600
        @test steps(F, 300) == 1800

        @test timeaxis(F) == range(-6, stop=6, length=600)
        @test timeaxis(F, 300) == range(-6, stop=6, length=1800)
        @test timeaxis(F, 0) == -6:12:-6
        @test timeaxis(F, (26,0.5)) == range(-6, length=26, step=0.5)
        @test timeaxis(F, (24,0.5)) == range(-6, length=24, step=0.5)
    end

    @testset "Strong-field properties" begin
        @field(F) do
            I₀ = 1.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
        end

        @test austrip(ponderomotive_potential(F)) ≈ 1/(4*π^2) rtol=1e-7
        @test austrip(ponderomotive_potential(F+F)) ≈ 2/(4*π^2) rtol=1e-7

        @test keldysh(F, 2u"Eh_au") ≈ 2π rtol=1e-7

        @test free_oscillation_amplitude(F) ≈ 1/π^2 rtol=1e-7
        @test free_oscillation_amplitude(delay(F, 100.0)) ≈ 1/π^2 rtol=1e-7
        @test free_oscillation_amplitude(F+F) ≈ 2/π^2 rtol=1e-7
    end

    include("analytic_continuation.jl")

    # include("sellmeier.jl")
end


text:
module Var

using Random: AbstractRNG
using ..Basis

"""
`recurse(f, ω)`
Recursively apply contextual execution to internals of `f`"""
function recurse end

include("variable.jl")          # Random / Parametric Variables
include("member.jl")            # Families
include("primparam.jl")         # Primitive Parameters
include("distributions.jl")     # Primitive Distributions
include("multivariate.jl")      # Multivariate Distributions
include("constant.jl")          # Constant distribution
include("pointwise.jl")         # Point wise variable application
include("dispatch.jl")          # Contextual application

end


text:
module FineShiftBenchmarks
using Printf
using InterpolationKernels
using BenchmarkTools
using FineShift
using FineShift.Impl: Adjoint, Direct
using FineShift.Impl: REFERENCE, JUMP, SPLIT, LOCALVALUES, SIMD
const PLOTTING = false
if PLOTTING
    using PyPlot, YPlot
    const plt = PyPlot
end
include("alt.jl")

function testalternatives()
    T = Float64
    ker = CatmullRomSpline(T)
    off = 1.7
    off,wgt = FineShift.shiftparams(ker,off)
    n1, n2 = 32, 33
    n1p, n2p = 35, 36
    arr = convert(Array{T,2}, sin.(.3 .* (1:n1)).*cos.(.4 .* (1:n2))');

    z1 = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false,Val(REFERENCE));
    z2 = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false,Val(REFERENCE));

    z1a = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false);
    z2a = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false);
    println("\nVal(0): ",extrema(z1a - z1)," ",extrema(z2a - z2))

    z1b = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false,Val(0));
    z2b = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false,Val(0));
    println("\nVal(0): ",extrema(z1b - z1)," ",extrema(z2b - z2))

    z1c = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false,Val(JUMP));
    z2c = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false,Val(JUMP));
    println("\nVal(JUMP): ",extrema(z1c - z1)," ",extrema(z2c - z2))

    z1d = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false,Val(LOCALVALUES));
    z2d = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false,Val(LOCALVALUES));
    println("\nVal(LOCALVALUES): ",extrema(z1d - z1)," ",extrema(z2d - z2))

    z1e = correlate!(Array{T,2}(undef,n1p,n2),arr,off,wgt,1,false,Val(SPLIT));
    z2e = correlate!(Array{T,2}(undef,n1,n2p),arr,off,wgt,2,false,Val(SPLIT));
    println("\nVal(SPLIT): ",extrema(z1e - z1)," ",extrema(z2e - z2))

    wrk1 = Array{T,2}(undef,n1p,n2)
    wrk2 = Array{T,2}(undef,n1,n2p)

    println("\nPackage version:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false);
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false);
    println("\nReference version:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false,Val(REFERENCE));
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false,Val(REFERENCE));
    println("\nDefault version:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false,Val(0));
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false,Val(0));
    println("\nDefault version with long jumps:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false,Val(JUMP));
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false,Val(JUMP));
    println("\nDefault version with local values:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false,Val(LOCALVALUES));
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false,Val(LOCALVALUES));
    println("\nSplit version:")
    @btime correlate!($wrk1,$arr,$off,$wgt,1,false,Val(SPLIT));
    @btime correlate!($wrk2,$arr,$off,$wgt,2,false,Val(SPLIT));
end

function test1d(T::Type{<:AbstractFloat}=Float64, len::Int=1000)
    off = T(-0.7)
    ker = CatmullRomSpline(T)
    l, w = FineShift.shiftparams(ker, off)

    x1 = rand(T,len)
    x2 = copy(x1)
    y1 = Array{T}(undef,len)
    y2 = Array{T}(undef,len)
    fineshift!(y1, x1, ker, off, 1, false)
    correlate!(y2, x2, l, w, 1, false)
    @printf("max. abs. dif. direct: %g\n", maximum(abs.(y2 - y1)))

    y1 = rand(T,len)
    y2 = copy(y1)
    x1 = Array{T}(undef,len)
    x2 = Array{T}(undef,len)
    fineshift!(x1, y1, ker, off, 1, true)
    correlate!(x2, y2, l, w, 1, true)
    @printf("max. abs. dif. direct: %g\n", maximum(abs.(x2 - x1)))

    @printf("direct  1D")
    @btime fineshift!($y1, $x1, $ker, $off, 1, false)

    @printf("adjoint 1D")
    @btime fineshift!($x1, $y1, $ker, $off, 1, true)

    @printf("correl  1D")
    @btime correlate!($y2, $x2, $l, $w)

    @printf("correl' 1D")
    @btime correlate!($x2, $y2, $l, $w, 1, true)
end

function test2d(T::Type{<:AbstractFloat}=Float64)
    m1,m2 = 33,33
    n1,n2 = 40,40
    ref = rand(T,n1,n2)
    img = Array{T}(undef,m1,m2)
    ws1 = Array{T}(undef,m1,n2)
    ws2 = Array{T}(undef,n1,m2)
    off = map(T, (0.2,-0.3))
    ker = CatmullRomSpline(T)
    l1, w1 = FineShift.shiftparams(ker, off[1])
    l2, w2 = FineShift.shiftparams(ker, off[2])

    img1 = similar(img)
    img2 = similar(img)
    fineshift!(img1, fineshift!(ws1, ref, ker, off[1], 1, false),
               ker, off[2], 2, false)
    fineshift!(img2, fineshift!(ws2, ref, ker, off[2], 2, false),
               ker, off[1], 1, false)
    @printf("(a) max. abs. dif. direct: %g\n", maximum(abs.(img2 - img1)))
    correlate!(img2, correlate!(ws1, ref, l1, w1, 1, false), l2, w2, 2, false)
    @printf("(b) max. abs. dif. direct: %g\n", maximum(abs.(img2 - img1)))
    correlate!(img1, correlate!(ws2, ref, l2, w2, 2, false), l1, w1, 1, false)
    @printf("(c) max. abs. dif. direct: %g\n", maximum(abs.(img2 - img1)))

    @printf("direct  2D (1st then 2nd dim.)")
    @btime fineshift!($img, fineshift!($ws1, $ref, $ker, $(off[1]), 1, false),
                      $ker, $(off[2]), 2, false)

    @printf("direct  2D (2nd then 1st dim.)")
    @btime fineshift!($img, fineshift!($ws2, $ref, $ker, $(off[2]), 2, false),
                      $ker, $(off[1]), 1, false)

    @printf("adjoint 2D (1st then 2nd dim.)")
    @btime fineshift!($ref, fineshift!($ws2, $img, $ker, $(off[1]), 1, true),
                      $ker, $(off[2]), 2, false)

    @printf("adjoint 2D (2nd then 1st dim.)")
    @btime fineshift!($ref, fineshift!($ws1, $img, $ker, $(off[2]), 2, true),
                      $ker, $(off[1]), 1, false)

    @printf("correlate 2D (1st then 2nd dim.)")
    @btime correlate!($img, correlate!($ws1, $ref, $l1, $w1, 1, false),
                      $l2, $w2, 2, false)

    @printf("correlate 2D (2nd then 1st dim.)")
    @btime correlate!($img, correlate!($ws2, $ref, $l2, $w2, 2, false),
                      $l1, $w1, 1, false)

end

end

FineShiftBenchmarks.testalternatives()
FineShiftBenchmarks.test1d()
FineShiftBenchmarks.test2d()


text:
@testset "Axis arrays" begin
    @test ndims(zeros_axisarray(Int64, false)) == 0

    @test ndims(zeros_axisarray(Int64, true)) == 1
    @test axisnames(zeros_axisarray(Int64, true)) == (:spin,)
    @test axisvalues(zeros_axisarray(Int64, true)) == ([:↑, :↓],)

    unpol_time = zeros_axisarray(Int64, (:time, 0u"s":1u"s":9u"s"))
    @test ndims(unpol_time) == 1
    @test axisnames(unpol_time) == (:time,)
    @test size(unpol_time) == (10,)
    @test axisvalues(unpol_time) == (0u"s":1u"s":9u"s",)

    pol_time = zeros_axisarray(Int64, true, (:time, 0u"s":1u"s":9u"s"))
    @test ndims(pol_time) == 2
    @test size(pol_time) == (10, 2)
    @test axisnames(pol_time) == (:time, :spin)
    @test axisvalues(pol_time) == (0u"s":1u"s":9u"s", [:↑, :↓])

    pol_2d = zeros_axisarray(Int64, true, z=0u"s":1u"s":5u"s", x=0u"m":1u"m":10u"m")
    @test ndims(pol_2d) == 3
    @test size(pol_2d) == (6, 11, 2)
    @test axisnames(pol_2d) == (:z, :x, :spin)
    @test axisvalues(pol_2d) == (0u"s":1u"s":5u"s", 0u"m":1u"m":10u"m", [:↑, :↓])
end

@testset "Density arrays" begin
    ρ = density_array(Int64, true, x=0u"m":1u"m":10u"m")
    @test ndims(ρ) == 2
    @test size(ρ) == (11, 2)
    @test axisnames(ρ) == (:x, :spin)
    @test axisvalues(ρ) == (0u"m":1u"m":10u"m", [:↑, :↓])
    @test eltype(ρ) <: AtomicDFT.DFTUnits.Ρ
    @test eltype(ρ) <: AtomicDFT.HartreeUnits.ρ
    @test eltype(ρ) <: AtomicDFT.HartreeUnits.ρ{Int64}

    ρ = density_array(typeof(1u"nm^-3"), true, x=0u"m":1u"m":10u"m")
    @test eltype(ρ) <: AtomicDFT.DFTUnits.Ρ
    @test !(eltype(ρ) <: AtomicDFT.HartreeUnits.ρ)

    @test_throws ArgumentError density_array(typeof(1u"C*nm^-3"), true, x=0u"m":1u"m":10u"m")
end

@testset "is spin polarized?" begin
    @test is_spin_polarized(density_array(Int64, true, x=0u"m":1u"m":10u"m"))
    @test !is_spin_polarized(density_array(Int64, false, x=0u"m":1u"m":10u"m"))
end