text:
using HorizonSideRobots

function move_to_border!(r, s)
    while !isborder(r, s)
        move!(r, s)
    end
end

function move_to_side(s, num)
    for i in 1:num
        move!(r, s)
    end
end

function move_mark!(r, s)
    while !isborder(r, s)
        move!(r, s)
        putmarker!(r)
    end
end

function move_by_marks!(r, s)
    while ismarker(r)
        move!(r, s)
    end
end

function inverse_side(s)
    HorizonSide((Int(s) + 2) % 4)
end

function rotate_left(s)
    HorizonSide((Int(s) + 1) % 4)
end

function rotate_right(s)
    return HorizonSide((Int(s) + 3) % 4)
end

function mark_move!(r, s)
    while !isborder(r, s)
        putmarker!(r)
        move!(r, s)
    end
end

function go_to_corner!(r, side1, side2)
    h = 0
    w = 0
    while !isborder(r, side1) || !isborder(r, side2)
        if !isborder(r, side1)
            move!(r, side1)
            h += 1
        elseif !isborder(r, side2)
            move!(r, side2)
            w += 1
        end
    end
    return (h, w)
end

sides = []
function go_to_corner_with_walls!(r, side1, side2)
    while !isborder(r, side1) || !isborder(r, side2)
        while !isborder(r, side1)
            move!(r, side1)
            push!(sides, side1)
        end
        while !isborder(r, side2)
            move!(r, side2)
            push!(sides, side2)
        end
    end
end

function putmarker_to_corner(side)
    putmarker!(r)
    while !isborder(r, side)
        move!(r, side)
    end
end

function move_some_steps!(r, s, steps)
    for i in 1:steps
        move!(r, s)
    end
end

function move_mark_some_steps!(r, s, steps)
    putmarker!(r)
    for i in 1:steps
        move!(r, s)
        putmarker!(r)
    end
end

function move_back(r, side1, side2, steps1, steps2)
    move_some_steps!(r, side1, steps1)
    move_some_steps!(r, side2, steps2)
end

function mark_line!(r, s)
    steps = 0
    putmarker!(r)
    while !isborder(r, s)
        move!(r, s)
        putmarker!(r)
        steps += 1
    end
    return steps
end

function mark_snake!(r, s)
    while true
        mark_line!(r, s)
        if isborder(r, Nord)
            break
        end
        move!(r, Nord)
        s = inverse_side(s)
    end
end

function line(side, flag)
    while !isborder(r, side)
        if flag == 1
            putmarker!(r)
            move!(r, side)
            if (!isborder(r, side))
                move!(r, side)
            end
        else
            move!(r, side)
            putmarker!(r)
            if (!isborder(r, side))
                move!(r, side)
            end
        end
    end
    if (flag == 1)
        putmarker!(r)
    end
end

function next(s, num)
    return HorizonSide((Int(s) + num) % 4)
end

function print_line(s, num)
    while !isborder(r, s)
        for i in 1:num
            putmarker!(r)
            if !isborder(r, s)
                move!(r, s)
            end
        end
        for i in 1:num
            if !isborder(r, s)
                move!(r, s)
            end
        end
    end
end

function print_line(s, n, count)
    while !isborder(r, s)
        if trunc(count / n) % 2 == 0
            putmarker!(r)
        end
        move!(r, s)
        count += 1
    end
    if trunc(count / n) % 2 == 0
        putmarker!(r)
    end
    count += 1
    while !isborder(r, inverse_side(s))
        move!(r, inverse_side(s))
    end
end

function go_around(s, dir::AbstractArray)
    num = 0
    flag = true
    while (isborder(r, s) && !isborder(r, next(s, 1)))
        move!(r, next(s, 1))
        num += 1
        push!(dir, next(s, 1))
    end
    if !isborder(r, s)
        move!(r, s)
        push!(dir, s)
    else
        flag = false
    end
    if flag
        while isborder(r, next(s, 3))
            move!(r, s)
            push!(dir, s)
        end
        for i in 1:num
            move!(r, next(s, 3))
            push!(dir, (next(s, 3)))
        end
    end
end

function go_to_border_and_draw(s)
    dir = []
    while true
        if !isborder(r, s)
            putmarker!(r)
            move!(r, s)
            push!(dir, s)
        else
            putmarker!(r)
            go_around(s, dir)
        end
        if (isborder(r, s) && isborder(r, next(s, 1)))
            break
        end
    end

    while length(dir) > 0
        move!(r, next(pop!(dir), 2))
    end
end

function go_forward(s)
    while !isborder(r, s)
        putmarker!(r)
        move!(r, s)
    end
end

function set(side1, side2, dir)
    while !isborder(r, side1) || !isborder(r, side2)
        while !isborder(r, side1)
            move!(r, side1)
            push!(dir, side1)
        end
        while !isborder(r, side2)
            move!(r, side2)
            push!(dir, side2)
        end
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
using Test, Random, Zygote, Distributions, SpecialFunctions, FiniteDifferences
import TruncatedNormal as TN

for a = -10:10
    d = truncated(Normal(), a, Inf)
    @test TN.tnmean(a) ≈ mean(d)
    @test TN.tnstd(a)  ≈ std(d)
    @test TN.tnvar(a)  ≈ var(d)
    @test a ≤ TN.tnmean(a) < Inf
    @test 0 ≤ TN.tnvar(a) ≤ 1
end
@test 1e20 ≤ TN.tnmean(1e20) < Inf
@test_broken 0 ≤ TN.tnvar(1e80) ≤ 1

for a = -2:2
    @test mean(TN.randnt'(a) for _ = 1:10^6) ≈ TN.tnmean'(a) rtol=0.1
    z, da = TN.∇randnt(Random.GLOBAL_RNG, a)
    @test da ≈ erfc(z/√2) / erfc(a/√2) * exp((z^2 - a^2)/2)
end

@testset "randnt_half gradient" begin
    for μ = -1:1, σ = 1:2
        dμ_ = mean(gradient(TN.randnt_half, μ, σ)[1] for _ = 1:10^6)
        dσ_ = mean(gradient(TN.randnt_half, μ, σ)[2] for _ = 1:10^6)
        dμ, dσ = gradient(μ, σ) do μ, σ
            μ + σ * TN.tnmean(-μ/σ)
        end
        @test dμ ≈ dμ_ rtol=0.1
        @test dσ ≈ dσ_ rtol=0.1
    end
end

@testset "sqrt1half" begin
    @test TN.sqrt1half(5) ≈ 5.1925824035672520156
    @test TN.sqrt1half(0) == 1
    @test TN.sqrt1half(-1) == TN.sqrt1half(1) ≈ 1.6180339887498948482
    @test isnan(TN.sqrt1half(NaN))
    @test TN.sqrt1half(Inf) == TN.sqrt1half(-Inf) == Inf
    @test TN.sqrt1half(1e300) ≈ 1e300
end

@testset "randnt" begin
    @test TN.randnt(0) > 0
    @test TN.randnt(1e300) == 1e300
    @test TN.randnt(Inf) == Inf
    @test isnan(TN.randnt(NaN))
    @test TN.randnt(floatmax(Float64)) == floatmax(Float64)
end

ext: jl


text:
"""
This script finds the probability that a new user
pair will find a path in a grid lattice with
a number of other competing end-users.
"""

using LightGraphs
using QuNet
using Statistics
using Plots

# Make grid network
# Get a list of random users
# For each user
    # Find the shortest path between them
        # If it exists, remove it, add a success to the user tally
        # If it doesn't, add a failure to the user tally
# Repeat for another set of end-users

function make_userpairs(num_nodes::Int, num_pairs::Int; node_list=nothing)::Vector{Tuple{Int64, Int64}}
    @assert num_nodes >= num_pairs*2 "Graph space too small for number of pairs"
    if node_list != nothing
        rand_space = node_list
    else
        rand_space = Array(collect(1:num_nodes))
    end
    pairs = Vector{Tuple}()
    i = 0
    while i < num_pairs
        idx = rand(1:length(rand_space))
        u = rand_space[idx]
        deleteat!(rand_space, idx)
        idx = rand(1:length(rand_space))
        v = rand_space[idx]
        deleteat!(rand_space, idx)
        chosen_pair = (u, v)
        push!(pairs, chosen_pair)
        i += 1
    end
    return pairs
end

NO_PATH_FOUND = 0
PATH_FOUND = 1

# Params
numPairs = 50
numSamples = 50000
dimX = 10; dimY = 10
makeNewData = true

# MAIN
if makeNewData == false
    pathStats = [[] for i in collect(1:numPairs)]

    for sample in 1:numSamples
        g = LightGraphs.grid([dimX, dimY])
        userPairs = make_userpairs(nv(g), numPairs)
        for (pairidx, pair) in enumerate(userPairs)
            shortestPath = a_star(g, pair[1], pair[2])
            if length(shortestPath) == 0
                push!(pathStats[pairidx], NO_PATH_FOUND)
            else
                push!(pathStats[pairidx], PATH_FOUND)
            end
            for edge in shortestPath
                rem_edge!(g, edge)
            end
        end
    end

    pathProbs = []

    for pairNumData in pathStats
        pathRate = mean(pairNumData)
        push!(pathProbs, pathRate)
    end
end

x = collect(1:numPairs)
plot(x, pathProbs)
title!("Probability nth user finds path in 10 x 10 grid")
xlabel!("Pair number")
ylabel!("Probability of finding path")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
using PyCall;

try
  pyimport("galgebra")
catch e
  run(PyCall.python_cmd(`-m pip install -e git+https://github.com/pygae/galgebra.git#egg=galgebra --user`))
  run(PyCall.python_cmd(`-m pip install sympy==1.3 --user`))
end


text:
"""
    interpret_p(p::Number; alpha=::Number=0.05)

Interpret the *p value* based on the `alpha` level.

P-values are the probability of obtaining an effect at least as extreme as the one in your sample data, assuming the truth of the null hypothesis. The significance level, also denoted as alpha or α, is the probability of rejecting the null hypothesis when it is true (*i.e.*, claiming "there is an effect" when there is not). For example, a significance level of 0.05 indicates a 5% risk of concluding that a difference exists when there is no actual difference. The traditional default for `alpha` is `.05`. However, a strong wave of criticism suggests to either **justify your alpha** ([Lakens et al., 2018](https://psyarxiv.com/9s3y6)) or lower the treshold (for instance to .005; [Benjamin et al., 2018](https://psyarxiv.com/mky9j/)).

# Arguments
- `p`: The *p* value.
- `alpha`: Significance treshold.

# Examples
```jldoctest
interpret_p(0.04)

# output

"significant"
```
"""
function interpret_p(p::Number; alpha::Number=0.05)

    rules = Rules(
        [alpha],
        ["significant", "not significant"],
        true
        )

    interpretation = interpret(p, rules)
    return interpretation
end

"""
    format_p(p::Number; stars::Bool=false)

Format the *p* value according to APA standards.

# Arguments
- `p`: The *p* value.
- `stars`: Add stars (*) when significant.

# Examples
```jldoctest
format_p(0.04, stars=true)

# output

"p < .05*"
```
"""
function format_p(p::Number; stars::Bool=false)
    if p < 0.001
        output = "p < .001"
        if stars == true
            output *= "***"
        end
    elseif p < 0.01
        output = "p < .01"
        if stars == true
            output *= "**"
        end
    elseif p < 0.05
        output = "p < .05"
        if stars == true
            output *= "*"
        end
    elseif p > 0.1
        output = "p > .1"
    else
        output = replace("p = $(round(p, digits=2))", "0." => ".")
        if stars == true
            output *= "°"
        end
    end
    return output
end


text:
for n = 1:10
-2 print(3) end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function reorder(s1::AbstractString, s2::AbstractString)
    len1 = length(s1)
    len2 = length(s2)
    if len2 > len1
        return s2, len2, s1, len1
    else
        return s1, len1, s2, len2
    end
end


text:
using NLopt
using Roots
#########################################
#computes predicted probabilities for multinomial regression, given X, beta, and number of categories
function multProb(X,beta,k)
  p=size(X)[2]
  n=size(X)[1]
  beta=reshape(beta,p,k)
  denominator=zeros(n)
  numerator=exp(X*beta)
  denominator=sum(numerator,2)
  prob=numerator./denominator
  beta=vec(beta)
  return(prob)
end

##################################
#computes predicted probabilities for multinomial regression, given X*beta
function multProbPrecompute(Xbeta)
  numerator=exp(Xbeta)
  denominator=sum(numerator,2)
  prob=numerator./denominator
  return(prob)
end

####################################
#computes gradient for multinomial regression
function multinomialGradient(X,beta,z)
  k=size(z)[2]
  grad=X'*(z-multProb(X,beta,k))
  grad=vec(grad)
  return(grad)
end

######################################
#computes multinomial log likelihood given X, beta, z
function loglikelihood(X,beta,z)
  p=size(X)[2]
  k=size(z)[2]
  beta=reshape(beta,p,k)
  probs=multProb(X,beta,k)
  val=0
  for i = 1:(size(X)[1])
    val=val+log(probs[i,find(z[i,:].==1)])
  end
  beta=vec(beta)
  return(val)
end

#######################################
#computes loglikelihood, given a precomputed X*beta and z
function loglikelihoodPrecompute(Xbeta,z)
  probs=multProbPrecompute(Xbeta)
  val=0
  for i = 1:(size(Xbeta)[1])
    val=val+log(probs[i,find(z[i,:].==1)])
  end
  return(val)
end

########################################
#not really used, does multinomial regression starting at beta, given model matrix X, response matrix z
function multinomialOptim(X,beta,z)
  function loglikeClosure(beta)
    return loglikelihood(X,beta,z)[1]
  end
  function gradientClosure(beta)
    return(multinomialGradient(X,beta,z))
  end
  opt = NLopt.Opt(:LD_LBFGS, length(beta))
  function f(beta,grad)
    grad[:]=-gradientClosure(beta)
    val=-loglikeClosure(beta)
    return(val)
  end
  NLopt.min_objective!(opt, f)
  (optf,optx,ret) = NLopt.optimize(opt,vec(beta))
  k=size(z)[2]
  p=size(X)[2]
  beta=reshape(optx,p,k)
  return(beta)
end

#######################################
#compute hessian for symmetric multinomial regression parameterization
function multinomialHessian(X,Xbeta)
  k=size(Xbeta)[2]
  p=size(X)[2]
  n=size(Xbeta)[1]

  mu=multProbPrecompute(Xbeta)
  hess=zeros(p*k,p*k)
  for i = 1:k
    start1=(i-1)*p+1
    end1=i*p

    for j = i:k
      start2=(j-1)*p+1
      end2=j*p
      temp=zeros(p)
      if i==j
        w1=mu[:,i]-mu[:,i].^2
        temp=X'*(scale(-w1,X))
      else
        w1=-mu[:,j].*mu[:,i]
        temp=X'*(scale(-w1,X))
      end
      hess[start1:end1,start2:end2]=temp
      hess[start2:end2,start1:end1]=temp
    end
  end
  return(hess)
end

#########################################
#main function, solve multinomial group lasso problem over a sequence of lambda using block coordinate descent
function multinomialBCD(X,z,groups,penaltyFactor)
  #convergence tolerance
  tol=10.0^-8
  #compute starting point for lambda
  #get number of groups
  categories=unique(groups)
  g=length(categories)
  k=size(z)[2]
  p=size(X)[2]
  n=size(z)[1]
  #get number of coefficients in each group
  groupLengths=zeros(g)
  for i = 1:g
    groupLengths[i]=sum(groups.==categories[i])
  end
#
#   #starting beta, get intercept term since this is unpenalized
  beta=zeros(p,k)
  factoredResponse=zeros(n)
  for i = 1:n
    factoredResponse[i]=find(z[i,:].==1)[1]
  end
#
#initialize Xbeta, update without doing full multiplications
  Xbeta=zeros(n,k)
#initialize penalty, rather than computing all the time
  groupPenalty=zeros(1)

#needed for gradient, can avoid computation
  z_mu=zeros(n,k)

#initialize gradient
  fullGradient=zeros(p*k)
  #get intercept
  intercept=reshape(ones(n),n,1)
  beta[1,:]=multinomialOptim(intercept,beta[1,:],z)
  beta[1,:]=beta[1,:]-mean(beta[1,:])
  Xbeta=Xbeta+intercept*beta[1,:]

  #update z_mu
  z_mu=z-multProbPrecompute(Xbeta)

  #majorized hessian
  fullHess=kron(1/2*eye(k),X'*X)/n

  #get eigenvalues for step from zero
  eigList=getEigs(X,categories,groups)
  invList=getInvs(X,categories,groups,k,n)
  #update gradient
  fullGradient=X'*z_mu
  #get first lambda and lambda sequence
  norms=zeros(sum(penaltyFactor.>0))
  count=0
  for i = 1:g
    if (penaltyFactor[i]>0)
      count=count+1
      cat=categories[i]
      inds=(groups.==cat)
      gradient=fullGradient[inds,:]/n
      norms[count]=(sum(gradient.^2)^0.5)
    end
  end
  lambdas=norms./((groupLengths[penaltyFactor.>0]*k).^0.5)
  lambdaStart=maximum(lambdas)
  num=100
  lambdaVec=lambdaStart*(collect(num:-1:(num-95))/num)
  #lambdaVec=userLambda
  betaList=zeros(p,k,length(lambdaVec))
  #the first set of coefficients only has the intercept
  betaList[:,:,1]=beta

  count=1
  if (length(lambdaVec)==1)
    return(betaList)
  end
  #lambdaVec=lambdaVec[1:30]

  #initialize active set
  activeCats=zeros(length(categories)).>1
  activeInds=zeros(size(X)[2]).>1

  for lambda = lambdaVec[2:length(lambdaVec)]
    count=count+1
    diff=tol+1

    #objective at start
    #adjust the group penalty for the new lambda
    groupPenalty=groupPenalty*lambda/lambdaVec[round(Int,(find(lambdaVec.==lambda))-1)][1]
    objEnd=0

#####related to active set
    #get active set
    lambda1=lambda
    lambda0=lambdaVec[round(Int,(find(lambdaVec.==lambda))-1)][1]
    activeCats=zeros(length(categories)).>1
    activeInds=zeros(size(X)[2]).>1
    getActiveCategories!(activeCats,activeInds,beta,Xbeta,categories,groups,X,z_mu,lambda,lambda0)

    #while kkt condition not satisfied
    converged=false
    while !converged
      while (diff>tol)
        betaStart=copy(beta)
        #update beta based on penalized quadratic approximation
        #modifies Xbeta, z_mu, and groupPenalty
        indsRep=repmat(activeInds,k)
        changedFlag=false
        beta[activeInds,:]=coordDescentQuadratic!(beta[activeInds,:],Xbeta,z_mu,groupPenalty,X[:,activeInds],z,groups[activeInds],categories[activeCats],lambda,penaltyFactor[activeCats],fullHess[indsRep,indsRep],eigList[activeCats],invList[activeCats])
        diff=maximum(diag(fullHess).*(vec(beta-betaStart).^2))
        #println(diff)
      end
      #objective=-loglikelihoodPrecompute(Xbeta,z)/n+groupPenalty
      #println(objective)
      #check kkt conditions
      converged=checkKKT(activeCats,activeInds,X,z_mu,groups,categories,lambda)
    end
    betaList[:,:,count]=beta
  end
  print(lambdaVec)
  return(betaList)
end

####################################
#optimize penalized quadratic group lasso problem
function coordDescentQuadratic!(beta,Xbeta,z_mu,groupPenalty,X,z,groups,categories,lambda,penaltyFactor,H,eigList,invList)
  tolMiddle=10.0^-8
  diffMiddle=tolMiddle+1
  quadraticObj=0
  totalObj=groupPenalty
  k=size(z_mu)[2]
  n=size(X)[1]
  fullGradient=-X'*z_mu/n
  betaDiffMiddle=zeros(size(beta))
  betaDiffTotal=zeros(size(beta))
  beta0=copy(beta)
  beta0_1=copy(beta)

  #H=-multinomialHessian(X,Xbeta)/n
  currentGrad=copy(fullGradient)
  #keep track of nonzero indices of beta, only use these for X*beta
  nonZeroInds=zeros(size(X)[2]).>1
  innerCount=1
  while diffMiddle>tolMiddle
    totalObj0=totalObj
    beta0_1[:]=copy(beta)
    nonZeroCats=zeros(length(categories))
    for i = 1:length(categories)
      cat=categories[i]

      inds=(groups.==cat)

      #starter quantities
      inds2=repmat(inds,k)
      #gradient=vec(fullGradient)[inds2]+H[inds2,:]*vec(betaDiffTotal)-H[inds2,inds2]*vec(beta)[inds2]
      gradient=vec(currentGrad[inds,:])-H[inds2,inds2]*vec(beta[inds,:])
      innerCount=innerCount+1
      hess=H[inds2,inds2]
      norm=sum(gradient.^2)^0.5

      #optimize group
      #get group related quantities so each group and penalty
      #can be updated one at a time
      penaltyStart=penalty(beta[inds,:],lambda*penaltyFactor[i])
      penaltyEnd=0
      betaStart=beta[inds,:]
      if norm>lambda*penaltyFactor[i]*sqrt(length(betaStart))
        nonZeroCats[i]=1
        #beta[inds,:]=innerOptimize(betaStart,lambda*penaltyFactor[i]*sqrt(length(betaStart)),gradient,hess,eigList[i])
        #beta[inds,:]=innerCyclicMM(betaStart,lambda*penaltyFactor[i]*sqrt(length(betaStart)),gradient,hess,eigList[i])
        #beta[inds,:]=innerFullMM(betaStart,lambda*penaltyFactor[i]*sqrt(length(betaStart)),gradient,hess,eigList[i],invList[i],k)
        beta[inds,:]=innerSoftThresholdMM(betaStart,lambda*penaltyFactor[i]*sqrt(length(betaStart)),gradient,hess,eigList[i])
        betaDiffMiddle[inds,:]=beta[inds,:]-beta0_1[inds,:]
        penaltyEnd=penalty(beta[inds,:],lambda*penaltyFactor[i])
        nonZeroInds[inds]=true
      end
      startEndZero=false
      if norm<=lambda*penaltyFactor[i]*sqrt(length(betaStart))
        if (all(beta[inds,:].==0))
          startEndZero=true
          betaDiffMiddle[inds,:]=0
        else
          beta[inds,:]=0
          nonZeroInds[inds]=false
          betaDiffMiddle[inds,:]=-beta0_1[inds,:]
          nonZeroCats[i]=0
        end
      end

      if (!startEndZero)
        currentGrad=currentGrad+reshape(H[:,inds2]*vec(betaDiffMiddle[inds,:]),size(currentGrad))
      end
      penaltyDiff=penaltyEnd-penaltyStart

      #update penalty
      groupPenalty[:]=groupPenalty+penaltyDiff
    end
    #@show sum(nonZeroCats)/length(categories)
    #@show sum(size(beta)[1])
    totalObj=quadraticObj+groupPenalty
    diffMiddle=maximum(diag(H).*(vec(betaDiffMiddle).^2))
  end
  Xbeta[:]=X[:,nonZeroInds]*beta[nonZeroInds,:]
  z_mu[:]=z-multProbPrecompute(Xbeta)
  return(beta)
end

##########################################
#group update function as implemented in another paper
function innerOptimize(betaK,lambda,gradient,hess,maxEig)
   tolInner=10.0^-8
   diffInner=tolInner+1
   if (all(betaK.==0))
     betaK=max(0,(1-lambda/(sum(gradient.^2)^0.5)))*(-1/maxEig*gradient)
   end

   while diffInner>tolInner
     startBetaK=copy(betaK)
     tempInds=zeros(length(betaK)).>1
     for i = 1:length(betaK)
       tempInds[i]=true
       #sum of other betas
       cons=sum(betaK[!tempInds].^2)
       grad=gradient[i]+hess[tempInds,!tempInds]*betaK[!tempInds]
       subHess=hess[tempInds,tempInds]

       grad=grad[1]
       rootSign=-abs(grad)/grad
       grad=-abs(grad)
       subHess=abs(subHess[1])
       cons=cons[1]
       function innerGradClosure(x)
         return grad+subHess*x+lambda*x/(sqrt(x^2+cons))
       end
       val=fzero(innerGradClosure,0,-grad/subHess)
       betaK[i]=rootSign*val
       tempInds[i]=false
     end
     endBetaK=betaK
     diffInner=maximum(abs(startBetaK-endBetaK))
   end
   return(betaK)
end

##########################################
#cyclic MM update within a group, using majorized loss
function innerCyclicMM(betaK,lambda,gradient,hess,maxEig)
  #avoid evaluation at zero
  if all(betaK.==0)
    betaK=max(0,(1-lambda/(sum(gradient.^2)^0.5)))*(-1/maxEig*gradient)
  end
  #vector of false's
  tempInds=zeros(length(betaK)).>1

  #precompute constant for penalty term
  consTot=sum(betaK.^2)
  for i = 1:length(betaK)
    tempInds[i]=true
    #sum of other beta
    grad=gradient[i]+hess[tempInds,!tempInds]*betaK[!tempInds]
    subHess=hess[tempInds,tempInds]
    grad=grad[1]
    subHess=subHess[1]
    val=betaK[i]
    val=val-(grad+subHess*val+lambda*val/sqrt(consTot))*(subHess+lambda*(consTot)^-0.5)^(-1)
    consTot=consTot-betaK[i]^2+val^2
    betaK[i]=val
    tempInds[i]=false
  end
  return(betaK)
end

##########################################
#MM whole-group update, using majorized penalty
function innerFullMM(betaK,lambda,gradient,hess,maxEig,eig_xx,k)
  #avoid evaluation at zero
  betaK=max(0,(1-lambda/(sum(gradient.^2)^0.5)))*(-1/maxEig*gradient)
  betaK=vec(betaK)
  consTot=sum(betaK.^2)
  gradTot=vec(gradient)+hess*betaK+lambda*vec(betaK)./consTot^0.5
  hessTot=hess+lambda*(1/consTot^0.5)*eye(size(hess)[1])
  test=kron(eye(k),eig_xx.vectors')'*inv(lambda/consTot^(0.5)*eye(size(hess)[1])+kron(1/2*eye(k),diagm(eig_xx.values)))*kron(eye(k),eig_xx.vectors')
  #betaK=betaK-hessTot\gradTot
  betaK=betaK-test*gradTot
  #@show maximum(abs(betaK betaK2))
  return(betaK)
end

############################################
#MM whole-group update using soft thresholding, using majorized loss
function innerSoftThresholdMM(betaK,lambda,gradient,hess,maxEig)
  if all(betaK.==0)
  end
  betaK=vec(betaK)
  betaR=-1/maxEig*(gradient+hess*betaK)+betaK
  betaK=max(0,(1-lambda/(maxEig*sum(betaR.^2)^0.5)))*(betaR)
  return(betaK)
end

#############################################
#helper function, computes l2 norm of a group and multiplies by lambda and sqrt(group length)
function penalty(beta,lambda)
  val=lambda*((length(beta))^0.5)*sum(beta.^2)^0.5
  return(val)
end

#############################################
#see which groups have a high chance of being non-zero, then we can do coordinate descent over these groups first
function getActiveCategories!(activeCats,activeInds,beta,Xbeta,categories,groups,X,z_mu,lambda1,lambda0)
  g=length(groups)
  fullGradient=-X'*z_mu/size(X)[1]
  for i=1:length(categories)
     cat=categories[i]
     inds=(groups.==cat)
     gradient=fullGradient[inds,:]
     norm=sum(gradient.^2)^0.5

     #for strong rule
     lambdaDiff=(2*lambda1-lambda0)*penaltyFactor[i]
     if (norm>(length(gradient))^(0.5)*lambdaDiff) | (penaltyFactor[i]==0)
       activeCats[i]=true
       activeInds[inds]=true
     else
       activeCats[i]=false
       activeInds[inds]=false
       if sum(beta[inds,:].!=0)>0
         Xbeta[:]=Xbeta-X[:,inds]*beta[inds,:]
         beta[inds,:]=0
       end
     end
   end
   z_mu[:]=z-multProbPrecompute(Xbeta)
end

######################################
#cycle through all inactive groups to see if any should be included in the model
function checkKKT(activeCats,activeInds,X,z_mu,groups,categories,lambda)
  fullGradient=-X'*z_mu/size(X)[1]
  converged=true
  for i = 1:length(categories)
    if (!activeCats[i])
      cat=categories[i]
      inds=(groups.==cat)
      gradient=fullGradient[inds,:]
      norm=sum(gradient.^2)^0.5
      if norm>length(gradient)^0.5*penaltyFactor[i]*lambda
        converged=false
        activeCats[i]=true
        activeInds[inds]=true
      end
    end
  end
  return(converged)
end

#################################
#get maximum eigenvalue of the hessian for each group for majorization step
function getEigs(X,categories,groups)
  eigs=zeros(length(categories))
  for i=1:length(categories)
    cat=categories[i]
    inds=(groups.==cat)
    eigs[i]=0.5*maximum(eigvals(X[:,inds]'*X[:,inds]))/size(X)[1]
  end
  return(eigs)
end

#################################
#precompute majorized hessian inverses
function getInvs(X,categories,groups,k,n)
  invList=Array(Any,0)
  for i=1:length(categories)
    cat=categories[i]
    inds=(groups.==cat)
    xx=X[:,inds]'*X[:,inds]/n
    push!(invList,eigfact(xx))
  end
  return(invList)
end


text:
# A-weighted sound pressure level
function f_aspl(data::PyObject, spl)

    f_AW = LinearInterpolation(data.aw_freq, data.aw_db)

    # Number of time steps
    n_t = size(spl)[1]

    aspl = zeros(eltype(spl), n_t)
    for i in 1:1:n_t
        weights = f_AW(data.f)
        aspl[i] = 10 * log10.( sum(10 .^((spl[i,:] .+ weights)./10.) ) )
    end

    return aspl
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
@testset "Soss" begin
    @testset "GP regression" begin
        k = SqExponentialKernel()
        y = randn(3)
        X = randn(3, 1)
        x = [rand(1) for _ in 1:3]

        gp_regression = Soss.@model X begin
            # Priors.
            α ~ LogNormal(0.0, 0.1)
            ρ ~ LogNormal(0.0, 1.0)
            σ² ~ LogNormal(0.0, 1.0)

            # Realized covariance function
            kernel = α * (SqExponentialKernel() ∘ ScaleTransform(1 / ρ))
            f = GP(kernel)

            # Sampling Distribution.
            y ~ f(X, σ² + 1e-9)
        end

        # Test for matrices
        m = gp_regression(; X=RowVecs(X))
        @test length(Soss.sample((m | (y=y,)), dynamichmc(), 5, 1)) == 5

        # Test for vectors of vector
        m = gp_regression(; X=x)
        @test length(Soss.sample((m | (y=y,)), dynamichmc(), 5, 1)) == 5
    end
    @testset "latent GP regression" begin
        X = randn(3, 1)
        x = [rand(1) for _ in 1:3]
        y = rand.(Poisson.(exp.(randn(3))))

        latent_gp_regression = Soss.@model X begin
            f = GP(Matern32Kernel())
            u ~ f(X)
            λ = exp.(u)
            y ~ For(eachindex(λ)) do i
                Poisson(λ[i])
            end
        end

        m = latent_gp_regression(; X=RowVecs(X))
        @test length(Soss.sample((m | (y=y,)), dynamichmc(), 5, 1)) == 5

        # Test for vectors of vector
        m = latent_gp_regression(; X=x)
        @test length(Soss.sample((m | (y=y,)), dynamichmc(), 5, 1)) == 5
    end
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using Test
using ExpectationStubs: DoNotCare

@testset "DoNotCare is equal to everything" begin
    @test DoNotCare{Int}()==4
    @test 5==DoNotCare{Integer}()
    @test 5==DoNotCare{Any}()
    @test !(5===DoNotCare{Any}())

    @test 5!=DoNotCare{String}()

    @test (Any, DoNotCare{Any}()) === (Any,DoNotCare{Any}())
    @test !((Any, 1) === (Any,DoNotCare{Int}()))
end


text:
module Patents

using Statistics
using StatsBase
using Dictionaries
using LightGraphs
using Requires
import Base: ==, hash

include("applications.jl")
include("classification.jl")
include("families.jl")
include("portfolios.jl")
include("citations.jl")
include("diversity.jl")
include("helpers.jl")
include("mainpath.jl")

# function __init__()
#     @require PlotlyJS="" include("plotting.jl")
# end

export Application, ApplicationID, Title, Abstract, Classification, NPLCitation

export
lensid, id, jurisdiction, status, docnr, kind, date,
title, abstract, lang, text,
applicants, inventors,
siblings,
classification, code, subgroup, maingroup, subclass, class, section,
cites, cites_npl, cites_count, cites_count_npl, citedby, citedby_count

export Family, applications, jurisdictions, dates, earliest_filing,
latest_filing, aggregate_families

export Portfolio, diversity, owner, families

export citationgraph, coclassification, normalized_citations

end


text:
using InteractionNetworkModels, Distributions

# On uniform data
# ---------------

S = [rand(1:4,rand(5:10)) for i in 1:7]

@time get_subseqs(S, 3)

data = [
    [rand(1:4,rand(5:30)) for i in 1:4] for i in 1:10
]

get_subseq_counts(data, 3)

x = get_subseq_proportions(data, 3)

dict_rank(x)

# On SIM model data
# -----------------

# The Model(s)

S = [[1,1], [1,2], [1,2,3], [1,4], [2,3]]
V = 1:13
d = MatchingDist(FastLCS(100))
model = SIM(S, 4.0, d, V, 50, 50)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=6, β=0.6, ν_td=1,
    len_dist=TrGeometric(0.8, 1, model.K_inner.u),
    lag=1,
    K=200)

@time out=mcmc_sampler(
    model,
    init=model.mode,
    desired_samples=100,
    burn_in=5000, lag=300
)

summaryplot(out)
out.sample
subseq_props = get_subseq_proportions(out.sample, 3)
dict_rank(subseq_props)

out.sample


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

using Markdown
using InteractiveUtils

# ╔═╡ f7e92500-194d-11eb-00c9-6f427ae012ac
md"# The Basel Problem

_Leonard Euler_ proved in 1741 that the following series

$$\frac{1}{1} + \frac{1}{4} + \frac{1}{9} + \cdots$$

converges to

$$\frac{\pi^2}{6} .$$
"

# ╔═╡ 287ee380-194e-11eb-01fd-517dabec4f05
n = 1 : 500000000

# ╔═╡ 31d135f0-194e-11eb-105b-51bb3ea5d3d8
sequence = n .^ -2

# ╔═╡ b6920b1e-194e-11eb-200d-17bf72fdff3a
π = sqrt(6 * sum(sequence))

# ╔═╡ Cell order:
# ╟─f7e92500-194d-11eb-00c9-6f427ae012ac
# ╠═b6920b1e-194e-11eb-200d-17bf72fdff3a
# ╠═31d135f0-194e-11eb-105b-51bb3ea5d3d8
# ╠═287ee380-194e-11eb-01fd-517dabec4f05

max_stars_repo_path: The Basel Problem.jl

max_stars_repo_name: tinfungster/Ways-to-get-PI

max_issues_repo_path: The Basel Problem.jl

max_issues_repo_name: tinfungster/Ways-to-get-PI

max_forks_repo_path: The Basel Problem.jl

max_forks_repo_name: tinfungster/Ways-to-get-PI


text:
using MeshCat: open_url
using MeshCat.ZMQServer: ZMQWebSocketBridge, web_url

function serve()
    zmq_url = nothing
    open_browser = false
    for arg in ARGS
        m = match(r"^\-\-zmq\-url=(.*)$", arg)
        if m !== nothing
            zmq_url = strip(m[1])
        end
        m = match(r"^\-\-open$", arg)
        if m !== nothing
            open_browser = true
        end
    end
    bridge = ZMQWebSocketBridge(zmq_url)
    if open_browser
        open_url(web_url(bridge))
    end
    run(bridge)
end

serve()


text:
export AddAttention

doc"""
    AddAttention(x1, x2)

```math
f(x_{i},x_{j}) = W^{T} \sigma (W_{1}x_{i} + W_{2}x_{j} + b) + b
```
"""
struct AddAttention
    l1
    l2
end

function AddAttention(::Type{T}, insize::Int, outsize::Int) where T
    l1 = Linear(T, 2insize, insize)
    l2 = Linear(T, insize, outsize)
    AddAttention(l1, l2)
end

function (f::AddAttention)(x1::Var, x2::Var)
    h = pairwise(x1, x2)
    h = f.l1(h)
    h = relu(h)
    f.l2(h)
end


text:
using DynamicGrids: init

"""
Optimisation mode for running replicate simulations inside an optimiser.
"""
abstract type Replicates end
"""
Use Thread.@threads to run replicate simulations
"""
struct ThreadedReplicates <: Replicates end
"""
Use @distributed to run replicate simulations
"""
struct DistributedReplicates <: Replicates end
"""
Use a single processor to run replicate simulations
"""
struct SingleCoreReplicates <: Replicates end

"""
    Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize,
                 tspan, threading)

Parametrizer functor to use with Optim.jl or similar.

# Arguments

- `ruleset::Ruleset`: simulation ruleset, with `init` array attached
- `output::Output`: optional output type. By default an ArrayOutput will be generated.
- `objective::Objective`: objective data
- `transform`: single argument function to transform targets and predictions before the loss function
- `loss`: LossFunctions.jl loss function
- `ngroups::Int`: number of replicate simulation
- `groupsize::Int`: number of simulations in a group. Larger groups may inprove distributed performance.
- `threading::Replicates`: Type to define threading mode for simulation replicates:
  `SingleCoreReplicates`, `ThreadedReplicates`, `DistributedReplicates`
"""
struct Parametriser{RU,OP,OB,TR,L,NR,GS,TH,D,TB,PB,RE}
    ruleset::RU
    output::OP
    objective::OB
    transform::TR
    loss::L
    ngroups::NR
    groupsize::GS
    threading::TH
    data::D
    targetbuffer::TB
    predictionbuffer::PB
    results::RE
end

DynamicGrids.ruleset(p::Parametriser) = p.ruleset

data(p::Parametriser) = p.data
output(p::Parametriser) = p.output
transform(p::Parametriser) = p.transform
objective(p::Parametriser) = p.objective
loss(p::Parametriser) = p.loss
ngroups(p::Parametriser) = p.ngroups
groupsize(p::Parametriser) = p.groupsize
threading(p::Parametriser) = p.threading
targetbuffer(p::Parametriser) = p.targetbuffer
predictionbuffer(p::Parametriser) = p.predictionbuffer
results(p::Parametriser) = p.results
DynamicGrids.tspan(p::Parametriser) = tspan(output(p))

Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize,
             threading::ThreadedReplicates) = begin
    targetbuffer = transform.(targets(objective))
    # Make copies of anything threads will write to
    predictionbuffer = [transform.(targets(objective)) for i in 1:Threads.nthreads()]
    outputs = [deepcopy(output) for i in 1:Threads.nthreads()]
    data = [DynamicGrids.SimData(deepcopy(extent(output)), deepcopy(ruleset)) for i in 1:Threads.nthreads()]
    results = [zeros(groupsize) for g in 1:ngroups]

    Parametriser(ruleset, outputs, objective, transform, loss, ngroups, groupsize,
                 threading, data, targetbuffer, predictionbuffer, results)
end

Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize,
             threading=SingleCoreReplicates()) = begin
    targetbuffer = transform.(targets(objective))
    predictionbuffer = transform.(targets(objective))
    data = DynamicGrids.SimData(deepcopy(extent(output)), ruleset)
    results = [zeros(groupsize) for g in 1:ngroups]

    Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize,
                 threading, data, targetbuffer, predictionbuffer, results)
end

"""
    (p::Parametriser)(params)

Provides an objective function for an optimiser like Optim.jl
"""
(p::Parametriser)(params) = begin
    # Rebuild the rules with the current parameters
    p.ruleset.rules = Flatten.reconstruct(p.ruleset.rules, params, Real)
    obj = objective(p)

    names = fieldnameflatten(rules(ruleset(p)), Real)
    println("\nParameters: \n", typeof(rules(ruleset(p))))
    display(collect(zip(names, params)))
    println()

    replicate!(threading(p), p, params)

    meanloss = sum(sum.(results(p))) / (ngroups(p) * groupsize(p))
    println("mean loss from $(p.ngroups * p.groupsize): ", meanloss, "\n")
    return meanloss
end

replicate!(::ThreadedReplicates, p::Parametriser, params) = begin
    obj = p.objective
    Threads.@threads for g in 1:p.ngroups
        id = Threads.threadid()
        o = output(p)[id]
        predbuffer = predictionbuffer(p)[id]
        for n in 1:p.groupsize
            sim!(o, ruleset(p); simdata=data(p)[id])
            predbuffer .= transform(p).(predictions(obj, o))
            results(p)[g][n] = value(loss(p), targetbuffer(p), predbuffer, AggMode.Sum())
        end
    end
end

# replicate!(::DistributedReplicates, p, params) = begin
#     obj = p.objective
#     cumsum = @distributed (+) for g in 1:p.ngroups
#         for n in 1:p.groupsize
#             sim!(p.output, p.ruleset; tstop=p.tstop)
#             p.predictionbuffer .= p.transform.(predictions(p.objective, p.output))
#             p.results[g][n] = value(p.loss, p.targetbuffer, p.predictionbuffer, AggMode.Sum())
#         end
#         # Core.println("group: ", g, " - reps ($(p.groupsize)) mean loss: ", grouploss/p.groupsize)
#     end
#     cumsum
# end

replicate!(::SingleCoreReplicates, p::Parametriser, params) = begin
    obj = p.objective
    for g in 1:p.ngroups
        grouploss = 0.0
        for n in 1:p.groupsize
            println(n)
            sim!(output(p), ruleset(p); simdata=data(p))
            p.predictionbuffer .= p.transform.(predictions(p.objective, p.output))
            loss = value(p.loss, p.targetbuffer, p.predictionbuffer, AggMode.Sum())
            grouploss += loss
            p.results[g][n] = loss
        end
        println("group: ", g, " - reps ($(p.groupsize)) mean loss: ", grouploss / p.groupsize)
    end
end


text:
#### Frank wolfe for matrix completion
using FirstOrderOptimization
using PyPlot
import Base: axpy!, scale!

## Generate data

# make it deterministic
srand(10)

m, n, r, rtrue = 100, 100, 10, 2
X = randn(m, rtrue) * randn(rtrue, n)
smooth_inf = log(min(m,n))
# smooth_inf and smooth_1 are conjugate
# 1/p + 1/q = 1 => p = 1/(1-1/q) = q / (q - 1)
smooth_1 = smooth_inf / (smooth_inf - 1)
u,s,v = svd(X)
tau = sum(s)
tau_smooth = sum(s.^(smooth_1))^(1/smooth_1)

pobs = .5
obs = sprand(m, n, pobs).>0
iobs = find(obs)
A = IndexingOperator(m, n, iobs)
b = A*X

srand(10)

objective_sketched(z) = .5*sum((z - b).^2)
grad_objective_sketched(z) = LowRankOperator(A, z-b, transpose = Symbol[:T, :N]) # concretely, should be A'*(z-b)

const_nucnorm(z) = tau # we'll always saturate the constraint, don't bother computing it
function min_lin_st_nucnorm_sketched(g, tau)
	u,s,v = svds(Array(g), nsv=1) # for this case, g is a sparse matrix so representing it is O(m)
	return LowRankOperator(-tau*u, v')
end
# I can't think of a pretty way to compute
# <g, tilde_x - x> = <A'(z-b), u tau v' - x> = <z-b, A(u tau v') - A(x)> = <z-b, A(u tau v') - z>
# dot_g_w_tx_minus_x(g, tilde_x, x) = dot(g.factors[2], g.factors[1]*tilde_x - x)
dot_g_w_tx_minus_x(g, Delta, z) = dot(vec(g.factors[2]), g.factors[1]*Delta - z)
# we end up computing A*Delta twice with this scheme
update_var!(z, Delta, a) = (scale!(z, 1-a); axpy!(a, A*Delta, z))

z = zeros(length(b))
ch = ConvergenceHistory("fw_sketched")
X_sketched = frank_wolfe_sketched(
	z,
	objective_sketched, grad_objective_sketched,
	const_nucnorm,
	tau,
	min_lin_st_nucnorm_sketched,
	AsymmetricSketch(m,n,r),
	FrankWolfeParams(100, 1e-2, DecreasingStepSize(2,1)),
	ch,
	dot_g_w_tx_minus_x = dot_g_w_tx_minus_x,
	update_var! = update_var!,
	#LB = 0,
	verbose=true
	)

figure()
# yscale("log")
plot(ch.objective, label="fw")
# plot(ch.objective, label="fw")
xlabel("Iterations")
ylabel("Objective")

## now smooth version

for (nsv, nsv_label) = [(r, "r"), (Int(ceil(smooth_inf)), "log(n)"), (min(m,n), "n")]

	function min_lin_st_nucnorm_sketched(g, tau,
			smooth_inf = smooth_inf,
			nsv = nsv)
		u,s,v = svds(Array(g), nsv=nsv) # for this case, g is a sparse matrix so representing it is O(m)
		s = s.^smooth_inf
		s = s / sum(s.^smooth_1).^(1/smooth_1) * -tau # normalize so the smooth_1 norm of the solution is tau
		return LowRankOperator(u, spdiagm(s), v')
	end
	const_nucnorm(z) = tau_smooth # we'll always saturate the constraint, don't bother computing it

	srand(10)

	ch_smooth = ConvergenceHistory("smooth_fw_sketched")
	z = zeros(length(b))
	X_sketched_smooth = frank_wolfe_sketched(
		z,
		objective_sketched, grad_objective_sketched,
		const_nucnorm,
		tau_smooth,
		min_lin_st_nucnorm_sketched,
		AsymmetricSketch(m,n,r),
		FrankWolfeParams(100, 1e-2, DecreasingStepSize(2,1)),
		ch_smooth,
		dot_g_w_tx_minus_x = dot_g_w_tx_minus_x,
		update_var! = update_var!,
		#LB = 0,
		verbose=true
		)

	plot(ch_smooth.objective, label="smoothed fw top "*nsv_label)
	#@show(vecnorm(Array(X_sketched_smooth) - X))
end
legend()

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
export CSRMatrix, insertGlobalValues

"""
An implementation of [`RowMatrix`](@ref) that uses CSR format
"""
mutable struct CSRMatrix{Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer} <: RowMatrix{Data, GID, PID, LID}
    rowMap::BlockMap{GID, PID, LID}
    colMap::Union{BlockMap{GID, PID, LID}, Nothing}

    importMV::Union{MultiVector{Data, GID, PID, LID}, Nothing}
    exportMV::Union{MultiVector{Data, GID, PID, LID}, Nothing}

    myGraph::CSRGraph{GID, PID, LID}

    localMatrix::LocalCSRMatrix{Data, LID}

    values2D::Array{Array{Data, 1}, 1}

    #pull storageStatus and fillComplete  from graph

    #Dict keys are row indices
    #first element of each tuple is a column index
    #second element of each tuple is the matching entry
    nonlocals::Dict{GID, Tuple{Array{Data, 1}, Array{GID, 1}}}

    plist::Dict

    function CSRMatrix{Data, GID, PID, LID}(rowMap::BlockMap{GID, PID, LID}, colMap::Union{BlockMap{GID, PID, LID}, Nothing}, myGraph::CSRGraph{GID, PID, LID}, localMatrix::LocalCSRMatrix{Data, LID}, plist::Dict) where {Data, GID, PID, LID}

        #allocate values
        localNumRows = getLocalNumRows(myGraph)
        if getProfileType(myGraph) == STATIC_PROFILE
            ptrs = myGraph.rowOffsets
            localTotalNumEntries = ptrs[localNumRows+1]

            resize!(localMatrix.values, localTotalNumEntries)

            values2D = Array{Array{Data, 1}, 1}(undef, 0)
        else #DYNAMIC_PROFILE
            if isLocallyIndexed(myGraph)
                graphIndices = myGraph.localIndices2D
            else
                graphIndices = myGraph.globalIndices2D
            end
            values2D = Array{Array{Data, 1}, 1}(undef, localNumRows)
            for r = 1:length(graphIndices)
                values2D[r] = Array{Data, 1}(undef, length(graphIndices[r]))
            end
        end

        new(rowMap,
            colMap,
            nothing,
            nothing,
            myGraph,
            localMatrix,
            values2D,
            Dict{GID, Tuple{Array{Data, 1}, Array{GID, 1}}}(),
            plist)
    end

end

#### Constructors ####
#TODO document Constructors

function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType; plist...) where {Data, GID, PID, LID}
    CSRMatrix{Data}(rowMap, maxNumEntriesPerRow, pftype, Dict(plist))
end
function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType, plist::Dict) where {Data, GID, PID, LID}
    CSRMatrix{Data}(rowMap, nothing,
        maxNumEntriesPerRow, pftype, plist)
end

#=  These should be covered by the following constructor
function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        colMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType; plist...) where {Data, GID, PID, LID}
    CSRMatrix{Data}(rowMap, colMap, maxNumEntriesPerRow, pftype, Dict(plist))
end
function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        colMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType, plist::Dict) where {Data, GID, PID, LID}
    CSRMatrix{Data}(rowMap, colMap, maxNumEntriesPerRow,
        pftype, plist)
end
=#

function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType; plist...) where {Data, GID, PID, LID}
    CSRMatrix{Data}(rowMap, colMap, maxNumEntriesPerRow, pftype, Dict(plist))
end
function CSRMatrix{Data}(rowMap::BlockMap{GID, PID, LID},
        colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        maxNumEntriesPerRow::Union{Integer, Array{<:Integer, 1}},
        pftype::ProfileType, plist::Dict) where {Data, GID, PID, LID}
    graph = CSRGraph(rowMap, maxNumEntriesPerRow, pftype, plist)

    matrix = CSRMatrix{Data, GID, PID, LID}(rowMap, colMap,
        graph, LocalCSRMatrix{Data, LID}(), plist)

    resumeFill(matrix, plist)

    matrix
end

function CSRMatrix{Data}(graph::CSRGraph{GID, PID, LID}; plist...
        ) where {Data, GID, PID, LID}
    CSRMatrix{Data}(graph, Dict(plist))
end
function CSRMatrix{Data}(graph::CSRGraph{GID, PID, LID},plist::Dict
        ) where {Data, GID, PID, LID}
    numCols = numMyElements(getColMap(graph))
    localGraph = getLocalGraph(graph)
    val = Array{Data, 1}(undef, length(localGraph.entries))
    localMatrix = LocalCSRMatrix(numCols, val, localGraph)

    CSRMatrix(graph.rowMap, graph.colMap, graph, localMatrix, plist)
end

function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        rowOffsets::AbstractArray{LID, 1}, colIndices::AbstractArray{LID, 1}, values::AbstractArray{Data, 1};
        plist...) where {Data, GID, PID, LID}
    CSRMatrix(rowMap, colMap, rowOffsets, colIndices, values, Dict(plist))
end
function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        rowOffsets::AbstractArray{LID, 1}, colIndices::AbstractArray{LID, 1}, values::AbstractArray{Data, 1},
        plist::Dict) where {Data, GID, PID, LID}

    #check user's input.  Might throw on only some processes, causing deadlock
    if length(values) != length(colIndices)
        throw(InvalidArgumentError("values and columnIndices must "
                * "have the same length"))
    end

    graph = CSRGraph(rowMap, colMap, rowOffsets, columnIndices, plist)
    localGraph = getLocalGraph(graph)

    numCols = numMyElements(getColMap(graph))
    localMatrix = LocalCSRMatrix(numCols, values, localGraph)

    CSRMatrix(rowMap, colMap, graph, localMatrix, plist)
end

function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localMatrix::LocalCSRMatrix{Data, LID}; plist...
        ) where {Data, GID, PID, LID}
    CSRMatrix(rowMap, colMap, localMatrix, Dict(plist))
end
function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localMatrix::LocalCSRMatrix{Data, LID}, plist::Dict
        ) where {Data, GID, PID, LID}

    graph = CSRGraph(rowMap, colMap, localMatrix.graph, plist)

    matrix = CSRMatrix(rowMap, colMap, graph, localMatrix, plist)

    computeGlobalConstants(matrix)
    matrix
end

function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localMatrix::AbstractArray{Data, 2}; plist...
        ) where {Data, GID, PID, LID}
    CSRMatrix(rowmap, colMap, localMatrix, Dict(plist))
end
function CSRMatrix(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localMatrix::AbstractArray{Data, 2}, plist::Dict
        ) where {Data, GID, PID, LID}
    linearIndices = find(x -> x!=0, localMatrix)
    rowIndicesIter, colIndicesIter, valuesIter = zip(
        sort!(collect(zip(ind2sub(size(localMatrix), linearIndices)...,
                          localMatrix[linearIndices])))...)
    rowIndices = collect(rowIndicesIter)
    rowOffsets = Array{LID, 1}(undef, size(localMatrix, 1)+1)
    row = 1
    j = 1
    for i in LID(1):LID(length(rowIndices))
        if rowIndices[i] > row
            row += 1
            rowOffsets[row] = i
        end
    end
    rowOffsets[length(rowOffsets)] = length(rowIndices)+1

    CSRMatrix(rowMap, colMap, rowOffsets,
        collect(colIndicesIter), collect(valuesIter), plist)
end

#### Internal methods ####
function combineGlobalValues(matrix::CSRMatrix{Data, GID, PID, LID},
        globalRow::GID, indices::AbstractArray{GID, 1},
        values::AbstractArray{Data, 1}, cm::CombineMode
        ) where {Data, GID, PID, LID}

    if cm == ADD || cm == INSERT
        insertGlobalValuesFiltered(globalRow, indices, values)
    else
        #TODO implement ABSMAX and REPLACE
        #not implmenented in TPetra, because its not a common use case and difficult (see FIXME on line 6225)
        throw(InvalidArgumentError("Not yet implemented for combine mode $cm"))
    end
end

#does nothing, exists only to be a parallel to CSRGraph
computeGlobalConstants(matrix::CSRMatrix) = nothing

#Tpetra's only clears forbNorm, exists only to be a parallel to CSRGraph
clearGlobalConstants(matrix::CSRMatrix) = nothing

function globalAssemble(matrix::CSRMatrix)
    comm = getComm(matrix)
    if !isFillActive(matrix)
        throw(InvalidStateError("Fill must be active to call this method"))
    end

    myNumNonLocalRows = length(matrix.nonlocals)
    nooneHasNonLocalRows = maxAll(comm, myNumNonLocalRows == 0)
    if nooneHasNonLocalRows
        #no process has nonlocal rows, so nothing to do
        return
    end

    #nonlocalRowMap = BlockMap{GID, PID, LID}()

    numEntPerNonlocalRow = Array{Integer, 1}(undef, numNonLocalRows)
    myNonlocalGlobalRows = Array{GID, 1}(undef, numNonLocalRows)

    curPos = 1
    for (key, val) in matrix.nonlocal
        myNonlocalGlobalRows[curPos] = key

        values = val[1]
        globalColumns = val[2]

        order = sortperm(globalColumns)
        permute!(globalColumns, order)
        permute!(values, order)

        curPos += 1
    end

    #merge2
    if length(globalColumns) > 0
        setIndex = 1
        valuesSum = 0
        currentIndex = globalColumns[1]
        for i = 1:length(globalColumns)
            if currentIndex != globalColumns[i]
                values[setIndex] = valuesSum
                globalColumns[setIndex] = currentIndex
                setIndex += 1
                valuesSum = 0
            end
            valuesSum += values[i]
            i += 1
        end
        values[setIndex] = valuesSum
        globalColumns[setIndex] = currentIndex

        resize!(values, setIndex)
        resize!(globalColumns, setIndex)
    end

    numEntPerNonloalRow[curPose] = length(globalColumns)

    #don't need to worry about finding the indexBase
    nonlocalRowMap = BlockMap(0, myNonlocalGlobalRows, comm)

    nonlocalMatrix = CSRMatrix(nonlocalRowMap, numEntPernonlocalRow, STATIC_PROFILE)

    curPos = 1
    for (key, val) in matrix.nonlocals
        globalRow = key

        vals = val[1]
        globalCols = val[2]

        insertGlobalValues(nonlocalMatrix, globalRow, globalCols, vals)
    end

    origRowMap = rowMap(matrix)
    origRowMapIsOneToOne = isOneToOne(origRowMap)

    if origRowMapIsOneToOne
        exportToOrig = Export(nonlocalRowMap, origRowMap)
        isLocallyComplete = isLocallyComplete(exportToOrig)
        doExport(nonlocalMatrix, matrix, exportToOrig, ADD)
    else
        oneToOneRowMap = createOneToOne(origRowMap)
        exportToOneToOne = Export(nonlocalRowMap, oneToOneRowMap)

        isLocallyComplete = isLocallyComplete(exportToOneToOne)

        oneToOneMatrix = CSRMatrix{Data}(oneToOneRowMap, 0)

        doExport(nonlocalMatrix, onToOneMatrix, exportToOneToOne, ADD)

        #nonlocalMatrix = null

        importToOrig = Import(oneToOneRowMap, origRowMap)
        doImport(oneToOneMatrix, matrix, importToOrig, ADD)
    end

    empty!(matrix.nonlocals)

    globallyComplete = minAll(comm, isLocallyComplete)
    if !globallyComplete
        throw(InvalidArgumentError("On at least one process, insertGlobalValues "
                * "was called with a global row index which is not in the matrix's "
                * "row map on any process in its communicator."))
    end
end

function fillLocalGraphAndMatrix(matrix::CSRMatrix{Data, GID, PID, LID},
        plist::Dict) where {Data, GID, PID, LID}
    localNumRows = getLocalNumRows(matrix)

    myGraph = matrix.myGraph
    localMatrix = matrix.localMatrix

    matrix.localMatrix.graph.entries = myGraph.localIndices1D

    #most of the debug sections were taken out
    if getProfileType(matrix) == DYNAMIC_PROFILE
        numRowEntries = myGraph.numRowEntries

        ptrs = Array{LID, 1}(undef, localNumRows+1)
        localTotalNumEntries = computeOffsets(ptrs, numRowEntries)

        inds = Array{LID, 1}(undef, localTotalNumEntries)
        #work around type instability required by localMatrix.values
        vals_concrete = Array{Data, 1}(undef, localTotalNumEntries)
        vals = vals_concrete

        localIndices2D = myGraph.localIndices2D
        for row = 1:localNumRows
            numEntries = numRowEntries[row]
            dest = range(ptrs[row], step=1, length=numEntries)

            inds[dest] = localIndices2D[row][1:numEntries]
            vals_concrete[dest] = matrix.values2D[row][1:numEntries]
        end
    elseif getProfileType(matrix) == STATIC_PROFILE
        curRowOffsets = myGraph.rowOffsets

        if myGraph.storageStatus == STORAGE_1D_UNPACKED
            #pack row offsets into ptrs

            localTotalNumEntries = 0

            ptrs = Array{LID, 1}(undef, localNumRows + 1)
            numRowEntries = myGraph.numRowEntries
            localTotalNumEntries = computeOffsets(ptrs, numRowEntries)

            inds = Array{LID, 1}(undef, localTotalNumEntries)
            #work around type instability required by localMatrix.values
            vals_concrete = Array{Data, 1}(undef, localTotalNumEntries)
            vals = vals_concrete

            #line 1234
            for row in 1:localNumRows
                srcPos = curRowOffsets[row]
                dstPos = ptrs[row]
                dstEnd = ptrs[row+1]-1
                dst = dstPos:dstEnd
                src = srcPos:srcPos+dstEnd-dstPos

                inds[dst] = myGraph.localIndices1D[src]
                vals_concrete[dst] = localMatrix.values[src]
            end
        else
            #dont have to pack, just set pointers
            ptrs = myGraph.rowOffsets
            inds = myGraph.localIndices1D
            vals = localMatrix.values
        end
    end

    if get(plist, :optimizeStorage, true)
        empty!(myGraph.localIndices2D)
        empty!(myGraph.numRowEntries)

        empty!(matrix.values2D)

        myGraph.rowOffsets = ptrs
        myGraph.localIndices1D = inds

        myGraph.pftype = STATIC_PROFILE
        myGraph.storageStatus = STORAGE_1D_PACKED
    end

    myGraph.localGraph = LocalCSRGraph(inds, ptrs)
    matrix.localMatrix = LocalCSRMatrix(myGraph.localGraph, vals, getLocalNumCols(matrix))
end

function insertNonownedGlobalValues(matrix::CSRMatrix{Data, GID, PID, LID},
        globalRow::GID, indices::AbstractArray{GID, 1}, values::AbstractArray{Data, 1}
        ) where {Data, GID, PID, LID}

    curRow = matrix.nonlocals[globalRow]
    curRowVals = curRow[1]
    curRowInds = curRow[2]

    newCapacity = length(curRowInds) + length(indices)

    append!(curRowVals, values)
    append!(curRowInds, indices)
end

function getView(matrix::CSRMatrix{Data, GID, PID, LID}, rowInfo::RowInfo{LID})::AbstractArray{Data, 1} where {Data, GID, PID, LID}
    if getProfileType(matrix) == STATIC_PROFILE && rowInfo.allocSize > 0
        range = rowInfo.offset1D:rowInfo.offset1D+rowInfo.allocSize-LID(1)
        baseArray = matrix.localMatrix.values
        if baseArray isa Vector{Data}
            view(matrix.localMatrix.values, range)
        else
            view(matrix.localMatrix.values, range)
        end
    elseif getProfileType(matrix) == DYNAMIC_PROFILE
        baseArray = matrix.values2D[rowInfo.localRow]
        view(baseArray, LID(1):LID(length(baseArray)))
    else
        Data[]
    end
end

function getDiagCopyWithoutOffsets(rowMap, colMap, A::CSRMatrix{Data}) where {Data}
    errCount = 0

    D = Array{Data, 1}(undef, getNumLocalRows(A))

    for localRowIndex = 1:length(D)
        D[localRowIndex] = 0
        globalIndex = gid(rowMap, localRowIndex)
        localColIndex = lid(colMap, globalIndex)
        if localColIndex != 0
            colInds, vals = getLocalRowView(A, localRowIndex)

            offset = 1
            numEnt = length(curRow)
            while offset <= numEnt
                if colInds[offset] == localColIndex
                    break;
                end
                offset += 1
            end

            if offset > numEnt
                errCount += 1
            else
                D[localRowIndex] = vals[offset]
            end
        end
    end
    D
end

function sortAndMergeIndicesAndValues(matrix::CSRMatrix{Data, GID, PID, LID},
        sorted, merged) where {Data, GID, PID, LID}
    graph = getGraph(matrix)
    localNumRows = getLocalNumRows(graph)
    totalNumDups = 0

    for localRow in LID(1):localNumRows
        rowInfo = getRowInfo(graph, localRow)
        if !sorted
            inds, vals = getLocalRowView(matrix, rowInfo)

            order = sortperm(inds)
            permute!(inds, order)
            permute!(vals, order)
        end
        if !merged
            totalNumDups += mergeRowIndicesAndValues(matrix, rowInfo)
        end

        recycleRowInfo(rowInfo)
    end

    if !sorted
        graph.indicesAreSorted = true
    end
    if !merged
        graph.nodeNumEntries -= totalNumDups
        graph.noRedundancies = true
    end
end

function mergeRowIndicesAndValues(matrix::CSRMatrix{Data, GID, PID, LID},
        rowInfo::RowInfo{LID})::LID where {Data, GID, PID, LID}

    graph = getGraph(matrix)
    indsView, valsView = getLocalRowView(matrix, rowInfo)

    if rowInfo.numEntries != 0
        newend = 1
        for cur in 2:rowInfo.numEntries
            if indsView[newend]::LID != indsView[cur]::LID
                #new entry, save it
                newend += 1
                indsView[newend] = indsView[cur]::LID
                valsView[newend] = valsView[cur]::Data
            else
                #old entry, merge it
                valsView[newend] += valsView[cur]::Data
            end
        end
    else
        newend = 0
    end

    graph.numRowEntries[rowInfo.localRow] = newend

    rowInfo.numEntries - newend
end

#### External methods ####
#TODO document external methods

function insertGlobalValues(matrix::CSRMatrix{Data, GID, PID, LID}, globalRow::Integer,
        indices::AbstractArray{GID, 1}, values::AbstractArray{Data, 1}
        ) where {Data, GID, PID, LID}
    myGraph = matrix.myGraph

    localRow = lid(getRowMap(matrix), globalRow)

    if localRow == 0
        insertNonownedGlobalValues(matrix, globalRow, indices, values)
    else
        numEntriesToInsert = length(indices)
        if hasColMap(matrix)
            colMap = getColMap(matrix)

            for k = 1:numEntriesToInsert
                if !myGID(colMap, indices[k])
                    throw(InvalidArgumentError("Attempted to insert entries into "
                            * "owned row $globalRow, at the following column indices "
                            * "$indices.  At least one of those indices ($(indices[k])"
                            * ") is not in the column map on this process"))
                end
            end
        end

        rowInfo = getRowInfo(myGraph, localRow)
        curNumEntries = rowInfo.numEntries
        newNumEntries = curNumEntries + numEntriesToInsert
        if newNumEntries > rowInfo.allocSize
            if(getProfileType(matrix) == STATIC_PROFILE
                    && newNumEntries > rowInfo.allocSize)
                throw(InvalidArgumentError("number of new indices exceed "
                        * "statically allocated graph structure"))
            end

            updateGlobalAllocAndValues(myGraph, rowInfo, newNumEntries,
                            matrix.values2D[localRow])

            recycleRowInfo(rowInfo)
            rowInfo = getRowInfo(myGraph, localRow);

        end

        insertIndicesAndValues(myGraph, rowInfo, indices, getView(matrix, rowInfo),
                values, GLOBAL_INDICES)

        recycleRowInfo(rowInfo)
    end
    nothing
end

function resumeFill(matrix::CSRMatrix, plist::Dict)
    resumeFill(matrix.myGraph, plist)

    clearGlobalConstants(matrix)
    #graph handles fillComplete variable
end

fillComplete(matrix::CSRMatrix; plist...) = fillComplete(matrix, Dict(plist))

function fillComplete(matrix::CSRMatrix, plist::Dict)
    #TODO figure out if the second arg should be getColMap(matrix)
    fillComplete(matrix, getRowMap(matrix), getRowMap(matrix), plist)
end

function fillComplete(matrix::CSRMatrix{Data, GID, PID, LID},
        domainMap::BlockMap{GID, PID, LID}, rangeMap::BlockMap{GID, PID, LID};
        plist...) where {Data, GID, PID, LID}
    fillComplete(matrix, domainMap, rangeMap, Dict(plist))
end

function fillComplete(matrix::CSRMatrix{Data, GID, PID, LID},
        domainMap::BlockMap{GID, PID, LID}, rangeMap::BlockMap{GID, PID, LID},
        plist::Dict) where {Data, GID, PID, LID}
    if isFillComplete(matrix)
        throw(InvalidStateError(
                "Matrix cannot be fill complete when fillComplete(...) is called"))
    end

    myGraph = matrix.myGraph

    assertNoNonlocalInserts = get(plist, :noNonlocalChanges, false)
    #skipping sort ghosts stuff

    numProcs = numProc(getComm(matrix))

    needGlobalAssemble = !assertNoNonlocalInserts && numProcs > 1
    if needGlobalAssemble
        globalAssemble(matrix)
    else
        if numProcs == 1 && length(matrix.nonlocals) != 0
            throw(InvalidStateError("Cannot have nonlocal entries on a serial run.  An invalid entry is present."))
        end
    end

    setDomainRangeMaps(myGraph, domainMap, rangeMap)
    if !hasColMap(myGraph)
        makeColMap(myGraph)
        matrix.colMap = myGraph.colMap
    end

    makeIndicesLocal(myGraph)

    sortAndMergeIndicesAndValues(matrix, isSorted(myGraph), isMerged(myGraph))

    makeImportExport(myGraph)
    computeGlobalConstants(myGraph)

    myGraph.fillComplete = true
    checkInternalState(myGraph)

    fillLocalGraphAndMatrix(matrix, plist)
end

getProfileType(mat::CSRMatrix) = getProfileType(mat.myGraph)
isStorageOptimized(mat::CSRMatrix) = isStorageOptimized(mat.myGraph)

function getLocalDiagOffsets(matrix::CSRMatrix{Data, GID, PID, LID})::AbstractArray{LID, 1} where {Data, GID, PID, LID}
    graph = matrix.myGraph
    localNumRows = getLocalNumRows(graph)
    getLocalDiagOffsets(graph)
end

#### Row Matrix functions ####

getRowMap(mat::CSRMatrix) = mat.rowMap
hasColMap(mat::CSRMatrix) = mat.colMap != nothing
getColMap(mat::CSRMatrix) = mat.colMap
getGraph(mat::CSRMatrix) = mat.myGraph

function setColumnMapMultiVector(mat::CSRMatrix{Data, GID, PID, LID}, mv::Union{MultiVector{Data, GID, PID, LID}, Nothing}) where {Data, GID, PID, LID}
    mat.importMV = mv
end

function setRowMapMultiVector(mat::CSRMatrix{Data, GID, PID, LID}, mv::Union{MultiVector{Data, GID, PID, LID}, Nothing}) where {Data, GID, PID, LID}
    mat.exportMV = mv
end

getColumnMapMultiVector(mat::CSRMatrix) = mat.importMV
getRowMapMultiVector(mat::CSRMatrix) = mat.exportMV

function getGlobalRowCopy(matrix::CSRMatrix{Data, GID, PID, LID},
        globalRow::Integer
        )::Tuple{Array{GID, 1}, Array{Data, 1}} where {Data, GID, PID, LID}
    myGraph = matrix.myGraph

    rowInfo = getRowInfoFromGlobalRow(myGraph, GID(globalRow))
    viewRange = 1:rowInfo.numEntries

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return (GID[], Data[])
    end

    @inbounds if isLocallyIndexed(myGraph)
        colMap = getColMap(myGraph)
        curLocalIndices = getLocalView(myGraph, rowInfo)[viewRange]
        curGlobalIndices = @. gid(colMap, curLocalIndices)
    else
        curGlobalIndices = getGlobalView(myGraph, rowInfo)[viewRange]
    end
    @inbounds curValues = getView(matrix, rowInfo)[viewRange]

    retVal = (curGlobalIndices, curValues)
    recycleRowInfo(rowInfo)
    retVal
end

function getGlobalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}},
        matrix::CSRMatrix{Data, GID, PID, LID}, localRow::GID)::LID where {Data, GID, PID, LID}
    myGraph = matrix.myGraph

    rowInfo = getRowInfoFromGlobalRow(myGraph, GID(globalRow))
    viewRange = 1:rowInfo.numEntries

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return LID(0)
    end

    numElts = rowInfo.numEntries

    @inbounds dataView = getView(matrix, rowInfo)[viewRange]
    @inbounds if isGloballyIndexed(myGraph)
        indsView = getGlobalView(myGraph, rowInfo[viewRange])
        for i in LID(1):numElts
            copy[1][i] = indsView[i]
            copy[2][i] = dataView[i]
        end
    else
        colMap = getColMap(myGraph)
        indsView = getLocalView(myGraph, rowInfo)[viewRange]
        for i in LID(1):numElts
            copy[1][i] = gid(colMap, indsView[i])
            copy[2][i] = dataView[i]
        end
    end

    recycleRowInfo(rowInfo)
    numElts
end

function getLocalRowCopy(matrix::CSRMatrix{Data, GID, PID, LID},
        localRow::Integer
        )::Tuple{AbstractArray{LID, 1}, AbstractArray{Data, 1}} where {
        Data, GID, PID, LID}
    myGraph = matrix.myGraph

    rowInfo = getRowInfo(myGraph, LID(localRow))

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return (LID[], Data[])
    end

    viewRange = 1:rowInfo.numEntries

    @inbounds if isLocallyIndexed(myGraph)
        curLocalIndices = Vector{LID}(getLocalView(myGraph, rowInfo)[viewRange])
    else
        colMap = getColMap(myGraph)
        curGlobalIndices = getGlobalView(myGraph, rowInfo)[viewRange]
        curLocalIndices = @. lid(colMap, curLocalIndices)
    end
    @inbounds curValues = Array{Data}(getView(matrix, rowInfo)[viewRange])

    retVal = (curLocalIndices, curValues)
    recycleRowInfo(rowInfo)
    retVal
end

function getLocalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}},
        matrix::CSRMatrix{Data, GID, PID, LID}, localRow::LID)::LID where {Data, GID, PID, LID}
    myGraph = matrix.myGraph

    rowInfo = getRowInfo(myGraph, LID(localRow))
    viewRange = 1:rowInfo.numEntries

    @boundscheck if rowInfo.localRow == 0
        recycleRowInfo(rowInfo)
        return LID(0)
    end

    numElts = rowInfo.numEntries

    @inbounds dataView = getView(matrix, rowInfo)[viewRange]
    @inbounds if isLocallyIndexed(myGraph)
        indsView = getLocalView(myGraph, rowInfo[viewRange])
        for i in LID(1):numElts
            copy[1][i] = indsView[i]
            copy[2][i] = dataView[i]
        end
    else
        colMap = getColMap(myGraph)
        indsView = getGlobalView(myGraph, rowInfo)[viewRange]

        for i in LID(1):numElts
            copy[1][i] = lid(colMap, indsView[i])
            copy[2][i] = dataView[i]
        end
    end

    recycleRowInfo(rowInfo)
    numElts
end

function getGlobalRowView(matrix::CSRMatrix{Data, GID, PID, LID},
        globalRow::Integer
        )::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}} where {
        Data, GID, PID, LID}
    if isLocallyIndexed(matrix)
        throw(InvalidStateError("The matrix is locally indexed, so cannot return a "
                * "view of the row with global column indices.  Use "
                * "getGlobalRowCopy(...) instead."))
    end
    myGraph = matrix.myGraph

    rowInfo = getRowInfoFromGlobalRow(myGraph, globalRow)
    if rowInfo.localRow != 0 && rowInfo.numEntries > 0
        viewRange = 1:rowInfo.numEntries
        indices = getGlobalView(myGraph, rowInfo)[viewRange]
        values = getView(matrix, rowInfo)[viewRange]
    else
        indices = GID[]
        values = Data[]
    end
    recycleRowInfo(rowInfo)
    (indices, values)
end

function getLocalRowView(matrix::CSRMatrix{Data, GID, PID, LID},
        localRow::Integer
        )::Tuple{AbstractArray{LID, 1}, AbstractArray{Data, 1}} where {
        Data, GID, PID, LID}
    rowInfo = getRowInfo(matrix.myGraph, LID(localRow))
    retVal = getLocalRowView(matrix, rowInfo)
    recycleRowInfo(rowInfo)

    retVal
end

function getLocalRowView(matrix::CSRMatrix{Data, GID, PID, LID},
        rowInfo::RowInfo{LID}
        )::Tuple{AbstractArray{LID, 1}, AbstractArray{Data, 1}} where {
        Data, GID, PID, LID}

    if isGloballyIndexed(matrix)
        throw(InvalidStateError("The matrix is globally indexed, so cannot return a "
                * "view of the row with local column indices.  Use "
                * "getLocalalRowCopy(...) instead."))
    end

    myGraph = matrix.myGraph

    if rowInfo.localRow != 0 && rowInfo.numEntries > 0
        viewRange = LID(1):rowInfo.numEntries
		indices = view(getLocalView(myGraph, rowInfo), viewRange)
        values = view(getView(matrix, rowInfo), viewRange)
    else
        indices = LID[]
        values = Data[]
    end
    (indices, values)
end

Base.@propagate_inbounds @inline function getLocalRowViewPtr(
        matrix::CSRMatrix{Data, GID, PID, LID}, localRow::LID
        )::Tuple{Ptr{LID}, Ptr{Data}, LID} where {Data, GID, PID, LID}
    row = LID(localRow)
    graph = matrix.myGraph

    if getProfileType(graph) == STATIC_PROFILE
        offset1D = graph.rowOffsets[row]
        numEntries = (length(graph.numRowEntries) == 0 ?
                    graph.rowOffsets[row+1] - offset1D
                    : graph.numRowEntries[row])
        if numEntries > 0
            indicesPtr = pointer(graph.localIndices1D, offset1D)
            rowValues = matrix.localMatrix.values
            if rowValues isa SubArray{Data, 1, Vector{Data}, Tuple{UnitRange{LID}}, true}
                valuesPtr = pointer(rowValues.parent, offset1D + rowValues.indexes[1].start)
            elseif rowValues isa Vector{Data}
                #else should be Vector, but assert anyways
                valuesPtr = pointer(rowValues, offset1D)
            else
                error("localMatrix.values is of unsupported type $(typeof(rowValues)).")
            end

            return (indicesPtr, valuesPtr, numEntries)
        else
            return (C_NULL, C_NULL, 0)
        end
    else #dynamic profile
        indices = graph.localIndices2D[row]
        values = matrix.values2D[row]

        return (pointer(indices, 0), pointer(values, 0), length(indices))
    end
end

function getLocalDiagCopy(matrix::CSRMatrix{Data, GID, PID, LID})::MultiVector{Data, GID, PID, LID} where {Data, GID, PID, LID}
    if !hasColMap(matrix)
        throw(InvalidStateError("This method requires a column map"))
    end

    rowMap = getRowMap(matrix)
    colMap = getColMap(matrix)

    numLocalRows = getLocalNumRows(matrix)

    if isFillComplete(matrix)
        diag = MultiVector{Data, GID, PID, LID}(rowMap, 1, false)

        diag1D = getVectorView(diag, 1)
        localRowMap = getLocalMap(rowMap)
        localColMap = getLocalMap(colMap)
        localMatrix = matrix.localMatrix

        diag1D[:] = getDiagCopyWithoutOffsets(matrix, localRowMap, localColMap, localMatrix)

        diag
    else
        getLocalDiagCopyWithoutOffsetsNotFillComplete(matrix)
    end
end

function leftScale!(matrix::CSRMatrix{Data}, X::AbstractArray{Data, 1}) where {Data <: Number}
    for row in 1:getLocalNumRows(matrix)
        _, vals = getLocalRowView(matrix, row)
        LinAlg.scale!(vals, X[row])
    end
end

function rightScale!(matrix::CSRMatrix{Data}, X::AbstractArray{Data, 1}) where {Data <: Number}
    for row in 1:getLocalNumRows(matrix)
        inds, vals = getLocalRowView(matrix, row)
        for entry in 1:length(inds)
            vals[entry] *= X[inds[entry]]
        end
    end
end

#### DistObject methods ####
function checkSizes(source::RowMatrix{Data, GID, PID, LID},
        target::CSRMatrix{Data, GID, PID, LID})::Bool where {Data, GID, PID, LID}
    true
end

function copyAndPermute(source::RowMatrix{Data, GID, PID, LID},
        target::CSRMatrix{Data, GID, PID, LID}, numSameIDs::LID,
        permuteToLIDs::AbstractArray{LID, 1}, permuteFromLIDs::AbstractArray{LID, 1}
        ) where {Data, GID, PID, LID}
    sourceIsLocallyIndexed = isLocallyIndexed(source)

    srcRowMap = getRowMap(source)
    tgtRowMap = getRowMap(target)

    sameGIDs = @. gid(srcRowMap, collect(1:numSameIDs))
    permuteFromGIDs = @. gid(srcRowMap, permuteFromLIDs)
    permuteToGIDs   = @. gid(srcRowMap, permuteToLIDs)

    for (sourceGID, targetGID) in zip(vcat(sameGIDs, permuteFromGIDs), vcat(sameGIDs, permuteToGIDs))
        if sourceIsLocallyIndexed
            rowInds, rowVals = getGlobalRowCopy(source, sourceGID)
        else
            rowInds, rowVals = getGlobalRowView(source, sourceGID)
        end
        combineGlobalValues(target, targetGID, rowInds, rowVals, INSERT)
    end
end

function pack(source::CSRMatrix{Data, GID, PID, LID}, exportLIDs::AbstractArray{LID, 1},
        distor::Distributor{GID, PID, LID})::AbstractArray where {Data, GID, PID, LID}
    numExportLIDs = length(exportLIDs)

    localMatrix = source.localMatrix
    localGraph = localMatrix.graph

    packed = Array{Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}}(numExportLIDs)
    result = 0
    for i in 1:numExportLIDs
        exportLID = exportLIDs[i]
        start = localGraph.rowOffsets[exportLID]
        last = localGraph.rowOffsets[exportLIDs+1]-1
        numEnt = last - start +1
        if numEnt == 0
            packed[i] = GID[], Data[]
        else
            values = view(localMatrix.values, start:last)
            lids = view(localGraph.entries, start:last)
            gids = @. gid(getColMap(source), lids)
            packed[i] = gids, values
        end
    end
    packed
end

function unpackAndCombine(target::CSRMatrix{Data, GID, PID, LID},
        importLIDs::AbstractArray{LID, 1}, imports::AbstractArray, distor::Distributor{GID, PID, LID},
        cm::CombineMode) where{Data, GID, PID, LID}

    numImportLIDs = length(importLIDs)

    for i = 1:numImportLIDs
        if length(imports[i] > 0) #ensure theres acutually something in the row
            combineGlobalValues(target, importLIDs[i], imports[i][1], imports[i][2], cm)
        end
    end
end

#### Operator methods ####

function localApply(Y::MultiVector{Data, GID, PID, LID},
        A::CSRMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID},
        mode::TransposeMode, alpha::Data, beta::Data) where {Data, GID, PID, LID}

    rawY = getLocalArray(Y)
    rawX = getLocalArray(X)

    #TODO implement this better, can BLAS be used?
    if !isTransposed(mode)
        numRows = getLocalNumRows(A)
        for vect = LID(1):numVectors(Y)
            for row = LID(1):numRows
                sum::Data = Data(0)
                @inbounds (indices, values, len) = getLocalRowViewPtr(A, row)
                for i in LID(1):LID(len)
                    ind::LID = unsafe_load(indices, i)
                    val::Data = unsafe_load(values, i)
                    @inbounds sum += val*rawX[ind, vect]
                end
                sum = applyConjugation(mode, sum*alpha)
                @inbounds rawY[row, vect] *= beta
                @inbounds rawY[row, vect] += sum
            end
        end
    else
        rawY[:, :] *= beta
        numRows = getLocalNumRows(A)
        for vect = LID(1):numVectors(Y)
            for mRow in LID(1):numRows
                @inbounds (indices, values, len) = getLocalRowViewPtr(A, mRow)
                for i in LID(1):LID(len)
                    ind::LID = unsafe_load(indices, i)
                    val::Data = unsafe_load(values, i)
                    @inbounds rawY[ind, vect] += applyConjugation(mode, alpha*rawX[mRow, vect]*val)
                end
            end
        end
    end

    Y
end


text:
using PyPlot
using CSV, DataFrames, DataFramesMeta
using HypothesisTests

""" For Fig.2A
Compute average path length across real or simulated (edit!!) participants
"""
function get_pathlength(;
                        filestr = "",
                        # filestr = "_Sim",
                        # savepath = "./",
                        )
    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # ------ Which data? Real or simulated? Edit!

    # Real participants' data:
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEIC_all_fMRI.csv"

    # Simulated data:
    # data_path = "./mcmc_rl_fit/projects/fmri/data_sim/SARSPEICZVG_all_fMRI_Sim_exporder_SurpAC_run2.csv"

    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    # --------------------------------------------------------------------------
    save_path = makehomesavepath("pathlength")
    mkdir(save_path)

    max_episodes_per_subj = 60
    (all_data, all_stats) = RL_Fit.load_SARSPEI(data_path;max_episodes_per_subj=max_episodes_per_subj)
    all_subject_ids = all_stats.all_subject_ids
    println("all_subject_ids: $all_subject_ids")
    nr_of_subjects = length(all_subject_ids)

    pathlength_all_subjects = zeros(nr_of_subjects, max_episodes_per_subj)
    row_counter=1

    for subj_id in all_subject_ids
        subj_data = RL_Fit.filter_by_person(all_data,subj_id)
        idx_reward = findall(subj_data[:,Int(RL_Fit.c_reward)].>0)
        pathlength = diff(vcat(0,idx_reward))
        nr_episodes = length(pathlength)
        pathlength_all_subjects[row_counter, 1:nr_episodes] = pathlength
        pathlength_all_subjects[row_counter, nr_episodes+1:end] .= NaN
        #     figure()
        #     plot(pathlength)
        row_counter+=1

    end
    m, s, serror = calc_stats_nan(pathlength_all_subjects')
    weights = nanlength(pathlength_all_subjects',2)

    dfresults = DataFrame()
    dfresults[!, Symbol("meanlength")] = m
    dfresults[!, Symbol("stdlength")] = s
    dfresults[!, Symbol("selength")] = serror
    dfresults[!, Symbol("weights")] = weights[:]

    CSV.write(joinpath(save_path, "meanpathlength"*filestr*".csv"), dfresults,
        delim = " ");

    dfresults_perepisode = DataFrame()
    for i_ep in 1:max_episodes_per_subj
        dfresults_perepisode[!, Symbol("ep_"*string(i_ep))] = pathlength_all_subjects[:,i_ep]
    end
    CSV.write(joinpath(save_path, "pathlengthperepisode"*filestr*".csv"), dfresults_perepisode,
        delim = " ");

    # For random search comparison:
    # pathlength_firstepisode = pathlength_all_subjects[:,1]
    # pathlength_firstepisode_minusrandom = pathlength_firstepisode .- 6.
    #
    # @show OneSampleTTest(pathlength_firstepisode_minusrandom)
    # pval = pvalue(OneSampleTTest(pathlength_firstepisode_minusrandom))
    # @show pval

    pathlength_all_subjects, m, s, serror
end


text:
@testset "RSSItem" begin
    rss = F.RSSItem(
        "title", "www.link.com", "description", "author@author.com", "category",
        "www.comments.com", "enclosure", Date(2012,12,12))
    @test rss.title == "title"
    @test rss.link == "www.link.com"
    @test rss.description == "description"
    @test rss.author == "author@author.com"
    @test rss.category == "category"
    @test rss.comments == "www.comments.com"
    @test rss.enclosure == "enclosure"
    @test rss.pubDate == Date(2012,12,12)
end

@testset "RSSbasics" begin
    fs1()
    empty!(F.RSS_DICT)
    F.def_GLOBAL_VARS!()
    F.set_var!(F.GLOBAL_VARS, "website_title", "Website title")
    F.set_var!(F.GLOBAL_VARS, "website_descr", "Website descr")
    F.set_var!(F.GLOBAL_VARS, "website_url", "https://github.com/tlienart/Franklin.jl/")
    F.def_LOCAL_VARS!()
    set_curpath("hey/ho.md")
    F.set_var!(F.LOCAL_VARS, "rss_title", "title")
    F.set_var!(F.LOCAL_VARS, "rss", "A **description** done.")
    F.set_var!(F.LOCAL_VARS, "rss_author", "chuck@norris.com")

    item = F.add_rss_item()
    @test item.title == "title"
    @test item.description == "A description done.\n"
    @test item.author == "chuck@norris.com"
    # unchanged bc all three fallbacks lead to Data(1)
    @test item.pubDate == Date(1)

    F.set_var!(F.LOCAL_VARS, "rss_title", "")
    @test @test_logs (:warn, "Found an RSS description but no title for page /hey/ho.html.") F.add_rss_item().title == ""

    @test F.RSS_DICT["/hey/ho.html"].description == item.description

    # Generation
    F.PATHS[:folder] = td
    F.rss_generator()
    feed = joinpath(F.PATHS[:folder], "feed.xml")
    @test isfile(feed)
    fc = prod(readlines(feed, keep=true))
    @test occursin("<description><![CDATA[A description done.", fc)
end

@testset "RSSbasics" begin
    fs2()
    empty!(F.RSS_DICT)
    F.def_GLOBAL_VARS!()
    F.set_var!(F.GLOBAL_VARS, "website_title", "Website title")
    F.set_var!(F.GLOBAL_VARS, "website_descr", "Website descr")
    F.set_var!(F.GLOBAL_VARS, "website_url", "https://github.com/tlienart/Franklin.jl/")
    F.def_LOCAL_VARS!()
    set_curpath("hey/ho.md")
    F.set_var!(F.LOCAL_VARS, "rss_title", "title")
    F.set_var!(F.LOCAL_VARS, "rss", "A **description** done.")
    F.set_var!(F.LOCAL_VARS, "rss_author", "chuck@norris.com")

    item = F.add_rss_item()
    @test item.title == "title"
    @test item.description == "A description done.\n"
    @test item.author == "chuck@norris.com"
    # unchanged bc all three fallbacks lead to Data(1)
    @test item.pubDate == Date(1)

    F.set_var!(F.LOCAL_VARS, "rss_title", "")
    @test @test_logs (:warn, "Found an RSS description but no title for page /hey/ho/index.html.") F.add_rss_item().title == ""

    @test F.RSS_DICT["/hey/ho/index.html"].description == item.description

    # Generation
    F.PATHS[:folder] = td
    F.rss_generator()
    feed = joinpath(F.PATHS[:site], "feed.xml")
    @test isfile(feed)
    fc = prod(readlines(feed, keep=true))
    @test occursin("<description><![CDATA[A description done.", fc)
end


text:
"""
This is a Julia translation of the schema expected by typescript.

Calling [`to_json`](@ref) on a SemagramSchema should give a schema that
typescript can understand.

There is also a macro here for creating a SemagramSchema: [`@semagramschema`](@ref),
which makes writing down a schema much more convenient.
"""
module Schema

export AttributeType, Numeric, Stringlike,
  BoxProperties, PortProperties, WireProperties,
  SemagramSchema, BoxDesc, PortDesc, WireDesc, DataDesc,
  @semagramschema

import ..Muesli: to_json, from_json
using ..Muesli
using ..SVG
using ..Boxes
using Catlab.Present, Catlab.Theories
using MLStyle
using JSExpr

@enum AttributeType begin
  Numeric
  Stringlike
end

to_json(x::AttributeType) = @match x begin
  Numeric => "Numeric"
  Stringlike => "Stringlike"
end

function from_json(d::Any, ::Type{AttributeType})
  @match d begin
    "Numeric" => Numeric
    "Stringlike" => Stringlike
  end
end

@enum EntityType begin
  Box
  Port
  Wire
end

to_json(x::EntityType) = @match x begin
  Box => "Box"
  Port => "Port"
  Wire => "Wire"
end

function from_json(d::Any, ::Type{EntityType})
  @match d begin
    "Box" => Box
    "Port" => Port
    "Wire" => Wire
  end
end

struct OutgoingHom
  name::Symbol
  codom::Symbol
  codom_ty::EntityType
end

to_json(x::OutgoingHom) = generic_to_json(x)
from_json(d::Dict{String, <:Any}, ::Type{OutgoingHom}) = generic_from_json(d, OutgoingHom)

struct BoxProperties
  weights::Vector{Tuple{AttributeType, Symbol}}
  homs::Vector{OutgoingHom}
  shape::String
  label::Union{Symbol,Nothing}
  style_fn::String
end

to_json(x::BoxProperties) = generic_to_json(x)
from_json(d::Dict{String,<:Any}, ::Type{BoxProperties}) = generic_from_json(d, BoxProperties)

struct PortProperties
  weights::Vector{Tuple{AttributeType, Symbol}}
  homs::Vector{OutgoingHom}
  box::Symbol
  box_map::Symbol
  style::String
  style_fn::String
end

to_json(x::PortProperties) = generic_to_json(x)
from_json(d::Dict{String,<:Any}, ::Type{PortProperties}) = generic_from_json(d, PortProperties)

struct WireProperties
  weights::Vector{Tuple{AttributeType, Symbol}}
  homs::Vector{OutgoingHom}
  src::Tuple{EntityType, Symbol}
  src_map::Symbol
  tgt::Tuple{EntityType, Symbol}
  tgt_map::Symbol
  style::String
  style_fn::String
end

to_json(x::WireProperties) = generic_to_json(x)
from_json(d::Dict{String,<:Any}, ::Type{WireProperties}) = generic_from_json(d, WireProperties)

struct SemagramSchema
  box_types::Dict{Symbol, BoxProperties}
  port_types::Dict{Symbol, PortProperties}
  wire_types::Dict{Symbol, WireProperties}
end

to_json(x::SemagramSchema) = generic_to_json(x)
from_json(d::Dict{String, <:Any}, ::Type{SemagramSchema}) = generic_from_json(d, SemagramSchema)

struct BoxDesc
  name::Symbol
  shape::SVGNode
  label::Union{Symbol,Nothing}
  style_fn::JSString
  function BoxDesc(name,shape=Circle; label=nothing, style_fn=@js _ -> $(Dict()))
    new(name,shape,label,style_fn)
  end
end

struct PortDesc
  name::Symbol
  box_map::Symbol
  style::String
  style_fn::JSString
  function PortDesc(name,box_map; style="Circular", style_fn=@js _ -> $(Dict()))
    new(name,box_map,style,style_fn)
  end
end

struct WireDesc
  name::Symbol
  src_map::Symbol
  tgt_map::Symbol
  style::String
  style_fn::JSString
  function WireDesc(name,src_map,tgt_map; style="DefaultWire", style_fn=@js _ -> $(Dict()))
    new(name,src_map,tgt_map,style,style_fn)
  end
end

struct DataDesc
  name::Symbol
  type::AttributeType
end

const EntityDesc = Union{BoxDesc, PortDesc, WireDesc, DataDesc}

function pres_to_semagramschema(p::Presentation, descs::Array)
  ws = SemagramSchema(Dict(),Dict(),Dict())

  datadescs = filter(desc -> typeof(desc) <: DataDesc, descs)
  datatypes = Dict(map(desc -> desc.name => desc.type, datadescs)...)

  obtypes = Dict{Symbol, EntityType}(
    desc.name =>
      if typeof(desc) <: BoxDesc
        Box
      elseif typeof(desc) <: PortDesc
        Port
      else
        Wire
      end
    for desc in descs)

  function weights(ob::Symbol)
    attrs = filter(attr -> nameof(dom(attr)) == ob, p.generators[:Attr])
    [(get(datatypes, nameof(codom(attr)), Stringlike), nameof(attr)) for attr in attrs]
  end

  function homs(ob::Symbol, excluding::Vector{Symbol})
    outgoing = filter(f -> nameof(dom(f)) == ob && !(nameof(f) ∈ excluding), p.generators[:Hom])
    [OutgoingHom(nameof(f), nameof(codom(f)), obtypes[nameof(codom(f))]) for f in outgoing]
  end

  for desc in descs
    ob = desc.name
    if typeof(desc)<:BoxDesc
      ws.box_types[ob] = BoxProperties(
        weights(ob),
        homs(ob,Symbol[]),
        write_svg(desc.shape),
        desc.label,
        string(desc.style_fn),
      )
    elseif typeof(desc)<:PortDesc
      box = nameof(codom(p[desc.box_map]))
      ws.port_types[ob] = PortProperties(
        weights(ob),
        homs(ob, [desc.box_map]),
        box,
        desc.box_map,
        desc.style,
        string(desc.style_fn),
      )
    elseif typeof(desc)<:WireDesc
      src = nameof(codom(p[desc.src_map]))
      tgt = nameof(codom(p[desc.tgt_map]))
      ws.wire_types[ob] = WireProperties(
        weights(ob),
        homs(ob, [desc.src_map, desc.tgt_map]),
        (obtypes[src], src), desc.src_map,
        (obtypes[tgt], tgt), desc.tgt_map,
        desc.style,
        string(desc.style_fn),
      )
    end
  end
  ws
end

q(x) = Expr(:quote,x)

equals_to_kw(x) = x
equals_to_kw(x::Expr) =
  if x.head == :(=)
    Expr(:kw, x.args...)
  else
    x
  end

"""
See the examples/general documentation for how to use this.
"""
macro semagramschema(head,body)
  descs = @match body begin
    Expr(:block,lines...) => begin
      map(lines) do line
        @match line begin
          Expr(:macrocall, mname, _, desc, args...) => begin
            (name,morphs) = @match desc begin
              name::Symbol => (name,[])
              Expr(:call, name, args...) => (name,args)
            end
            constructor = if mname == Symbol("@box")
              BoxDesc
            elseif mname == Symbol("@port")
              PortDesc
            elseif mname == Symbol("@wire")
              WireDesc
            elseif mname == Symbol("@data")
              DataDesc
            end
            Expr(:call, constructor, q(name), q.(morphs)..., equals_to_kw.(args)...)
          end
          _ => missing
        end
      end
    end
    _ => error("the body must be a block")
  end
  descs = filter(x -> !ismissing(x), descs)
  name, pres = @match head begin
    Expr(:call, name, pres) => (name, pres)
    _ => error("the head must be a name and presentation")
  end
  :($(esc(name)) = pres_to_semagramschema($(esc(pres)), $(esc(Expr(:vect, descs...)))))
end

end


text:
using Flux
import Flux: trainable

using AbstractTrees
using AbstractTrees: children

"""
    Model(layer)

Wrapper for Flux layer that recursively wraps children layers.

Use `model.layer` to access wrapped layer

Use `children(model)` to access children wrappers
"""
struct Model{T}
    layer::T
    children::AbstractVector{Model}
    paramattrs::AbstractVector{Symbol}
    settingattrs::AbstractVector{Symbol}
    childrenattrs::AbstractVector{Symbol}
end

Model(layer::T, args...) where T = Model{T}(layer, args...)

(model::Model)(x) = model.layer(x)
trainable(model::Model) = (model.layer,)

totuple(t::Tuple) = t
totuple(t::NTuple) = t
totuple(x) = (x,)

function Model(layer::T) where T
    fields = Dict(
        ParamAttr => [],
        SettingAttr => [],
        ChildAttr => [],
    )
    for name in fieldnames(T)
        # dispatch on only layer type and field name
        at = attrtype(layer, Val(name))
        push!(fields[at], name)
    end

    chs = Model[]
    for field in fields[ChildAttr]
        for child in totuple(getfield(layer, field))
            push!(chs, Model(child))
        end
    end

    return Model{T}(
        layer,
        chs,
        fields[ParamAttr],
        fields[SettingAttr],
        fields[ChildAttr],
    )
end

# Iteration

AbstractTrees.children(model::Model) = model.children

"""
    IterModels(model::Model)

Pre-order (parents before children) iterator over all `Model`s.
If you want to iterate over the wrapped layers instead, use
[`IterLayers`](@ref).

Example
```julia
> model = Model(Chain(Conv((3, 3), 3 => 16), BatchNorm(16)))
> collect(IterModels(model))

[Model{Chain{...}}, Model{Conv{...}}, Model{BatchNorm{...}}]
```
"""
IterModels(model::Model) = PreOrderDFS(model)

"""
    IterLayers(model::Model)

Pre-order (parents before children) iterator over all wrapped layers.

Example
```julia
> model = Model(Chain(Conv((3, 3), 3 => 16), BatchNorm(16)))
> collect(IterModels(model))

[Chain{...}, Conv{...}, BatchNorm{...}]
```
"""
IterLayers(model::Model) = (m.layer for m in IterModels(model))
IterLayers(model) = (m.layer for m in IterModels(Model(model)))

numleaves(tree) = filter(m -> m !== tree, collect(Leaves(tree))) |> length

# Printing

function Base.show(io::IO, li::Model{T}) where T
    s = "$(T.name)"
    if length(li.paramattrs) > 0
        s = string(s, ", params $(Tuple(li.paramattrs))")
    end
    if length(li.settingattrs) > 0
        s = string(s, ", settings $(Tuple(li.settingattrs))")
    end
    if length(li.childrenattrs) > 0
        nch = length(li.children)
        s = string(s, ", $(nch) direct children")
    end
    print(io, s)
end

"""
    printmodel(model)

Print `model` as a tree
"""
printmodel(model::Model) = print_tree(model)
printmodel(model) = print_tree(Model(model))


text:
include("../src/Kinetics.jl");
include("../src/Control.jl");

# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function Balances(t,x,dxdt_vector,data_dictionary)
# ---------------------------------------------------------------------- #
# Balances.jl was generated using the Kwatee code generation system.
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 08-08-2016 20:01:41
#
# Arguments:
# t  - current time
# x  - state vector
# dxdt_vector - right hand side vector
# data_dictionary  - Data dictionary instance (holds model parameters)
# ---------------------------------------------------------------------- #

# Correct nagative x's = throws errors in control even if small -
idx = find(x->(x<0),x);
x[idx] = 0.0;

# Call the kinetics function -
(rate_vector) = Kinetics(t,x,data_dictionary);

# Call the control function -
(rate_vector) = Control(t,x,rate_vector,data_dictionary);

# Encode the balance equations as a matrix vector product -
# Get the required matricies -
maximum_specific_growth_rate = data_dictionary["MAXIMUM_SPECIFIC_GROWTH_RATE"];
S = data_dictionary["STOICHIOMETRIC_MATRIX"];
dilution_selection_array = data_dictionary["DILUTION_SELECTION_ARRAY"];
tau_array = data_dictionary["TIME_CONSTANT_ARRAY"];

# Calculate and correct dxdt -
tmp_vector = S*rate_vector;
number_of_states = length(tmp_vector);
for state_index in collect(1:number_of_states)
	dxdt_vector[state_index] = tau_array[state_index]*(tmp_vector[state_index] - maximum_specific_growth_rate*dilution_selection_array[state_index]*(x[state_index]));
end

# set ribosome and RNAP rates to 0
dxdt_vector[65] = 0;
dxdt_vector[66] = 0;
# set ppase rate to 0
dxdt_vector[60] = 0;

# Correct extracellular species -
cellmass = Cellmass(t,x,data_dictionary);

end

# ---------------------------------------------------------------------- #
# User defined cellmass function -
# ---------------------------------------------------------------------- #
function Cellmass(t,x,data_dictionary)
	# User function to calculate the cellmass abundance -
	# Default behavior is cellmass = 1.0
	cellmass = 1.0;
	return cellmass;
end