text:
module Macros
using Graph, Recode, CodeMatch, Toivo
export @pattern

# ==== MethodTable ============================================================

type Method
    sig::Pattern
    f::Function
end

type MethodTable
    name::Symbol
    methods::Vector{Method}
    MethodTable(name::Symbol) = new(name, Method[])
end

function add(mt::MethodTable, m::Method)
    # insert the pattern in ascending topological order, as late as possible
    i = length(mt.methods)+1
    for (k, mk) in enumerate(mt.methods)
        if m.sig <= mk.sig
            if m.sig >= mk.sig
                # equal signature ==> replace
                mt.methods[k] = m
                return
            else
                i = k
                break
            end
        end
    end
    insert(mt.methods, i, m)

    for mk in mt.methods
        lb = m.sig & mk.sig
        if lb === nullpat; continue; end
        if any([ml.sig == lb for ml in mt.methods]) continue; end

        println("Warning: New @pattern method ", mt.name, m.sig)
        println("         is ambiguous with   ", mt.name, mk.sig)
        println("         Make sure ", mt.name, lb, " is defined first")
    end
end

function dispatch(mt::MethodTable, args::Tuple)
    for m in mt.methods
        matched, result = m.f(args)
        if matched; return result; end
    end
    error("No matching method found for pattern function $(mt.name)")
end

# function create_method(p::Pattern, body)
#     code = code_match(p)
#     f = @eval $argsym->begin
#         $(code)
#         (true, $body)
#     end
#     Method(p, f)
# end

function code_metod(p::Pattern, body)
    pred, bind = code_pat(p)
    fdef = :($argsym->begin
        if $pred
            $bind
            (true, $body)
        else
            (false, nothing)
        end
    end)
    fdef
end
create_method(p::Pattern, body) = Method(p, eval(code_metod(p,body)))

# ==== @pattern ===============================================================

const method_tables = Dict{Function, MethodTable}()

macro pattern(ex)
    code_pattern(ex)
end

function code_pattern(ex)
    sig, body = split_fdef(ex)
    @expect is_expr(sig, :call)
    fname, args = sig.args[1], sig.args[2:end]
    psig = :($(args...),)
    p_ex = recode(psig)

    f = esc(fname)
    quote
        wasbound = try
            f = $f
            true
        catch e
            false
        end

        if !wasbound
            mt = MethodTable($(quot(fname)))
            const $f = (args...)->dispatch(mt, args)
            method_tables[$f] = mt
        else
            if !has(method_tables, $f)
                error($("$fname is not a pattern function"))
            end
            mt = method_tables[$f]
        end

        method = create_method($p_ex, $(quot(body)))
        add(mt, method)
    end
end

end # module


text:
# Run all analyses and make all plots

using CSV, MixedModels, DataFrames, Statistics
include("zscore.jl")
include("read_data.jl")

# Set constants
surp_col_name = :SRN_1000_0_surprisal    # Name of the first data frame column that contains surprisal values
nrrep  = 6                               # Number of networks trained per network type
nrtest = 9                               # Number of test points over incremental training

# Read subject and surprisal data, and add previous-word surprisal values
print("Reading subject data and surprisal files\n")
subjects  = CSV.read("subjects.csv", delim='\t')
surprisal = CSV.read("surprisal.csv",delim='\t')

surp_col     = findall(names(surprisal).==surp_col_name)[1]            # Find column index of first network's surprisal value
prevsurp_col = size(surprisal,2)+1                                     # Columns for previous-word surprisal will start here

for c = surp_col:size(surprisal,2)
    col_name = String(names(surprisal)[c])
    surprisal[!,Symbol(col_name*"_prev")] = [missing; surprisal[1:end-1,c]]
end

include("do_lmm.jl")                      # Fit LMMs to explain human data from RNN surprisal estimates
include("compare_mean_surp.jl")           # Pairwise comparisons between network types on surprisal averaged over repetitions
include("do_gam.jl")                      # Fit GAMs to explain LMM-fit from RNN language model accuracy
include("make_plots.jl")                  # Make scatterplots and GAM curve plots for each of the human data sets

ext: jl

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct TensorShapeProto_Dim <: ProtoType
    size::Int64
    name::AbstractString
    TensorShapeProto_Dim(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorShapeProto_Dim

mutable struct TensorShapeProto <: ProtoType
    dim::Base.Vector{TensorShapeProto_Dim}
    unknown_rank::Bool
    TensorShapeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorShapeProto
const __fnum_TensorShapeProto = Int[2,3]
meta(t::Type{TensorShapeProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_TensorShapeProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export TensorShapeProto_Dim, TensorShapeProto


text:
# residual_evaluation.jl: some functions used by NonlinearSolvers relating to
# residual evaluation

import Utils.array1DTo3D

"""
  This function computes the vector form of of the residual from the vector
  form of the solution, ie. q_vec -> rhs_vec, for a given physics.
  This is one of the two functions required by [`newtonInner`](@ref).

  Users of newtonInner that are not physics modules (for example, implicit
  time marching schemes) will need to implement their own version of this
  function.  See the Implimentation Notes section below.

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an AbstractOperator
   * eqn: an AbstractSolutionData (may be modified during this function)
   * opts: the options dictionary
   * ctx_residual: a tuple of values.  ctx_residual[1] must be a function
                   that computes (q -> res).  Typically this is evalResidual.
                   The other entries of the tuple (if any) are not used.

                   The purpose of this argument is to make the signature of
                   the function generic enough so that implict time marching
                   methods can use it.  (If you are confused about this
                   programming pattern, google how callback are implemented
                   in C, this ctx is like a void* in C).

   * t: the time at which to evalute the residual

  **Inputs/Outputs**

   * rhs_vec: vector to put the residual in

  **Output**

   * norm of the residual vector

  **Implementation Notes:**

  This function is really a wrapper around an evalResidual-like function.
  It has to do 5 things:

    1. scatter eqn.q_vec -> eqn.q
    2. start parallel communication if needed
    3. call the evalResidual-like function to compute eqn.q -> eqn.res
    4. assemble the residual into the output vector, ie. eqn.res -> rhs_vec
    5. compute the norm of the rhs_vec

  Any implementation of this function (used with newtonInner) must have the
  following properties:

    1. on exit, eqn.q and eqn.q_vec must be consistent
    2. this function start parallel communication if needed
    3. allow for the possibility that rhs_vec and eqn.res_vec alias

  Is is recommended to use [`calcNorm`](@ref Utils.calcNorm) to compute the norm.

  Other implementations of this function are encouraged to use this function
  to help construct their rhs_vec.
"""
function physicsRhs(mesh, sbp, eqn, opts, rhs_vec, ctx_residual, t=0.0)

  # q_vec -> q
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec)

  # start parallel communication if needed
  time = eqn.params.time
  time.t_send += @elapsed if opts["parallel_type"] == 2 && mesh.npeers > 0
    startSolutionExchange(mesh, sbp, eqn, opts)
  end

  # evaluate (q -> res)
  func = ctx_residual[1]
  func(mesh, sbp, eqn, opts, t)

  # res -> rhs_vec
  assembleResidual(mesh, sbp, eqn, opts, rhs_vec, assemble_edgeres=false)

  # calculate the norm
  rhs_0_norm = calcNorm(eqn, rhs_vec, strongres=true)

  return rhs_0_norm
end

"""
### NonlinearSolvers.assembleResidual

  This function takes the residual in eqn.res and performs an additive reduction
  into res_vec.
  This function wraps array3DTo1D.

  Inputs:
    mesh: AbstractMesh
    sbp:  SBP operator
    eqn:  AbstractEquation object
    opts: options dictionary
    res_vec: residual vector to put the residual into

  Outputs:
    none

  Aliasing restrictions: none

"""
function assembleResidual(mesh, sbp, eqn, opts, res_vec::AbstractArray{T, 1};
                          assemble_edgeres=true, zero_resvec=true) where T
# assembles all of the residuals into res_vec
# no aliasing concerns
#TODO: clarify zero_resvec (why not used for both calls?)

  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, res_vec)

  if assemble_edgeres

    for i=1:size(eqn.res_edge, 4)
      array3DTo1D(mesh, sbp, eqn, opts, sview(eqn.res_edge, :, :, :, i),
                           res_vec, zero_resvec=zero_resvec)
    end
  end

  return nothing
end

@doc """
### NonlinearSolvers.array1DTo3D

  This function performs the scatter q_vec -> eqn.q

  Inputs:
    mesh: AbstractMesh
    sbp:  SBP operator
    eqn:  AbstractEquation object
    opts: options dictionary
    q_vec: vector containing solution variables

  Outputs:
    none

  Aliasing Restrictions: none

"""->
function array1DTo3D(mesh, sbp, eqn, opts, q_vec::AbstractArray{T, 1}) where T
# scatters the q_vec to the 3d array eqn.q
# no aliasing concerns here
  array1DTo3D(mesh, sbp, eqn, opts, q_vec, eqn.q)

  return nothing
end


text:
module ParseTrees

import Base.Iterators: Generator, Filter
import LightXML: parse_file, root, name, child_elements, find_element, content, attribute
import LightGraphs: DiGraph, add_edge!, indegree, has_path, outneighbors, rem_edge!, vertices
import MetaGraphs: MetaDiGraph, set_props!, get_prop

idx(clause) = parse(Int, attribute(clause, "idx"))

export sentences
"""
    sentences(file)

Return an iterator over sentences of an `coreNLP` rules XML file.

You can run `coreNLP` through the [`command-line`](https://stanfordnlp.github.io/CoreNLP/cmdline.html).
This package only uses the [`parse`](https://stanfordnlp.github.io/CoreNLP/parse.html)
annotator.

```jldoctest
julia> using ParseTrees

julia> sentence = sentences("hammurabi.txt.xml")[2]
{29, 28} directed Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)

julia> flat(sentence)
"If any one ensnare another , putting a ban upon him , but he can not prove it , then he that ensnared him shall be put to death"

julia> flat(sentence, 4)
"If any one ensnare another , putting a ban upon him"
```
"""
sentences(file) = map(
    xml_sentence -> begin
        words = collect(child_elements(first(Filter(
            annotation ->
                name(annotation) == "dependencies" &&
                attribute(annotation, "type") == "basic-dependencies",
            child_elements(xml_sentence)
        ))))
        sentence = MetaDiGraph(DiGraph(length(words) - 1))
        foreach(
            word -> begin
                governor = idx(find_element(word, "governor"))
                dependent_word = find_element(word, "dependent")
                dependent = idx(dependent_word)
                if governor != 0
                    add_edge!(sentence, governor, dependent)
                end
                set_props!(sentence, dependent, Dict(
                    :type => attribute(word, "type"),
                    :content => content(dependent_word)
                ))
            end,
            words
        )
        sentence
    end,
        child_elements(find_element(
            find_element(root(parse_file(file)), "document"),
            "sentences"
        ))
)

flat_inner(sentence, words) = join(Generator(
    word -> get_prop(sentence, word, :content),
    words
), " ")

"""
    flat(sentence)

The whole sentence

    flat(sentence, root)

Just the parts of the `sentence` connected to `root`.
"""
flat(sentence) = flat_inner(sentence, vertices(sentence))
flat(sentence, root) = flat_inner(sentence, Filter(
    word -> has_path(sentence, root, word),
    vertices(sentence)
))
export flat

"""
    const deontic

Default pattern to recognize `deontic`s.
"""
const deontic = r"(should)|(will)|(shall)|(must)|(may)|(can)"
export deontic

"""
    dictionary

A dictionary based on Elinor Ostrom's [`institutional grammar`](https://www.jstor.org/stable/2082975).
"""
dictionary = Dict(
    "acl" => :not_in_root, # adjectival clause
    "advcl" => :condition,  # adverbial clause modifier
    "advmod" => :condition, # adverbial modifier
    "amod" => :not_in_root, # adjectival modifier
    "apos" => :not_in_root, # appositional modifier
    "aux" => :deontic, # auxiliary
    "auxpass" => :root, # passive auxiliary
    "case" => :not_in_root, # case marking
    "cc" => :remove, # coordinating conjunction
    "ccomp" => :object, # clausal component
    "compound" => :root,
    "conj" => :recur, # conjunct
    "cop" => :root, # copula
    "csubj" => :attribute, # clausal subject
    "csubjpass" => :attribute, # clausal subject passive
    "dep" => :not_sure, # unspecified dependency
    "det" => :not_in_root, # determiner
    "discourse" => :remove,
    "dislocated" => :not_sure,
    "dobj" => :object, # direct object
    "expl" => :root, # expletive
    "foreign" => :not_sure, # foreign words
    "goeswith" => :root, # goes with
    "iobj" => :root, # indirect object
    "list" => :recur,
    "mark" => :not_in_root, # marker
    "mwe" => :root, # multiword expression
    "name" => :not_in_root,
    "neg" => :root, # negation modifier
    "nmod" => :condition, # nominal modifier
    "nsubj" => :attribute, # nominal subject
    "nsubjpass" => :attribute, # nominal subject passive
    "nummod" => :not_in_root, # numeric modifier
    "parataxis" => :recur,
    "punct" => :remove, # punctuation
    "remnant" => :not_sure, # remnant in ellipisis
    "reparandum" => :not_sure, # overridden disfluency
    "root" => :root,
    "vocative" => :remove,
    "xcomp" => :object, # open clausal component
)
export dictionary

function clauses!(rules, sentence, dictionary, deontic, root)
    clauses = Vector{Pair{Int, Symbol}}()
    is_rule = false
    for clause in copy(outneighbors(sentence, root))
        component = dictionary[split(get_prop(sentence, clause, :type), ':')[1]]
        if component == :recur
            clauses!(rules, sentence, dictionary, deontic, clause)
            rem_edge!(sentence, root, clause)
        elseif component == :deontic
            if occursin(deontic, flat(sentence, clause))
                is_rule = true
                rem_edge!(sentence, root, clause)
                push!(clauses, clause => component)
            else
                component = :root
            end
        elseif component == :remove
            rem_edge!(sentence, root, clause)
        elseif component != :root
            rem_edge!(sentence, root, clause)
            push!(clauses, clause => component)
        end
    end
    if is_rule
        push!(clauses, root => :root)
        push!(rules, map(
            clause -> clause.second => flat(sentence, clause.first),
            clauses
        ))
    end
end

"""
    split_rules(file; dictionary = dictionary, deontic = deontic)

Split `file` into rules. Split each rule into components based on `dictionary`,
starting at the root. The `dictionary` should be a map from [Universal
Dependencies (v1)](http://universaldependencies.org/docsv1/u/dep/all.html) to
clause components. `:deontic`s must additionally match the `deontic` pattern.
Rules must contain a deontic. `:recur` will look for a rule inside the rule,
`:remove` will remove the clause, and `:root` will include the root.

```jldoctest
julia> using ParseTrees

julia> parsed = split_rules("hammurabi.txt.xml");

julia> parsed[25]
6-element Array{Pair{Symbol,String},1}:
 :condition => "If the slave will not give the name of the master"
 :attribute => "the finder"
   :deontic => "shall"
    :object => "him"
 :condition => "to the palace"
      :root => "bring"
```
"""
function split_rules(file; dictionary = dictionary, deontic = deontic)
    rules = Vector{Vector{Pair{Symbol, String}}}()
    for sentence in sentences(file)
        roots = filter(word -> indegree(sentence, word) == 0, vertices(sentence))
        if length(roots) >= 1
            clauses!(rules, sentence, dictionary, deontic, first(roots))
        end
    end
    rules
end
export split_rules

end


text:
# Depth-first visit / traversal

#################################################
#
#  Depth-first visit
#
#################################################

mutable struct DepthFirst <: AbstractGraphVisitAlgorithm
end

function depth_first_visit_impl!(
    graph::AbstractGraph{V,E},      # the graph
    stack,                          # an (initialized) stack of vertex
    vertexcolormap::Vector{Int},    # an (initialized) color-map to indicate status of vertices
    edgecolormap::Vector{Int},      # an (initialized) color-map to indicate status of edges
    visitor::AbstractGraphVisitor) where {V,E} # the visitor

    while !isempty(stack)
        u, uegs, tstate = pop!(stack)
        found_new_vertex = false
        while tstate <= length(uegs) && !found_new_vertex #!done(uegs, tstate)
            v_edge, tstate = iterate(uegs, tstate)
            # v_edge, tstate = next(uegs, tstate)
            v = v_edge.target
            v_color = vertexcolormap[vertex_index(v, graph)]
            e_color = edgecolormap[edge_index(v_edge, graph)]
            examine_neighbor!(visitor, u, v, v_color, e_color)

            if e_color == 0
                edgecolormap[edge_index(v_edge, graph)] = 1
            end

            if v_color == 0
                found_new_vertex = true
                vertexcolormap[vertex_index(v, graph)] = 1
                if !discover_vertex!(visitor, v)
                    return
                end
                push!(stack, (u, uegs, tstate))

                open_vertex!(visitor, v)
                vegs = out_edges(v, graph)
                push!(stack, (v, vegs, 1)) # push!(stack, (v, vegs, start(vegs)))
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vertexcolormap[vertex_index(u, graph)] = 2
        end
    end
end

function traverse_graph(
    graph::G,
    alg::DepthFirst,
    s::V,
    visitor::AbstractGraphVisitor;
    vertexcolormap = zeros(Int, num_vertices(graph)),
    edgecolormap = zeros(Int, num_edges(graph))) where {V,G <: AbstractGraph}

    @graph_requires graph incidence_list vertex_map

    vertexcolormap[vertex_index(s, graph)] = 1
    if !discover_vertex!(visitor, s)
        return
    end

    segs = out_edges(s, graph)
    sstate = 1 #sstate = start(segs)
    stack = [(s, segs, sstate)]

    depth_first_visit_impl!(graph, stack, vertexcolormap, edgecolormap, visitor)
end

#################################################
#
#  Useful applications
#
#################################################

# Test whether a graph is cyclic

mutable struct DFSCyclicTestVisitor <: AbstractGraphVisitor
    found_cycle::Bool

    DFSCyclicTestVisitor() = new(false)
end

function examine_neighbor!(
    vis::DFSCyclicTestVisitor,
    u::V,
    v::V,
    vcolor::Int,
    ecolor::Int) where {V}

    if vcolor == 1 && ecolor == 0
        vis.found_cycle = true
    end
end

discover_vertex!(vis::DFSCyclicTestVisitor, v) = !vis.found_cycle

function test_cyclic_by_dfs(graph::G) where {G <: AbstractGraph}
    @graph_requires graph vertex_list incidence_list vertex_map

    cmap = zeros(Int, num_vertices(graph))
    visitor = DFSCyclicTestVisitor()

    for s in vertices(graph)
        if cmap[vertex_index(s, graph)] == 0
            traverse_graph(graph, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end

        if visitor.found_cycle
            return true
        end
    end
    return false
end

# Topological sort using DFS

mutable struct TopologicalSortVisitor{V} <: AbstractGraphVisitor
    vertices::Vector{V}
    function TopologicalSortVisitor{V}(n::Int) where V
      vs = Array{Int}(undef, 0)
      sizehint!(vs, n)
      new{V}(vs)
    end
end

function examine_neighbor!(visitor::TopologicalSortVisitor{V}, u::V, v::V, vcolor::Int, ecolor::Int) where {V}
    if vcolor == 1 && ecolor == 0
        throw(ArgumentError("The input graph contains at least one loop."))
    end
end

function close_vertex!(visitor::TopologicalSortVisitor{V}, v::V) where {V}
    push!(visitor.vertices, v)
end

function topological_sort_by_dfs(graph::AbstractGraph{V}) where {V}
    @graph_requires graph vertex_list incidence_list vertex_map

    cmap = zeros(Int, num_vertices(graph))
    visitor = TopologicalSortVisitor{V}(num_vertices(graph))
    for s in vertices(graph)
        if cmap[vertex_index(s, graph)] == 0
            traverse_graph(graph, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end
    end

    reverse(visitor.vertices)
end


text:
# """ REWRITE:
# initback initializes the fields used by Net.back:
# - dif0, tmp: allocated or size checked.
# - reg.incr: depends on seq.
# - reg.back: depends on which dx args specified.
# - op.save: read-only, used for popping only if seq.
# """
function initback(f::Net, ygold, loss, getdx, seq)
    seq == !stack_isempty(f) || error("")
    nextback = (f.lastforw, getdx, seq)
    if isvoid(f,:lastback) || f.lastback != nextback
        initgrad(f, getdx)
        initincr(f)
        inittype2(f)
        initdif0(f)
        inittmp(f)
        f.lastback = nextback
        gpusync()
    end
end

backregs(f::Net)=filter(r->getp(r,:grad), regs(f))

initdif0(f::Net)=(for p in reverse(backregs(f)); initdif0(f,p); end)

function initdif0(f::Net,p::Reg)
    at, et, sz = getp(p,:diftype), getp(p,:eltype), getp(p,:size)
    (isa(at,DataType) && isa(et,DataType) && isa(sz,Dims)) || error("Reg has missing info: $((p.name,at,et,sz))")
    if !checkarray(p, :dif0, at, et, sz)
        p.dif0 = newarray(at, et, sz)
        getp(p,:incr) && fillsync!(p.dif0,0)
    end
    if canoverwrite(p.op) && getp(p,:backoverwrite,true) && !getp(p,:incr)
        for i in reverse(p.argv) # consider overwriting the last input first
            q = reg(f,i)
            if (isvoid(q,:dif0) && !getp(q,:incr) && (at,et,sz)==(getp(q,:diftype),getp(q,:eltype),getp(q,:size)))
                q.dif0 = p.dif0
                break
            end
        end
    end
end

function inittmp(f::Net)
    for p in backregs(f)
        if getp(p,:incr)
            at, et, sz = getp(p,:tmptype), getp(p,:eltype), getp(p,:size)
            checkarray(p,:tmp,at,et,sz) || (p.tmp = newarray(at,et,sz))
        end
    end
end

# """ REWRITE:
# set_toback(f::Net) sets f.toback[n] which is true if dif[n] should be
# calculated for op[n] during back calculation.  This is only needed if
# op[n] is a par node or a par node descendent.  Or if the caller asked
# for dx for network inputs, those and their descendents.
# """
function initgrad(f::Net, getdx)
    for p in regs(f)
        setp(p,:grad,false)
    end
    lastinput = 0
    for p in regs(f)
        if ((isa(p.op, Par) && p.op.initialized) || (isa(p.op, Input) && getdx[lastinput+=1]))
            setp(p,:grad,true)  # :grad means we need the gradient of the output of this operation, does not mean go back on this operation.
        end
    end
    count_grad(f::Net)=mapreduce(r->getp(r,:grad), +, 0, regs(f))
    nback = count_grad(f)
    while true
        for p in regs(f)
            getp(p,:grad) && continue
            for i in inputregs(f,p)
                if getp(i,:grad)
                    setp(p,:grad,true)
                    break
                end
            end
        end
        nb = count_grad(f)
        nb == nback ? break : nback = nb
    end
end

# REWRITE:
# set_toincr(f::Net) sets f.toincr[n] which is true if dif[n] should be
# incrementally updated.  This is necessary if op[n] has multiple
# outputs, or it is a Par and we are processing a sequence.  This is the
# one place we can't seem to get rid of the seq flag.
#
# We can make all Par registers :incr.
# - TODO: efficiency cost for fnn.
# - who resets them to zero for rnn or for fnn.

function initincr(f::Net)
    seq = (stack_length(f) > 0)
    for p in backregs(f)
        setp(p,:incr,seq && isa(p.op,Par))
        setp(p,:fanout, isreturn(p) ? 1 : 0)
    end
    for p in backregs(f)
        for i in inputregs(f,p)
            if getp(i,:grad) && incp(i,:fanout) > 1
                setp(i,:incr,true)
            end
        end
    end
end

# Note about sparse arrays:
# If input x is sparse (example: matrix of one-hot word columns)
# and the first op is dot(w,x):
# w will be dense (example: word embedding matrix)
# y = w * x and all following forw outputs will be dense.
# Going back we have dw = dy * x' and dx = w' * dy
# We would like to keep dw sparse.
# The only sparse matrices are the dw and iw for w that dot sparse inputs.
# Regular sparse x is a right csc input.
# NCE adds noise and ygold matrices as left csr inputs.
# TODO: cleanup this mess: more generic code? more versatile sparse type?

function inittype2(f::Net)
    for p in backregs(f)
        setp(p,:diftype, gpu() ? CudaArray : Array)
        setp(p,:tmptype, gpu() ? CudaArray : Array)
    end
    for p in backregs(f)
        isa(p.op, Dot) || continue # We only know how to do sparse dot
        (a, b) = inputregs(f, p)
        if issparse(a.out0) && issparse(b.out0)
            error("Dot of two sparse matrices")
        elseif issparse(b.out0) && getp(a,:grad) # y = w * x  with sparse x in rnnlm
            if getp(a,:incr)
                setp(a,:diftype, gpu() ? CudaArray : Array)
                setp(a,:tmptype, gpu() ? CudaSparseMatrixCSRU : SparseMatrixCSC)
            else
                setp(a,:diftype, gpu() ? CudaSparseMatrixCSR : SparseMatrixCSC)
            end
        elseif issparse(a.out0) && getp(b,:grad) # rw = r * w  with sparse r in nce
            if getp(b,:incr)
                setp(b,:diftype, gpu() ? CudaArray : Array)
                setp(b,:tmptype, gpu() ? CudaSparseMatrixCSCU : SparseMatrixCSC)
            else
                setp(b,:diftype, gpu() ? CudaSparseMatrixCSC : SparseMatrixCSC)
            end
        end
    end
end

### DEAD CODE

# function inittype2(f::Net)
#     fillsync!(f.sparse, nothing)
#     for o=1:length(f.op)
#         if isa(f.op[o], Dot)
#             (i,j) = f.inputs[o]
#             if issparse(f.out0[i]) && issparse(f.out0[j])
#                 error("Dot of two sparse matrices")
#             elseif issparse(f.out0[j])
#                 f.sparse[i] = f.toincr[i] ? CudaSparseMatrixCSRU : CudaSparseMatrixCSR # todo: clean up this mess!
#             elseif issparse(f.out0[i])
#                 f.sparse[j] = f.toincr[i] ? CudaSparseMatrixCSCU : CudaSparseMatrixCSC
#             end
#         end
#     end
# end

# function finddif(f::Net, n)
#     dif0 = nothing
#     if !f.toincr[n]
#         if length(f.outputs[n]) == 0
#             # There are no outputs to try overwriting, this may happen with childless ops
#         elseif length(f.outputs[n]) == 1
#             o = f.outputs[n][1]     # first try overwriting the output dif
#             if (o > n
#                 && !f.toincr[o]
#                 && overwrites(f.op[o])
#                 && (f.dif0[o]!=nothing)
#                 && size(f.dif0[o]) == size(f.out0[n]))
#                 dif0 = f.dif0[o]
#             end
#         else
#             error("toincr==false for multi output op")
#         end
#     end
#     if dif0 == nothing
#         et = eltype(f.out0[n])
#         sz = size(f.out0[n])
#         dif0 = (
#                 ### This is sparse non-incremental dw: can't really have dense without
#                 # rewriting CUSPARSE.csrmm to take sparse matrix in second position.  As
#                 # it stands, we'd have to transpose all three matrices: dw = dy * x' -> dw' = x * dy'
#                 (f.sparse[n]!=nothing) && !f.toincr[n] && !gpu() ? spzeros(et, sz...) :
#                 (f.sparse[n]!=nothing) && !f.toincr[n] && gpu()  ? (f.sparse[n])(spzeros(et, sz...)) : # t:12.38
#                 ### CSRU speed 20% slower than CSR on mnist (atomicAdd conflicts?), also cannot compute vecnorm.
#                 # f.sparse[n] && !f.toincr[n] && gpu()  ? CudaSparseMatrixCSRU(et, sz...) : # t:14.69
#                 ### Uncomment this if you want sparse incremental dw:
#                 ### Speed similar to dense on rnnlm, less memory, cannot compute vecnorm.
#                 # f.sparse[n] && f.toincr[n] && gpu()   ? ArrayAccumulator(et, sz) :
#                 # f.sparse[n] && f.toincr[n] && !gpu()  ? ArrayAccumulator(et, sz) :
#                 gpu() ? fillsync!(CudaArray(et, sz),0) : zeros(et, sz))
#     end
#     return dif0
# end

# function findtmp(f::Net, n)
#     tmp = nothing
#     for i=n+1:length(f.op)
#         if ((f.tmp[i]!=nothing) &&
#             (size(f.tmp[i]) == size(f.dif0[n])) &&
#             (f.sparse[n]==nothing || typeof(f.tmp[i]) <: f.sparse[n]))
#             tmp = f.tmp[i]
#             break
#         end
#     end
#     if tmp == nothing
#         et = eltype(f.out0[n])
#         sz = size(f.out0[n])
#         tmp = (
#                # CSRU is 5% faster if no atomic op conflicts (rnnlm),
#                # but significantly slower when there are lots of conflicts (mnist)
#                # Not worth the risk until I implement uniq for CSRU
#                # Seems significantly faster on s2s, putting csru back on
#                gpu() && f.sparse[n]!=nothing ? (f.sparse[n])(et, sz...) :
#                # gpu() && f.sparse[n] ? CudaSparseMatrixCSR(spzeros(et, sz...)) :
#                !gpu() && f.sparse[n]!=nothing ? spzeros(et, sz...) :
#                gpu() ? CudaArray(et, sz) :
#                Array(et, sz))
#     end
#     return tmp
# end

# if isa(f.op[i], Input) && issparse(f.out0[i])
#     for o in f.outputs[i]
#         if isa(f.op[o], Dot)
#             w = f.inputs[o][1]
#             @assert isa(f.op[w], Par)
#             f.sparse[w] = true
#         end
#     end
# end

# # instead of changing the ops, record the eltype in net and use it when generating stuff.
#     # for o in f.op
#     #     t = eltype(o)
#     #     if t == nothing
#     #         continue
#     #     elseif it == nothing
#     #         it = t
#     #     elseif t != it
#     #         seteltype!(o, it)
#     #     end
#     # end

# #####

# function netinit(f::Net)
#     initout(f, inputs...)
#     initdif(f, inputs...)
#     inittmp(f, inputs...)
# end

# """

# initout(f::Net, inputs...) initializes f.out[n] which points to the
# array holding the output of op[n].  f.out0[n] holds the actual array,
# allowing f.out[n] point to stack elements when necessary during back
# calculation.  We optimize memory use by sharing arrays when we can.

# """
# function initout(f::Net)
#     N = length(f.op)
#     index = zeros(Int,N)          			# index==0 represents need for new register
#     for n=1:N
#         overwrites(f.op[n]) || continue                 # see if we can reuse our input for output
#         f.tosave[n] && continue                         # a saved register should only be written by op[n]
#         i = f.inputs[n][1]                              # see if we can overwrite the first input
#         f.tosave[i] && continue                         # do not overwrite if you are going to save for back
#         ow = true; k = n                                # todo: checkout other inputs if this one fails
#         while true                                      # see if anybody else uses i before its next update
#             k = mod1(k+1, N)                            # todo: This is suboptimal, gives 13 regs for LSTM
#             in(i, f.inputs[k]) && (ow = false; break)   # Should look for existing regs no longer used
#             k == i && break                             # this is nontrivial because the sizes are unknown
#             k == N && i > N && break                    # but can possibly be inferred
#         end                                             # fencepost check: N==1:OK, i==n:OK, i>N:OK
#         ow && (index[n]=i)
#     end
#     f.out0 = cell(N)
#     for n=1:N                                           # index==0 ops get their own registers
#         index[n]==0 || continue
#         f.out0[n] = Any[]
#     end
#     for n=1:N                                             # other ops will overwrite existing registers
#         index[n]==0 && continue
#         k = index[n]
#         while index[k]!=0; k=index[k]; end
#         f.out0[n] = f.out0[k]
#     end
#     f.out = fillsync!(cell(N), nothing)                       # use nothing to represent identity element: 0 for add, 1 for mul
# end

# """
# initdif(f::Net) initializes f.dif[n], which is the loss gradient of the last output of op[n].
# f.dif0[n] holds the actual array, f.dif[n] can point to this or be 'nothing' representing the zero matrix.
# """
# function initdif(f::Net)
#     # todo: this is not optimal for add: no need for two copies of dy when neither input is overwriting.
#     # todo: tmp register to replace dif1 and serve other needs

#     index = zeros(Int,nops(f))                                  # index==0 represents need for new register
#     for n=1:nops(f)                                             # find out if back(op[n]) can overwrite dif[n]
#         overwrites(f.op[n]) || continue                         # overwrites means potentially both forw x<-y and back dy<-dx
#         f.toincr[n] && continue                                  # don't share dif[n] with multi-output
#         n == nops(f) && continue                                # don't share last output, dy may come in sparse
#         for i in f.inputs[n]
#             f.toincr[i] && continue                              # don't share dif[n] with multi-output
#             index[i]==0 || error("should not happen")
#             index[i]=n; break                                   # op[n] will overwrite dif[n] to get dif[i]
#         end
#     end
#     f.dif0 = cell(nops(f))
#     for n=1:nops(f)
#         index[n]==0 || continue
#         f.dif0[n] = Any[]
#     end
#     for n=1:nops(f)
#         index[n]==0 && continue
#         k = index[n]
#         while index[k]!=0; k=index[k]; end
#         f.dif0[n] = f.dif0[k]
#     end
#     f.dif = fillsync!(cell(nops(f)), nothing)
# end

# """
# inittmp(f::Net) initializes f.tmp[n] which is an extra array for:
# (1) incremental updates if f.toincr[n].
# (2) for forw/back ops that need extra space (todo)
# """
# function inittmp(f::Net)
#     # todo: this is inefficient, could use a single tmp for each size.
#     f.tmp = cell(nops(f))
#     for n=1:length(f.tmp)
#         f.tmp[n] = f.toincr[n] ? Any[] : nothing
#     end
# end

# ### DEAD CODE:

#     # initop(f, a...)
#     # initinputs(f, a...)
#     # @assert length(f.op)==length(f.inputs)
#     # initninputs(f)
#     # initparams(f)
# # function initstack(f::Net)
# #     f.stack = Any[]
# #     f.sp = 0
# # end

# #    initstack(f)

# # f.op[n] is the n'th operation in the net
# # The user specifies the operations in the MLP constructor arguments
# # If an argument is another MLP, its operations are spliced in

# # function initop(f::Net, a...)
# #     f.op = Op[]
# #     for ai in a
# #         isa(ai,Tuple) && (ai=ai[1])
# #         isa(ai,Op) ?  push!(f.op, ai) :
# #         isa(ai,Net)   ?  append!(f.op, ai.op) :
# #         error("Bad op: $ai")
# #     end
# # end

# # f.inputs[n] is an array of k indices for the inputs of the n'th op.
# # k is typically 1 but can be 2 for Add2 and Mul2 operations
# # index i<=nops(f) indicates output of f.op[i], i>nops(f) indicates network inputs
# # By default an op takes the results of the previous k op outputs
# # (or network inputs for initial ops)
# # The user can override this by providing a tuple argument for an op
# # The first element of the tuple is the op, the rest are user indices for inputs
# # userindex j>0 indicates output of userarg[j], j<=0 indicates network input 1-j.

# # function initinputs(f::Net, a...)
# #     newindex = Array(Int, length(a))
# #     lastindex = 0
# #     for i=1:length(a)
# #         ai = isa(a[i],Tuple) ? a[i][1] : a[i]
# #         lastindex += (isa(ai,Op) ? 1 : length(ai.op))
# #         newindex[i] = lastindex
# #     end
# #     f.inputs = Any[]
# #     for i=1:length(a)
# #         ai = isa(a[i],Tuple) ? a[i][1] : a[i]
# #         bi = isa(a[i],Tuple) ? a[i][2:end] : ((i-ninputs(ai)):(i-1))
# #         length(bi) == ninputs(ai) || error("Wrong number of inputs for $i:$(typeof(ai))")
# #         if isa(ai, Op)
# #             push!(f.inputs, map(j->(j>0 ? newindex[j] : nops(f)+1-j), Int[bi...]))
# #         else
# #             j0 = length(f.inputs)
# #             for aii in ai.inputs
# #                 push!(f.inputs, map(j->(j<=nops(ai) ? j+j0 : (j=bi[j-nops(ai)]; j>0 ? newindex[j] : nops(f)+1-j)), aii))
# #             end
# #         end
# #     end
# # end

# # f.ninputs is the number of inputs the whole MLP expects
# # indices i>nops(f) in f.inputs refer to network inputs

# # function initninputs(f::Net)
# #     n = 0
# #     for ai in f.inputs
# #         for aij in ai
# #             aij - nops(f) > n && (n = aij - nops(f))
# #         end
# #     end
# #     f.ninputs = n
# # end

# # f.params points to all op parameters

# # function initparams(f::Net)
# #     f.params = Any[]
# #     for o in f.op
# #         if isa(o,par) && isupdated(o)
# #             push!(f.params, o)
# #         end
# #     end
# # end

# # isupdated(o)=true # todo: should be false for rand and const

#     #     @assert !haskey(dict, target)			# each target has to be unique
#     #     for x in args
#     #         @assert haskey(dict, x)                     #
#     #     end
#     #     if func == :input && inputs != nothing          # subroutines take their input registers from caller
#     #         @assert length(inputs) >= nextinput
#     #         dict[target] = inputs[nextinput]
#     #         nextinput += 1
#     #     end
#     #     op = eval(Expr(:call, func, pars...))
#     # end
#     #

# # initforw(f::Net,x...) is called at the beginning of a sequence or
# # before processing a stand-alone item.  It is not called between
# # elements of a sequence.  It allocates and/or resizes f.out0.

# # function initforw2(f::Net, inputs...; keepstate=false, a...)
# #     f.dbg && display((:initforw0,keepstate,vecnorm0(f.out),vecnorm0(f.stack[1:f.sp])))
# #     f.sp == 0 || error("Stack corruption")
# #     length(inputs) == ninputs(f) || error("Wrong number of inputs")
# #     out = fillsync!(cell(length(f.out0)), nothing)          # todo: (minor) can we get rid of alloc
# #     lastinput = 0
# #     while findfirst(out,nothing) > 0                            # allow multiple passes for size inference
# #         nnothing = count(x->(x==nothing), out)
# #         for n = 1:nops(f)
# #             if out[n] != nothing
# #                 continue
# #             elseif isa(f.op[n], Input)
# #                 out[n] = initarray(f.out0, n, inputs[lastinput += 1])
# #             elseif isempty(f.inputs[n])                         # par, rnd, con
# #                 s = psize(f, n)                                 # todo: write psize inference
# #                 s == nothing && continue                        # may happen first pass
# #                 out[n] = initarray(f.out0, n, inputs[1], s; dense=true)
# #             else
# #                 s = ysize(f.op[n], out[f.inputs[n]]...)
# #                 s == nothing && continue # may happen with recurrent connections
# #                 out[n] = initarray(f.out0, n, out[f.inputs[n]][1], s; dense=true)
# #             end
# #         end
# #         nnothing == count(x->(x==nothing), out) && error("Cannot determine size of array $(findfirst(out,nothing))")
# #     end
# #     # We recover or reset f.out:
# #     keepstate ? copysync!(f.out, f.out0) : fillsync!(f.out, nothing) # todo: f.out[i] = arr(f.out0[i]) elsewhere..
# #     f.dbg && display((:initforw1,keepstate,vecnorm0(f.out),vecnorm0(f.stack[1:f.sp])))
# # end

# # function initarray(a, i, x, dims=size(x); dense=false)
# #     (display((:initarray0,summary(a[i]),i,summary(x),dims,dense,summary(a)));println())
# #     if isempty(a[i])
# #         oldai = a[i]
# #         if !dense && issparse(x)
# #             a[i] = spzeros(eltype(x), dims...)
# #             gpu() && (a[i] = CudaSparseMatrixCSC(a[i]))
# #         else
# #             a[i] = fillsync!(KUdense(gpu()?CudaArray:Array, eltype(x), dims), 0)
# #         end
# #         for j=1:length(a); a[j]===oldai && (a[j]=a[i]); end # preserve array sharing
# #     elseif eltype(a[i]) != eltype(x)
# #         error("Element type mismatch")
# #     elseif size(a[i]) != dims
# #         warn("Resizing $(size(a[i]))->$dims")
# #         fillsync!(resize!(a[i], dims), 0) # todo: this does not work for sparse
# #     end
# #     (display((:initarray1,summary(a[i]),i,summary(x),dims,dense,summary(a)));println())
# #     return a[i]
# # end

#             # elseif in(nothing, dims[f.inputs[n]])
#             #     continue
#             # else
#             # end

#             # if isempty(f.inputs[n])                     # happens with par, rnd, con, input
#             #     elseif prod(size(f.op[n])) > 0
#             #         dims[n] = size(f.op[n])
#             #     end
#             # end
#             # if isa(op, Input)
#             #     dims[n] = size(inputs[lastinput += 1])
#             # elseif isa(op, Par) && prod(op.dims) > 0
#             #     dims[n] = op.dims
#             # elseif isa(op, Dot)
#             #     (i1,i2) = f.inputs[n]
#             #     (d1,d2) = dims[i1,i2]
#             #     if d1 == nothing || d2 == nothing
#             #         continue
#             #     end
#             #     m1 = prod(d1[1:end-1])                   # treat [d...,dn] as [prod(d...), dn] matrix
#             #     n1 = d1[end]
#             #     m2 = prod(d2[1:end-1])
#             #     n2 = d2[end]
#             #     if m1 == 0 || n2 == 0
#             #         continue
#             #     end
#             #     if n1 == m2 == 0
#             #         continue
#             #     elseif n1 == 0
#             #         setdims!(f.op[i1], (d1[1:end-1]..., m2))
#             #     elseif m2 == 0
#             #         length(d2) == 2 || error("Cannot infer size of tensor")
#             #         setdims!(f.op[i2], (n1, d2[2]))
#             #     elseif n1 == m2
#             #         dims[n] = (m1,n2)
#             #     else
#             #         error("Dot size mismatch")
#             #     end
#             # end
#             # # input, par, dot, add, mul, actf, loss, conv, pool

# function initback0(f::Net, dy; a...)
#     N = length(f.op)
#     firstwrite = zeros(Int, N)                          # dif[n] is read at t=n, find out when it is first written.
#     # todo: dif[N] is special!
#     for n=N:-1:1
#         f.toback[n] || continue
#         k = n
#         while true
#             in(n, f.inputs[k]) && (firstwrite[n] = k; break)
#             k = mod1(k-1, N)
#         end
#     end
#     return firstwrite
#     for n=N:-1:1
#         nsparse = false # todo
#         f.toback[n] || continue
#         f.dif[n] = finddif(f, n, nsparse)
#         f.toincr[n] && (f.tmp[n] = findtmp(f, n, nsparse))
#     end
# end

# function finddif(f::Net, n, nsparse)
#     isa(f.op[n], Par) && return                         # par difs should persist across iterations
#     N = length(f.op)
#     k = n; nw1 = 0                                      # dif[n] is read at t=n, find out when it is first written.
#     while true
#         in(n, f.inputs[k]) && (nw1 = k; break)          # it could potentially be written at t=n right after reading
#         k = mod1(k-1, N)
#     end
#     @assert nw1 > 0                                     # dif[n] will be busy from nw1 to n going back cyclical
#     free = nothing                                      # we need dif[i] that is free during this period
#     for i = n+1:N                                       # note that dif[i] may already be shared to some dif[j]
#         isa(f.op[i], Par) && return
#         size(f.dif[i]) == size(f.out0[n]) || continue
#         issparse(f.dif[i]) == nsparse || continue

#     end
#     return free
# end

# initback(f::Net, dy) called at the beginning of a sequence or a
# stand-alone item, never between elements of a sequence.

# function initback1(f::Net, dy; seq=false, a...)
#     fillsync!(f.dif, nothing)                           # why? (todo)
#     for n=1:length(f.dif0)
#         y = (n==nops(f) && dy!=nothing ? dy : f.out0[n])
#         initarray(f.dif0, n, y; dense=true) # x and dw may be sparse, dx and w always dense
#         if f.toincr[n]
#             initarray(f.tmp, n, f.dif0[n])
#         end
#     end
#     for n=1:length(f.dif0)
#         f.toincr[n] && fillsync!(f.dif0[n], 0)           # zeroed by back at every item in sequence
#     end
#     if seq
#         for w in params(f)       # zeroed only once at the beginning of the sequence
#             isdefined(w,:diff) && isdefined(w,:inc) && fillsync!(w.diff,0)
#         end
#     end
# end

# function initback0(f::Net, dy, dx...; a...)
#     N = length(f.op)
#     dxback = falses(N)
#     if !isempty(dx)                                     # todo: what if this changes after net init?
#         lastinput = 0
#         for n=1:N
#             isa(f.op[n], Input) || continue
#             dx[lastinput += 1] == nothing || (dxback[n] = true)
#         end
#     end
#     for n=N:-1:1
#         nsparse = difsparse(f, dy, n)
#         f.toback[n] || dxback[n] || continue
#         # f.dif0[n] = finddif(f, n, nsparse)  # todo-OPTIMIZATION
#         f.dif0[n] = newarray(gpu(), nsparse, eltype(f.out0[n]), size(f.out0[n]))
#         if f.toincr[n]    # todo: what if this changes after init?
#             # f.tmp[n] = findtmp(f, n, nsparse) # todo-OPTIMIZATION
#             f.tmp[n] = newarray(gpu(), nsparse, eltype(f.out0[n]), size(f.out0[n]))
#         end
#     end
# end

        # else                    # otherwise find one with matching size
        #     for k=o+1:length(f.op)
        #         if (!f.toincr[k]
        #             && (f.dif0[ k]!=nothing)
        #             && size(f.dif0[k]) == size(f.out0[n])
        #             && stype(f.dif0[k]) == st
        #             && f.outputs[k][1] > k)
        #             dif0 = f.dif0[k]
        #             # todo: However we need to check and see if this has been used between n..o
        #             # Also we don't know whether o > n for sure
        #             # This all needs more thinking
        #             break
        #         end
        #     end
# should be part of reset:
    # fillsync!(f.dif, nothing)
    # for n=1:length(f.op)
    #     (f.dif0[ n]!=nothing) && f.toincr[n] && fillsync!(f.dif0[n], 0)
    # end

# for n=length(f.op):-1:1
#     f.toback[n] || continue # we mix initback and initback0 here because if getdx or seq changes we may need to alloc new arrays.
#     if (f.dif0[n]!=nothing) # todo: implement batch size change
#         @assert (issimilar2(f.dif0[n], f.out0[n]) ) #todo: && issparse(f.dif0[n])==(f.sparse[n]!=nothing)
#     else
#         f.dif0[n] = finddif(f, n)
#         # f.dif[n] = nothing # Leave this to reset! otherwise s2s does not work
#     end
#     if f.toincr[n]
#         if (f.tmp[n]!=nothing)
#             @assert (issimilar2(f.tmp[n], f.out0[n]) && issparse(f.tmp[n])==(f.sparse[n]!=nothing))
#         else
#             f.tmp[n] = findtmp(f, n)
#             # fillsync!(f.dif0[n], 0)  # This will break s2s best if finddif zeros during alloc
#         end
#     end
# end

    # This is not true for nce; ygold and f.out0[end] have different dimensions.
    # @assert ygold == nothing || issimilar2(ygold, f.out0[end])


text:
using BenchmarkTools
using ExtendableSparse
using GFlops
using IterativeSolvers
using LinearAlgebra
using Random
using SparseArrays
using Suppressor
@suppress using MatrixDepot

include("utils/myrand.jl")
include("utils/utils.jl")

import Base: sum

# Command line arguments
# @show ARGS
nARGS = length(ARGS)
csv = ARGS[1]
solver = ARGS[2]
bench = ARGS[3]
matrix = ARGS[4]
if bench == "nthreads" || bench == "matrix"
	sizes = [parse(Int, ARGS[5])]
elseif bench == "size"
	N0 = parse(Int, ARGS[5])
	ppomag = parse(Int, ARGS[6])
	nrun = parse(Int, ARGS[7])
	sizes = get_vsizes(N0, ppomag, nrun)
else
	throw(ArgumentError(bench))
end
kernel = ARGS[nARGS-7]
format = ARGS[nARGS-6]
orth = ARGS[nARGS-5]
precon = ARGS[nARGS-4]
nthreads = parse(Int, ARGS[nARGS-3])
if nthreads != Threads.nthreads()
	throw(ArgumentError(nthreads))
end
nsamples = parse(Int, ARGS[nARGS-2])
niter = ARGS[nARGS-1]
small_tol = parse(Bool, ARGS[nARGS])

# Use packages with spmv kernel
# (Did not find a better way)
if kernel == "par"
	using MtSpMV
elseif kernel == "mkl"
	using MKLSparse
elseif kernel == "test"
	include("mymul.jl")
elseif kernel == "hist"
	using MtSpMV
elseif kernel == "flop"
	include("utils/custom_blas.jl")
end

# Benchmark parameters
BenchmarkTools.DEFAULT_PARAMETERS.samples = nsamples # 10000 by default
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 24*60*60 # 5 by default

# Returns solver function
function get_solver()
	usablesolvers = Dict(
		"cg" => cg,
	    "minres" => minres,
	    "gmres" => gmres,
	    "bicgstabl" => bicgstabl
	)
	if haskey(usablesolvers, solver)
		solverfunc = usablesolvers[solver]
	else
		throw(ArgumentError(solver))
	end
	return solverfunc
end

# Returns type of orthogonalization method (for gmres)
function get_orth()
	usableorth = Dict(
		"mgs" => ModifiedGramSchmidt(),
	    "cgs" => ClassicalGramSchmidt(),
	    "dgks" => DGKS()
	)
	if haskey(usableorth, orth)
		orth_meth = usableorth[orth]
	else
		throw(ArgumentError(orth))
	end
	return orth_meth
end

# Returns preconditioner
function get_precon()
	usableprecon = Dict(
		"jacobi" => JacobiPreconditioner,
	    "ilu0" => ILU0Preconditioner,
	    "parjacobi" => ParallelJacobiPreconditioner,
	    "parilu0" => ParallelILU0Preconditioner
	)
	if haskey(usableprecon, precon)
		preconfunc = usableprecon[precon]
	else
		throw(ArgumentError(precon))
	end
	return preconfunc
end

# Sum all types of flop of a counter
function sum(cnt::GFlops.Counter)
	cnt32 = cnt.fma32 + cnt.add32 + cnt.sub32 + cnt.mul32 + cnt.div32 + cnt.sqrt32
	cnt64 = cnt.fma64 + cnt.add64 + cnt.sub64 + cnt.mul64 + cnt.div64 + cnt.sqrt64
	return cnt32 + cnt64
end

# Use costum BLAS kernels to compute gflop
function flop_count()

	# Solver
	solverfunc = get_solver()

	# Orthogonalization method (for gmres)
	if solver == "gmres"
		orth_meth = get_orth()
	end

	for n in sizes

		# Linear system
		A, transA, b = get_linsys(matrix, n)

		# Number of iterations and tolerance
		maxiter = size(A, 2)
		if niter != "default"
		 	maxiter = parse(Int, niter)
		end
		reltol = sqrt(eps(real(eltype(b))))
		if small_tol == true
			reltol = 1e-99
		end

		# Bench
		cnt_precon = 0
		cnt_reorder = 0
		cnt_solver = 0
		if precon == "none"
			if solver == "gmres"
				cnt_iter = @count_ops solverfunc($transA, $b, maxiter=$maxiter,
					reltol=$reltol, orth_meth=$orth_meth)[2]
			else
				cnt_iter = @count_ops solverfunc($transA, $b, maxiter=$maxiter,
					reltol=$reltol)[2]
			end
		else
			if precon == "parilu0"
				Random.seed!(1)
				cnt_reorder = @count_ops reorderlinsys(A, b)
				Random.seed!(1)
				A, b = reorderlinsys(A, b)
			end
			preconfunc = get_precon()
			P = preconfunc(A)
			cnt_precon = @count_ops preconfunc($A)
			if solver == "gmres"
				cnt_iter = @count_ops solverfunc($transA, $b, Pl=$P,
					maxiter=$maxiter, reltol=$reltol, orth_meth=$orth_meth)[2]
			else
				cnt_iter = @count_ops solverfunc($transA, $b, Pl=$P,
					maxiter=$maxiter, reltol=$reltol)[2]
			end
		end

		# Write flop of solver to csv file
		csvv = replace(csv, "flop" => "flop-solver")
		s = string(A.m, " ", A.n, " ", length(A.nzval), " ", flop_solver, "\n")
		open(csvv, "a") do io
			write(io, s)
		end

		# Write flop of preconditioner creation to csv file
		if precon != "none"
			csvv = replace(csv, "flop" => "flop-precon")
			s = string(A.m, " ", A.n, " ", length(A.nzval), " ", flop_precon, "\n")
			open(csvv, "a") do io
				write(io, s)
			end
		end

		# # Write flop of reordering to csv file
		# if precon == "parilu0"
		# 	csvv = replace(csv, "flop" => "flop-reorder")
		# 	s = string(A.m, " ", A.n, " ", length(A.nzval), " ", flop_reorder, "\n")
		# 	open(csvv, "a") do io
		# 		write(io, s)
		# 	end
		# end

	end
end

# Single solve for convergence history
function conv_hist()

	# Solver
	solverfunc = get_solver()

	# Orthogonalization method (for gmres)
	if solver == "gmres"
		orth_meth = get_orth()
	end

	for n in sizes

		# Linear system
		A, transA, b = get_linsys(matrix, n)

		# Number of iterations and tolerance
		maxiter = size(A, 2)
		if niter != "default"
		 	maxiter = parse(Int, niter)
		end
		reltol = sqrt(eps(real(eltype(b))))
		if small_tol == true
			reltol = 1e-99
		end

		# Bench
		if precon == "none"
			if solver == "gmres"
				ch = solverfunc(transA, b, maxiter=maxiter, reltol=reltol,
					orth_meth=orth_meth, log=true)[2]
			else
				ch = solverfunc(transA, b, maxiter=maxiter, reltol=reltol,
					log=true)[2]
			end
		else
			if precon == "parilu0"
				Random.seed!(1)
				P = ParallelILU0Preconditioner(A)
				coloring = P.coloring
				c = collect(Iterators.flatten(coloring))
				A = A[c,:][:,c]
				b = b[c]
				transA = transpose(sparse(transpose(A)))
			else
				preconfunc = get_precon()
				P = preconfunc(A)
			end
			if solver == "gmres"
				ch = solverfunc(transA, b, Pl=P, maxiter=maxiter, reltol=reltol,
					orth_meth=orth_meth, log=true)[2]
			else
				ch = solverfunc(transA, b, Pl=P, maxiter=maxiter, reltol=reltol,
					log=true)[2]
			end
		end

		# Write results to csv file
		s = string(A.m, " ", A.n, " ", length(A.nzval), " ",
				   ch.isconverged, " ", ch.iters, " ", nrests(ch), " ", nprods(ch),
				   " ", ch[:abstol], " ", ch[:reltol], "\n")
		open(csv, "a") do io
			write(io, s)
		end

	end

end

# Run benchmarks
function run_benchmark()

	# Write header to output file
	if bench == "size"
		write(csv, "nthreads nrows ncols nnz min med mean max\n")
	end

	# Terminal output
	if bench == "size"
		print("size = ")
	end

	# Run benchmarks for all sizes
	for n in sizes

		# Terminal output
		if bench == "nthreads"
			if kernel != "hist"
				print(nthreads, "..")
			end
		elseif bench == "size"
			print(n, "..")
		end

		# Get linear system
		A, transA, b = get_linsys(matrix, n)

		# Number of iterations and tolerance
		maxiter = size(A, 2)
		if niter != "default"
		 	maxiter = parse(Int, niter)
		end
		reltol = sqrt(eps(real(eltype(b))))
		if small_tol == true
			reltol = 1e-99
		end

		# Orthogonalization method (for gmres)
		if solver == "gmres"
			orth_meth = get_orth()
		end

		# Run benchmarks and get runtimes
		solverfunc = get_solver()

		# Bench
		btimes_solver = zeros(4)
		btimes_precon = zeros(4)
		btimes_reorder = zeros(4)
		if precon == "none"
			if solver == "gmres"
				if format == "csr"
					x, btimes_solver = @mybtimes2 $solverfunc($transA, $b,
						maxiter=$maxiter, reltol=$reltol, orth_meth=$orth_meth)
				elseif format == "csc"
					x, btimes_solver = @mybtimes2 $solverfunc($A, $b,
						maxiter=$maxiter, reltol=$reltol, orth_meth=$orth_meth)
				else
					throw(ArgumentError(format))
				end
			else
				if format == "csr"
					x, btimes_solver = @mybtimes2 $solverfunc($transA, $b,
						maxiter=$maxiter, reltol=$reltol)
				elseif format == "csc"
					x, btimes_solver = @mybtimes2 $solverfunc($A, $b,
						maxiter=$maxiter, reltol=$reltol)
				else
					throw(ArgumentError(format))
				end
			end
		else
			if precon == "parilu0"
				Random.seed!(1)
				P, btimes_precon = @mybtimes2 ParallelILU0Preconditioner($A)
				coloring = P.coloring
				c = collect(Iterators.flatten(coloring))
				A = A[c,:][:,c]
				b = b[c]
				transA = transpose(sparse(transpose(A)))
			else
				preconfunc = get_precon()
				P, btimes_precon = @mybtimes2 $preconfunc($A)
			end
			if solver == "gmres"
				if format == "csr"
					x, btimes_solver = @mybtimes2 $solverfunc($transA, $b, Pl=$P,
						maxiter=$maxiter, reltol=$reltol, orth_meth=$orth_meth)
				elseif format == "csc"
					x, btimes_solver = @mybtimes2 $solverfunc($A, $b, Pl=$P,
						maxiter=$maxiter, reltol=$reltol, orth_meth=$orth_meth)
				else
					throw(ArgumentError(format))
				end
			else
				if format == "csr"
					x, btimes_solver = @mybtimes2 $solverfunc($transA, $b, Pl=$P,
						maxiter=$maxiter, reltol=$reltol)
				elseif format == "csc"
					x, btimes_solver = @mybtimes2 $solverfunc($A, $b, Pl=$P,
						maxiter=$maxiter, reltol=$reltol)
				else
					throw(ArgumentError(format))
				end
			end
		end

		# Set first column for csv files
		if bench == "matrix"
			firstcol = matrix
		else
			firstcol = nthreads
		end

		# Write runtimes of solver to csv file
		btimes = btimes_solver
		csvv = replace(csv, "runtimes" => "runtimes-solver")
		open(csvv, "a") do io
			s = string(firstcol, " ", A.m, " ", A.n, " ", length(A.nzval),
				" ", btimes[1], " ", btimes[2], " ", btimes[3], " ", btimes[4], "\n")
			write(io, s)
		end

		# Write runtimes of preconditioner creation to csv file
		if precon != "none"
			btimes = btimes_precon
			csvv = replace(csv, "runtimes" => "runtimes-precon")
			open(csvv, "a") do io
				s = string(firstcol, " ", A.m, " ", A.n, " ", length(A.nzval),
					" ", btimes[1], " ", btimes[2], " ", btimes[3], " ", btimes[4], "\n")
				write(io, s)
			end
		end

		# # Write runtimes of reordering to csv file
		# if precon == "parilu0"
		# 	btimes = btimes_reorder
		# 	csvv = replace(csv, "runtimes" => "runtimes-reorder")
		# 	open(csvv, "a") do io
		# 		s = string(nthreads, " ", A.m, " ", A.n, " ", length(A.nzval),
		# 			" ", btimes[1], " ", btimes[2], " ", btimes[3], " ", btimes[4], "\n")
		# 		write(io, s)
		# 	end
		# end

	end

	# Terminal output
	if bench == "size"
		println()
	end
end

# Run from command line
function main()
	if kernel == "hist"
		conv_hist()
	elseif kernel == "flop"
		flop_count()
	elseif any(x->x==kernel, ["ser", "par", "mkl"])
		run_benchmark()
	else
		throw(ArgumentError(kernel))
	end
end
main()

max_stars_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_issues_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_forks_repo_name: jkrch/JuliaParGMRES-BenchCLI


text:
"
    onesite dmrg
"
@with_kw struct Dmrg{F} <: Algorithm
    tol::Float64 = Defaults.tol
    maxiter::Int = Defaults.maxiter
    verbose::Bool = Defaults.verbose
    finalize::F = Defaults._finalize
end

find_groundstate(state,H,alg::Dmrg,envs...) = find_groundstate!(copy(state),H,alg,envs...)
function find_groundstate!(state::Union{FiniteMPS,MPSComoving}, H::Hamiltonian,alg::Dmrg,envs = environments(state,H))
    tol=alg.tol;maxiter=alg.maxiter
    iter = 0; delta::Float64 = 2*tol

    while iter < maxiter && delta > tol
        delta=0.0

        for pos = [1:(length(state)-1);length(state):-1:2]
            (eigvals,vecs) = @closure eigsolve(state.AC[pos],1,:LR,Arnoldi(krylovdim=60)) do x
                ac_prime(x,pos,state,envs)
            end
            delta = max(delta,calc_galerkin(state,pos,envs))

            state.AC[pos] = vecs[1]
        end

        alg.verbose && @info "Iteraton $(iter) error $(delta)"
        flush(stdout)

        iter += 1

        #finalize
        (state,envs) = alg.finalize(iter,state,H,envs)::Tuple{typeof(state),typeof(envs)};
    end

    return state,envs,delta
end

"twosite dmrg"
@with_kw struct Dmrg2{F} <: Algorithm
    tol = Defaults.tol;
    etol = 1e-6;
    maxiter = Defaults.maxiter;
    trscheme = truncerr(1e-6);
    krylovdim = 30;
    verbose = Defaults.verbose
    finalize::F = Defaults._finalize
end

find_groundstate(state,H,alg::Dmrg2,envs...) = find_groundstate!(copy(state),H,alg,envs...)
function find_groundstate!(state::Union{FiniteMPS,MPSComoving}, H::Hamiltonian,alg::Dmrg2,envs = environments(state,H))
    tol=alg.tol;maxiter=alg.maxiter
    iter = 0; delta::Float64 = 2*tol
    prev_energy = 10000000;
    energy_diff = 10000000;

    while iter < maxiter && delta > tol && energy_diff > alg.etol
        delta=0.0

        ealg = Arnoldi(krylovdim=alg.krylovdim)

        #left to right sweep
        for pos= 1:(length(state)-1)
            @tensor ac2[-1 -2; -3 -4]:=state.AC[pos][-1,-2,1]*state.AR[pos+1][1,-3,-4]
            (eigvals,vecs) = @closure eigsolve(ac2,1,:LR,ealg) do x
                ac2_prime(x,pos,state,envs)
            end
            newA2center = vecs[1]

            (al,c,ar,ϵ) = tsvd(newA2center,trunc=alg.trscheme,alg=TensorKit.SVD())
            normalize!(c);
            v = @tensor ac2[1,2,3,4]*conj(al[1,2,5])*conj(c[5,6])*conj(ar[6,3,4])
            delta = max(delta,abs(1-abs(v)));

            state.AC[pos] = (al,complex(c))
            state.AC[pos+1] = (complex(c),_permute_front(ar))
        end

        for pos = length(state)-2:-1:1
            @tensor ac2[-1 -2; -3 -4]:=state.AL[pos][-1,-2,1]*state.AC[pos+1][1,-3,-4]
            (eigvals,vecs) = @closure eigsolve(ac2,1,:LR,ealg) do x
                ac2_prime(x,pos,state,envs)
            end
            newA2center = vecs[1]

            (al,c,ar,ϵ) = tsvd(newA2center,trunc=alg.trscheme,alg=TensorKit.SVD())
            normalize!(c);
            v = @tensor ac2[1,2,3,4]*conj(al[1,2,5])*conj(c[5,6])*conj(ar[6,3,4])
            delta = max(delta,abs(1-abs(v)));

            state.AC[pos+1] = (complex(c),_permute_front(ar))
            state.AC[pos] = (al,complex(c))
        end

        #finalize
        (state,envs) = alg.finalize(iter,state,H,envs)::Tuple{typeof(state),typeof(envs)};
        energy = real(sum(expectation_value(state,H)))
        energy_diff = abs(energy - prev_energy)
        prev_energy = energy

        alg.verbose && @info "Iteraton $(iter) error $(delta) energy_diff $(energy_diff)"
        flush(stdout)
        iter += 1
    end

    return state,envs,delta
end


text:
using Test

using FluidSynth

settings = Settings()
setstr(settings, "audio.driver", "pulseaudio")

synth = Synth(settings)
sfload(synth, "violin_sample.sf2")

player = Player(synth)
audio_driver = AudioDriver(settings, synth)

add(player, "test.mid")
play(player)

FluidSynth.join(player)


text:
module CGinJulia
#
export Vector2, Point2, Point2C, Point2H, Line2
export Vector3, Point3, Point3C, Point3H, Plane3, Point3L
export toPoint2H, toPoint3H, lift
export dot, perpdot
export det2, det3, det4
export abovePl, belowPl
export plane
export ccw,leftturn,rightturn
export leftAhead, rightAhead, leftOnAhead, rightOnAhead
export leftofseg, rightofseg, leftward, rightward
export line, orientptln, leftof, rightof

export CH, ptsAbove
export Bl_Export,Export

include("LA_Dets.jl")
include("CG_2dTypes.jl")
include("CG_3dTypes.jl")
include("Export.jl")
include("Blender_Export.jl")
include("Convex_Hull.jl")

end  # module CGinJulia

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
# Data Crunching Debugger

module Dcd

	type Session
		enable::Bool
		iters::Array{Dict{String,Any}, 1}

		Session(enable = false) = ( x = new(); x.enable = enable; x.iters = Array(Dict{String,Any}, 0); return x; )
	end

	macro set(sess, key, val)
		return :(
			 $(esc(sess)).iters[end][$key] = deepcopy( $(esc(val)) )
		)
	end

	macro it(sess)
		return :(
			push!( $(esc(sess)).iters, Dict{String,Any}())
		)
	end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using CrossDecomposition
using Base.Test

# write your own tests here
include("canoncor.jl")


text:
# Given some model parameters, synthesize and plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    # default = 0
    "--norad"
    help = "Use the image alread here."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    required = true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using constants
using image
using model
using HDF5

import PyPlot.plt
using LaTeXStrings

# Plot the channel maps using sky convention
function plot_chmaps(img::image.SkyImage)

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    # CO 2-1 rest frame
    lam0 = cc/230.538e9 * 1e4 # [microns]
    nlam = length(img.lams)

    # Specify velocities

    # convert wavelengths to velocities
    vels = c_kms * (img.lams .- lam0)/lam0

    # Figure out how many plots we'll have.
    ncols = 5
    nrows = 1

    fig, ax = plt.subplots(nrows=1, ncols=ncols, figsize=(5, 1.5))

    for col=1:ncols
        iframe = col

        if iframe > nlam
            # Stop if we run out of channels
            # Plot a blank square
            ax[col][:imshow](zeros((im_ny, im_nx)), cmap=plt.get_cmap("PuBu"), vmin=0, vmax=20, extent=ext, origin="lower")

        else
            #Flip the frame for Sky convention
            frame = fliplr(img.data[:,:,iframe])
            frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
            lframe = log10(frame)
            max = maximum(lframe)
            ix,iy = ind2sub(size(lframe), indmax(lframe))
            ax[col][:imshow](lframe, extent=ext, vmin=max - 6, vmax=max, interpolation="none", origin="lower", cmap=plt.get_cmap("PuBu"))
            levels = linspace(max - 0.8, max, 8)
        end

        if col != 1
            ax[col][:xaxis][:set_ticklabels]([])
            ax[col][:yaxis][:set_ticklabels]([])
            ax[col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
        else
            ax[col][:set_xlabel](L"$\Delta \alpha$ ('')")
            ax[col][:set_ylabel](L"$\Delta \delta$ ('')")
            ax[col][:annotate](@sprintf("%.1f km/s", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
        end

    end

    fig[:subplots_adjust](hspace=0.06, wspace=0.01, top=0.98, bottom=0.13, left=0.05, right=0.95)

    plt.savefig("plots/channel_maps_sky.svg")

end

# Plot the spatially-integrated spectrum
function plot_spectrum(img::image.SkyImage)

    fig = plt.figure(figsize=(3.4,2.5))
    ax = fig[:add_subplot](111)

    spec = imToSpec(img)

    # CO 2-1 rest frame
    lam0 = cc/230.538e9 * 1e4 # [microns]

    # convert wavelengths to velocities
    vels = c_kms * (img.lams .- lam0)/lam0 # [km/s]

    ax[:plot](vels, spec[:,2], ls="steps-mid")
    # ax[:plot](vels, reverse(spec[:,2]), ls="steps-mid")
    ax[:set_ylabel](L"$f_\nu$ [Jy]")
    ax[:set_xlabel](L"$v$ [km/s]")
    ax[:set_xlim](minimum(vels), maximum(vels))

    fig[:subplots_adjust](left=0.2, bottom=0.2, right=0.8)

    plt.savefig("plots/spectrum_talk.svg")
end

# read the wavelengths for all 23 channels
fid = h5open(config["data_file"], "r")
lams = read(fid["lams"]) # [μm]
close(fid)

pp = config["parameters"]
params = ["M_star", "r_c", "T_10", "q", "gamma", "logM_CO", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)

for i=1:nparam
    starting_param[i] = pp[params[i]][1]
end

# Convert logM_CO to M_CO
starting_param[6] = 10^starting_param[6]

pars = Parameters(starting_param...)

vel = pars.vel # [km/s]
# RADMC conventions for inclination and PA
incl = pars.incl # [deg]
PA = pars.PA # [deg] Position angle runs counter clockwise
npix = config["npix"] # number of pixels

# Doppler shift the dataset wavelength to rest-frame wavelength
# beta = vel/c_kms # relativistic Doppler formula
# shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]

grd = config["grid"]
grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# Comment out afterward
# For channel maps
# nchan = 5
# global const vels = linspace(-4., 4, nchan) # [km/s]

# For spectrum
nchan = 60
global const vels = linspace(-16., 16, nchan) # [km/s]
# # CO 2-1 rest frame
lam0 = cc/230.538e9 * 1e4 # [microns]
# # convert velocities to wavelengths
shift_lams = lam0 * (vels/c_kms + 1)
#
# write_grid("", grid)
# write_model(pars, "", grid)
# write_lambda(shift_lams, "")
#
# if !parsed_args["norad"]
#     run(`radmc3d image incl $incl posang $PA npix $npix loadlambda`)
# end

im = imread()

# plot_chmaps(im)

skim = imToSky(im, pars.dpc)

# plot_chmaps(skim)
plot_spectrum(skim)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using LowRankModels
using Plotly

test_losses = Loss[
QuadLoss(),
L1Loss(),
HuberLoss(),
PeriodicLoss(1),
OrdinalHingeLoss(1,10),
WeightedHingeLoss()
LogisticLoss(),
]

#for test_iteration = 1:5
	# Create the configuration for the model (random losses)
	config = round(Int, abs(round(5*rand(length(test_losses)))));
	# config = [1,1,1,1,1,1,10]
	losses, doms = Array(Loss,1), Array(Domain,1);
	for (n,l) in zip(config, test_losses)
		for i=1:n
			push!(losses, l);
			push!(doms, l.domain);
		end
	end
	losses, doms = losses[2:end], doms[2:end]; # this is because the initialization leaves us with an #undef
	# losses = Array(Loss, 20)
	# fill!(losses, QuadLoss())
	# doms = Domain[l.domain for l in losses]

	# Make a low rank matrix as our data precursor
	m, n, true_k = 1000, length(doms), round(Int, round(length(losses)/2));
	X_real, Y_real = 2*randn(m,true_k), 2*randn(true_k,n);
	A_real = X_real*Y_real;

	# Impute over the low rank-precursor to make our heterogenous dataset
	A = impute(doms, losses, A_real);				# our imputed data

	p = Params(1, max_iter=1000, abs_tol=0.000001, min_stepsize=0.001);
	rx, ry = ZeroReg(), ZeroReg();

	skip = 5
	k0=skip
	model = GLRM(A, losses, rx, ry, k0, scale=false, offset=false);
	X_fit, Y_fit, ch = fit!(model, params=p, verbose=false);


text:
"""
    noise_uu(n::Int, lo = - 1, hi = 1)

Generate a signal consisting of `n` steps of uncorrelated uniform noise from
a uniform distribution on `[lo, hi]`.
"""
function noise_uu(n::Int; lo = - 1, hi = 1)
    u = Uniform(-lo, hi)
    rand(u, n)
end

"""
    noise_ug(n::Int; μ = 0, σ = 1)

Generate a signal consisting of `n` steps of uncorrelated Gaussian noise from
a normal distribution with mean `μ` and standard deviation `σ`.
"""
function noise_ug(n::Int; μ = 0, σ = 1)
    d = Normal(μ, σ)
    rand(d, n)
end

"""
    noise_brownian(n::Int; lo = - 1, hi = 1)
    noise_brownian(d::Distribution, n::Int)

Generate a signal consisting of `n` steps of Brownian noise, generated as
the zero-mean and unit standard deviation normalised cumulative sum of noise
generated from a uniform distribution on `[lo, hi]`. Optionally, a distribution
`d` from which to sample can be provided.

## Examples

```julia
# Based on uncorrelated uniform noise
noise_brownian(100)
noise_brownian(100, lo = -2, hi = 2)
noise_brownian(Uniform(-3, 3), 100)

# Based on uncorrelated Gaussian noise
μ, σ = 0, 2
noise_brownian(Normal(μ, σ), 100)
```
"""
function noise_brownian(n::Int; lo = - 1, hi = 1)
    u = Uniform(lo, hi)
    xs = cumsum(rand(u, n))
    (xs .- mean(xs)) ./ std(xs)
end

function noise_brownian(d::Distribution, n::Int)
    xs = cumsum(rand(d, n))
    (xs .- mean(xs)) ./ (std(xs))
end

export noise_uu, noise_ug, noise_brownian


text:
# This script provides an example of simulating changepoint data,
# using PELT to find these changepoints, and plotting the results
#
# Requires Plots package

using Changepoints, Distributions
using Plots

#########################
# Simulate changepoints #
#########################

n = 1000           # Sample size
λ = 100            # Frequencey of changes

# Sigma parameter sampled from a Uniform(2,15) distribution for each changepoint
μ, σ = 1.0, Uniform(2.0, 15.0)
# Generate a sample with from Normal distribution with changing variance
y, cps = @changepoint_sampler n λ Normal(μ, σ)

########################
# Finding changepoints #
########################

pelt_cps, cost = @PELT y Normal(μ, :?)

println("True changepoints:", cps)
println("Changepoint found by PELT:", pelt_cps)

changepoint_plot(y, pelt_cps)
gui()


text:
############################################################
## joAbstractFosterLinearOperator - outer constructors
############################################################

# common exports
export joLoosen

############################################################
## joLooseMatrix - outer constructors

"""
joLooseMatrix outer constructor

    joLooseMatrix(array::AbstractMatrix;
             DDT::DataType=eltype(array),
             RDT::DataType=promote_type(eltype(array),DDT),
             name::String="joLooseMatrix")

Look up argument names in help to joLooseMatrix type.

# Example
- joLooseMatrix(rand(4,3)) # implicit domain and range
- joLooseMatrix(rand(4,3);DDT=Float32) # implicit range
- joLooseMatrix(rand(4,3);DDT=Float32,RDT=Float64)
- joLooseMatrix(rand(4,3);name="my matrix") # adding name

# Notes
- if DDT:<Real for complex matrix then imaginary part will be neglected for transpose/adjoint operator
- if RDT:<Real for complex matrix then imaginary part will be neglected for forward/conjugate operator

"""
function joLooseMatrix(array::AbstractMatrix{EDT};
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),name::String="joLooseMatrix") where {EDT}

        (typeof(array)<:DArray || typeof(array)<:SharedArray) && @warn "Creating joLooseMatrix from non-local array like $(typeof(array)) is likely going to have adverse impact on JOLI's health. Please, avoid it."

        return joLooseMatrix{DDT,RDT}(name,size(array,1),size(array,2),
            v1->jo_convert(RDT,array*v1,false),
            v2->jo_convert(DDT,transpose(array)*v2,false),
            v3->jo_convert(DDT,adjoint(array)*v3,false),
            v4->jo_convert(RDT,conj(array)*v4,false),
            v5->jo_convert(DDT,array\v5,false),
            v6->jo_convert(RDT,transpose(array)\v6,false),
            v7->jo_convert(RDT,adjoint(array)\v7,false),
            v8->jo_convert(DDT,conj(array)\v8,false)
            )
end

joLoosen(A::joMatrix{DDT,RDT}) where {DDT,RDT} =
    joLooseMatrix{DDT,RDT}("joLoosen("*A.name*")",A.m,A.n,
        v1->A.fop(v1),
        v2->A.fop_T(v2),
        v3->A.fop_A(v3),
        v4->A.fop_C(v4),
        v5->get(A.iop)(v5),
        v6->get(A.iop_T)(v6),
        v7->get(A.iop_A)(v7),
        v8->get(A.iop_C)(v8)
        )

############################################################
## joLooseLinearFunction - outer constructors

export joLooseLinearFunctionAll, joLooseLinearFunction_T, joLooseLinearFunction_A,
       joLooseLinearFunctionFwd, joLooseLinearFunctionFwd_T, joLooseLinearFunctionFwd_A

"""
joLooseLinearFunction outer constructor

    joLooseLinearFunctionAll(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,iMVok::Bool=false,
        name::String="joLooseLinearFunctionAll")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunctionAll(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    iop::Function,iop_T::Function,iop_A::Function,iop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,iMVok::Bool=false,
    name::String="joLooseLinearFunctionAll") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,fMVok,
            iop,iop_T,iop_A,iop_C,iMVok
            )
"""
joLooseLinearFunction outer constructor

    joLooseLinearFunction_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function, iop::Function,iop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,iMVok::Bool=false,
        name::String="joLooseLinearFunction_T")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunction_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function, iop::Function,iop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,iMVok::Bool=false,
    name::String="joLooseLinearFunction_T") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,
            fop_T,
            v3->conj(fop_T(conj(v3))),
            v4->conj(fop(conj(v4))),
            fMVok,
            iop,
            iop_T,
            v7->conj(iop_T(conj(v7))),
            v8->conj(iop(conj(v8))),
            iMVok
            )
"""
joLooseLinearFunction outer constructor

    joLooseLinearFunction_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function, iop::Function,iop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,iMVok::Bool=false,
        name::String="joLooseLinearFunction_A")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunction_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function, iop::Function,iop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,iMVok::Bool=false,
    name::String="joLooseLinearFunction_A") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,
            v2->conj(fop_A(conj(v2))),
            fop_A,
            v4->conj(fop(conj(v4))),
            fMVok,
            iop,
            v6->conj(iop_A(conj(v6))),
            iop_A,
            v8->conj(iop(conj(v8))),
            iMVok
            )
"""
joLooseLinearFunction outer constructor

    joLooseLinearFunctionFwd(m::Integer,n::Integer,
        fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,
        name::String="joLooseLinearFunctionAll")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunctionFwd(m::Integer,n::Integer,
    fop::Function,fop_T::Function,fop_A::Function,fop_C::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,
    name::String="joLooseLinearFunctionAll") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,fop_T,fop_A,fop_C,fMVok,
            @joNF, @joNF, @joNF, @joNF, false
            )
"""
joLooseLinearFunction outer constructor

    joLooseLinearFunctionFwd_T(m::Integer,n::Integer,
        fop::Function,fop_T::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,
        name::String="joLooseLinearFunctionFwd_T")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunctionFwd_T(m::Integer,n::Integer,
    fop::Function,fop_T::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,
    name::String="joLooseLinearFunctionFwd_T") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,
            fop_T,
            v3->conj(fop_T(conj(v3))),
            v4->conj(fop(conj(v4))),
            fMVok,
            @joNF, @joNF, @joNF, @joNF, false
            )
"""
joLooseLinearFunction outer constructor

    joLooseLinearFunctionFwd_A(m::Integer,n::Integer,
        fop::Function,fop_A::Function,
        DDT::DataType,RDT::DataType=DDT;
        fMVok::Bool=false,
        name::String="joLooseLinearFunctionFwd_A")

Look up argument names in help to joLooseLinearFunction type.

# Notes
- the developer is responsible for ensuring that used functions take/return correct DDT/RDT

"""
joLooseLinearFunctionFwd_A(m::Integer,n::Integer,
    fop::Function,fop_A::Function,
    DDT::DataType,RDT::DataType=DDT;
    fMVok::Bool=false,
    name::String="joLooseLinearFunctionFwd_A") =
        joLooseLinearFunction{DDT,RDT}(name,m,n,
            fop,
            v2->conj(fop_A(conj(v2))),
            fop_A,
            v4->conj(fop(conj(v4))),
            fMVok,
            @joNF, @joNF, @joNF, @joNF, false
            )

joLoosen(A::joLinearFunction{DDT,RDT}) where {DDT,RDT} =
    joLooseLinearFunction{DDT,RDT}("joLoosen("*A.name*")",A.m,A.n,
        v1->A.fop(v1),
        v2->get(A.fop_T)(v2),
        v3->get(A.fop_A)(v3),
        v4->get(A.fop_C)(v4),
        A.fMVok,
        v5->get(A.iop)(v5),
        v6->get(A.iop_T)(v6),
        v7->get(A.iop_A)(v7),
        v8->get(A.iop_C)(v8),
        A.iMVok
        )


text:
using SQLiteFs
using Test

@testset "SQLiteFs.jl" begin
    include("cstruct.jl")
    include("layout.jl")
end


text: include("backends/imagetile.jl")


text:
@testset "Flux variability analysis" begin
    cp = test_simpleLP()
    optimizer = Tulip.Optimizer
    fluxes = flux_variability_analysis(cp, optimizer)

    @test size(fluxes) == (2, 2)
    @test fluxes ≈ [
        1.0 1.0
        2.0 2.0
    ]

    fluxes = flux_variability_analysis(cp, [2], optimizer)

    @test size(fluxes) == (1, 2)
    @test isapprox(fluxes, [2 2], atol = TEST_TOLERANCE)

    # a special testcase for slightly sub-optimal FVA (gamma<1)
    cp = CoreModel(
        [-1.0 -1.0 -1.0],
        [0.0],
        [1.0, 0.0, 0.0],
        [0.0, 0.0, -1.0],
        1.0 * ones(3),
        ["r$x" for x = 1:3],
        ["m1"],
    )
    fluxes = flux_variability_analysis(cp, optimizer)
    @test isapprox(
        fluxes,
        [
            1.0 1.0
            0.0 0.0
            -1.0 -1.0
        ],
        atol = TEST_TOLERANCE,
    )
    fluxes = flux_variability_analysis(cp, optimizer; bounds = gamma_bounds(0.5))
    @test isapprox(
        fluxes,
        [
            0.5 1.0
            0.0 0.5
            -1.0 -0.5
        ],
        atol = TEST_TOLERANCE,
    )
    fluxes = flux_variability_analysis(cp, optimizer; bounds = _ -> (0, Inf))
    @test isapprox(
        fluxes,
        [
            0.0 1.0
            0.0 1.0
            -1.0 0.0
        ],
        atol = TEST_TOLERANCE,
    )

    @test isempty(flux_variability_analysis(cp, Vector{Int}(), Tulip.Optimizer))
    @test_throws DomainError flux_variability_analysis(cp, [-1], Tulip.Optimizer)
    @test_throws DomainError flux_variability_analysis(cp, [99999999], Tulip.Optimizer)
end

@testset "Parallel FVA" begin
    cp = test_simpleLP()

    fluxes = flux_variability_analysis(cp, [1, 2], Tulip.Optimizer; workers = W)
    @test isapprox(
        fluxes,
        [
            1.0 1.0
            2.0 2.0
        ],
        atol = TEST_TOLERANCE,
    )
end

@testset "Flux variability analysis with StandardModel" begin
    model = load_model(StandardModel, model_paths["e_coli_core.json"])
    mins, maxs = flux_variability_analysis_dict(
        model,
        Tulip.Optimizer;
        bounds = objective_bounds(0.99),
        modifications = [
            change_optimizer_attribute("IPM_IterationsLimit", 500),
            change_constraint("EX_glc__D_e", -10, -10),
            change_constraint("EX_o2_e", 0.0, 0.0),
        ],
    )

    @test isapprox(maxs["EX_ac_e"]["EX_ac_e"], 8.5185494, atol = TEST_TOLERANCE)
    @test isapprox(mins["EX_ac_e"]["EX_ac_e"], 7.4483887, atol = TEST_TOLERANCE)
end


text:
abstract type AbstractBoosting <: AbstractPolicy end

# function boost!(boosting::AbstractBoosting,gprev::AbstractVector,gcurr::AbstractVector)
#     error("No defined boost! function for boosting policy ", typeof(boosting))
# end

function boost!() end

function boost_wrapper(wid::Cint,klocal::Cint,kglobal::Cint,gprev_b::Ptr{Cdouble},gprev_e::Ptr{Cdouble},gcurr_b::Ptr{Cdouble},boost_data::Ptr{Cvoid})
    boost_policy = unsafe_pointer_to_objref(boost_data)::AbstractBoosting
    ptrdiff = Int(gprev_e - gprev_b)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    gprev = unsafe_wrap(Array, gprev_b, N)
    gcurr = unsafe_wrap(Array, gcurr_b, N)
    boost!(boost_policy, wid, klocal, kglobal, gprev, gcurr)
    return gcurr_b + ptrdiff
end

module Boosting

using Parameters
using POLO: AbstractBoosting, AbstractPolicyParameters
import POLO: initialize!, boost!

function initialize!(policy::AbstractBoosting,x₀::Vector{Float64})
    nothing
end

include("none.jl")
include("momentum.jl")

end


text:
module Senuk
    export
    Solvers,
    Regression

    include("./printutil.jl")
    include("./solvers.jl")
    include("./regression.jl")

end

max_stars_repo_name: upul/Senuk

max_issues_repo_name: upul/Senuk

max_forks_repo_name: upul/Senuk


text:
module CamJulia
export IdrisList, Runtime

export IdrisIntegerT
IdrisIntegerT = Int64

include("IdrisList.jl")
include("IdrisForeignCollections.jl")
include("Runtime.jl")
include("CommonAbstractMachine.jl")
include("ReadIR.jl")

macro load_cam(path)
    aeson = ReadIR.load_aeson(path);
    ir = ReadIR.aeson_to_ir(aeson)
    x = CAM.ir_to_julia(ir)
    esc(x)
end

end # module

max_stars_repo_name: thautwarm/idris-cam

max_issues_repo_name: thautwarm/idris-cam

max_forks_repo_name: thautwarm/idris-cam


text:
# This file was generated, do not modify it. # hide
data_repo = "https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data"

url_power   = data_repo * "/power_syst/DE_power_hourly.csv"
url_weather = data_repo * "/power_syst/DE_weather_data.csv"

power   = DataFrame(urldownload(url_power))
weather = DataFrame(urldownload(url_weather));


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

# slightly optimized version of DataStructure.jl priority queue
# swapped index field from a dict to an array
# removed checking for presence of keys when inserting

# Binary heap indexing
heapleft(i::Int)   = i << 1         #2i
heapright(i::Int)  = (i << 1) + 1   #2i + 1
heapparent(i::Int) = i >> 1         #div(i, 2)

using StaticArrays
mutable struct ArrayPQ{K,V,O<:Ordering} <: AbstractDict{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Vector{Pair{K,V}}
    o::O

    # Map elements to their index in xs
    index::Vector{Int}

    function ArrayPQ{K,V,O}(o::O, itr) where {K,V,O<:Ordering}
        xs    = Vector{Pair{K,V}}(undef, length(itr))
        index = Vector{Int}(length(itr))
        @inbounds for (i,p) in enumerate(itr)
            xs[i]    = p
            index[i] = i
        end
        pq = new{K,V,O}(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        pq
    end
end
ArrayPQ(kv::Vector{T}, o::Ordering=Forward) where {K,V,T <: Pair{K,V}} = ArrayPQ{K,V,typeof(o)}(o, kv)
ArrayPQ(kv, o::Ordering=Forward) = ArrayPQ(o, kv)

length(pq::ArrayPQ)  = length(pq.xs)
isempty(pq::ArrayPQ) = isempty(pq.xs)
peek(pq::ArrayPQ)    = (@inbounds return pq.xs[1])

function percolate_down!(pq::ArrayPQ, i::Integer)
    @inbounds x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r
        if lt(pq.o, pq.xs[j].second, x.second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    @inbounds pq.index[x.first] = i
    @inbounds pq.xs[i] = x
    nothing
end

function percolate_up!(pq::ArrayPQ, i::Integer)
    @inbounds x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x.second, pq.xs[j].second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    @inbounds pq.index[x.first] = i
    @inbounds pq.xs[i] = x
    nothing
end

function getindex(pq::ArrayPQ{K,V}, key) where {K,V}
    @inbounds return pq.xs[pq.index[key]].second
end

# Change the priority of an existing element
# does not check if the element is present
function setindex!(pq::ArrayPQ{K, V}, value, key) where {K,V}
    @inbounds i = pq.index[key]
    @inbounds oldvalue = pq.xs[i].second
    @inbounds pq.xs[i] = Pair{K,V}(key, value)
    if lt(pq.o, oldvalue, value)
        percolate_down!(pq, i)
    else
        percolate_up!(pq, i)
    end
    value
end

# Unordered iteration through key value pairs in a ArrayPQ
start(pq::ArrayPQ) = start(pq.index)
done(pq::ArrayPQ, i) = done(pq.index, i)
function next(pq::ArrayPQ{K,V}, i) where {K,V}
    (k, idx), i = next(pq.index, i)
    return (pq.xs[idx], i)
end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
# run `export JULIA_LOAD_PATH=.:` to enable loading of LinearRegression module

# using Plots
using LinearRegression

import Random

Random.seed!(0)

# optimal solution using a linear model for comparison

x_b = [1, 2, 3]
y_b = [1, 1, 2]

function fitline(x::Vector{<:Number}, y::Vector{<:Number})::Function
    n = size(x)[1]
    a = (sum(y) * sum(x.^2) - sum(x) * sum(x .* y)) / (n * sum(x.^2) - sum(x)^2) # t
    b = (n * sum(x .* y) - sum(x) * sum(y)) / (n * sum(x.^2) - sum(x)^2)# m
    x->b * x + a
end

optimal_linear_model = fitline(x_b, y_b)

# end of optimal solution

# a few different basis functions for testing
Φ1(j, 𝐱) =
    if j == 0
    1
else
    𝐱[1]
end

Φ2(j, 𝐱) = 𝐱[1]^j
Φ3(j, 𝐱) = sin(1 / j * 𝐱[1])
σ(a) = 1 / (1 + exp(-a))
function Φ4(j, 𝐱)
    μ = 0.2
    s = 0.2
    σ((𝐱[1] - μ) / s)
end

"""
# test1

𝐗 = [[1], [2], [3]]
t = [1, 1, 2]

model1 = gd(Φ1, 𝐗, t, 0.001, 2, 20000, 10e-12)

x = 1:0.1:5

p = scatter(map(x->x[1], 𝐗), t, label = "training");
plot!(x, model1.(map(x->[x], x)), label = "prediction")
plot!(x, optimal_linear_model.(x), label = "optimal", line = :dot)
display(p)
readline()
"""

"""
# test2

𝐗 = [[0], [1], [2], [3], [4], [5]]
t = [0, 1, 4, 9, 16, 25]
t += randn(size(t)[1]) * 3

model1 = gd(Φ3, 𝐗, t, 0.00001, 5, 200000, 10e-12)

x = 0:0.1:5

p = scatter(map(x->x[1], 𝐗), t, label = "training");
plot!(x, model1.(map(x->[x], x)), label = "prediction")
display(p)
readline()
"""

# test3
𝐗 = hcat([0; 1; 2; 3; 4; 5]) # hcat to convert to matrix because julia is weird like that
t = [0, 1, 4, 9, 16, 25]
t += randn(size(t)[1]) * 3

for i = 0:100
    (model1, residual_error) = @time fit_linear_model(Φ2, 𝐗, t, 0.00005, 2, 300000, 10e-3, 0.9, 0) # Φ2, 𝐗, t, 0.00005, 2, 300000, 10e-5, 0.2)
end
x = 0:0.1:5

# p = scatter(map(x->x[1], 𝐗), t, label = "training");
# plot!(x, model1.([[x] for x in x]), label = "prediction")
# display(p)
# readline()

"""
# test4
using StatsBase

f(x) = sin(3*x + 10)
x = -1:0.1:2
y = f.(x)
𝐗 = hcat([x for x in sample(x, 20)]) # hcat to convert to matrix because julia is weird like that
t = f.(𝐗)[:]
t += randn(size(t)[1]) * 0.3

(model1, residual_error) = fit_linear_model(Φ2, 𝐗, t, 0.00008, 5, 200000, 10e-8, 0.95, 0.008)

x = -1:0.1:2

p = scatter(map(x->x[1], 𝐗), t, label = "training");
plot!(x, model1.([[x] for x in x]), label = "prediction")
plot!(x, y, label = "real")
display(p)
readline()
"""


text:
# Run Nested Sampling on the model imported from here
include("models/SpikeSlab.jl")

# Tuning parameters
num_particles = 100
mcmc_steps = 1000

# Depth in nats
max_depth = 100.0

# Do an NS run
include("Sampler.jl")
do_nested_sampling(num_particles, mcmc_steps, max_depth)


text:
## HMM.jl, contains constructors for HMM objects and 'generate()' function

# HMM -- HMM with emission probability distribution 'C'
type HMM{C<:Distribution}   <: StatisticalModel
	n::Int             # Number of hidden states
	A::Matrix{Float64} # Estimated state-transition matrix A[i,j] = Pr[i->j]
	B::Vector{C}       # Estimated emission probability distributions
	p::Vector{Float64} # Estimiated initial state probabilities

	# Notes:
	#   "A" is a NxN matrix, rows sum to one

	# To do:
	#    Allow B to depend on other observables, for observation o and param k, B(o|k)
end

# TO DO: add a constructor each B[i] is drawn from some distribution specified by hyperparams

function HMM(n::Int,C::Distribution)
	# Randomize state-transition matrix
	A = rand(n,n)
	A ./= sum(A,2) # normalize rows

	# Specify a distribution of type C for each state
	B = (Distribution)[]
	for i = 1:n
		push!(B,deepcopy(C))
	end

	# Randomize initial state probabilities
	p = rand(n)
	p ./= sum(p)

	return HMM(n,A,B,p)
end

function HMM(A::Matrix{Float64},C::Distribution)
	# determine number of states
	@assert size(A,1) == size(A,2)
	n = size(A,1)

	# Specify a distribution of type C for each state
	B = (Distribution)[]
	for i = 1:n
		push!(B,deepcopy(C))
	end

	# Randomize initial state probabilities
	p = rand(n)
	p ./= sum(p)

	return HMM(n,A,B,p)
end

function HMM(A::Matrix{Float64},Bmat::Matrix{Float64})
	# Initialize a discrete HMM, with Categorical() emissions
	n = size(A,1)

	# Each row of Bmat specifies a Categorical pdf
	B = (Distribution)[]
	for i = 1:n
		push!(B,Categorical(vec(Bmat[i,:])))
	end

	return HMM(A,B)
end

function HMM(A::Matrix{Float64},B::Vector{Distribution})
	# determine number of states
	@assert size(A,1) == size(A,2)
	n = size(A,1)
	@assert length(B) == n

	# Randomize initial state probabilities
	p = rand(n)
	p ./= sum(p)

	return HMM(n,A,B,p)
end

function HMM(A::Matrix{Float64},B,p::Vector{Float64})
	@assert sum(p)==1
	hmm = HMM(A,B)
	hmm.p = p # reset initial probabilities
	return hmm
end

function rand(hmm::HMM, n_obs::Int)
	# Generate a sequence of n_obs observations from an HMM.

	# Sequence of states and observations
	s = zeros(Int,n_obs) # states
	o = zeros(n_obs)     # observations

	# Choose initial state with probabilities weighted by "init_state"
	s[1] = rand(Categorical(hmm.p))  # hmm.p are the initial state probabilities
	o[1] = rand(hmm.B[s[1]])         # draw observation given initial state

	# Construct categorical distributions from each row of A
	Ac = (Categorical)[]
	for i = 1:hmm.n
		push!(Ac,Categorical(vec(hmm.A[i,:])))
	end

	# Iterate drawing observations and updating state
	for t = 2:n_obs
		s[t] = rand(Ac[s[t-1]])   # sample from appropriate row of state-transition matrix
		o[t] = rand(hmm.B[s[t]])  # sample from emission probability distribution
	end

	# return sequence of states and observations
	return (s,o)
end


text:
function virial(A::AbstractSystem, p::EmpiricalPotential)
    v = virial_stress(A, p)
    v[1] + v[2] + v[3]
end

function virial_stress(A::AbstractSystem, p::EmpiricalPotential)
    v = zeros(6)
    nnlist = neighborlist(A, p.rcutoff)

    for ii in 1:length(A)
        r = nnlist.r[ii]
        R = nnlist.R[ii]
        for (rj, Rj) in zip(r, R)
            fo = force(Rj, rj, p)
            vi = rj * fo'
            v += [vi[1, 1], vi[2, 2], vi[3, 3], vi[3, 2], vi[3, 1], vi[2, 1]]
        end
    end
    SVector{6}(v)
end


text:
module HealthHandler

import ..State
import HTTP

struct Handler
    state::State.SharedState
end

function (h::Handler)(req::HTTP.Request)
    if h.state.healthy
        HTTP.Response(200)
    else
        HTTP.Response(500)
    end
end

end # module

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
module Janitor

# Write your package code here.

using DataFrames

export describe_class

include("describe_class.jl")

end