text:
# Undergraduate Student: Arturo Burgos
# Professor: João Rodrigo Andrade
# Federal University of Uberlândia - UFU, Fluid Mechanics Laboratory - MFLab, Block 5P, Uberlândia, MG, Brazil

using Printf # Just the same as import a module
using Statistics
using CPUTime

A = rand(1:50,5000,5000,3)

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("The A[1,1,2] is: ",A[1,1,2])
println("The A[1,1,3] is: ",A[1,1,3])

@time @CPUtime A[:, :, 1] = A[:, :, 2]
@time @CPUtime A[:, :, 3] = A[:, :, 1]

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("The A[1,1,2] is: ",A[1,1,2])
println("The A[1,1,3] is: ",A[1,1,3])

max_stars_repo_name: bangyen/Comparing-Languages

max_issues_repo_name: bangyen/Comparing-Languages

max_forks_repo_name: bangyen/Comparing-Languages


text:
#

using GRAMPC

####################################################
# Dynamics functions
####################################################
function dyn( dx::Vector{Float64}, x::Vector{Float64}, u::Vector{Float64},
              p::Vector{Float64}, t::Float64, param::Vector{Float64} )
    dx[1] = -0.04*u[1] + x[2]
    dx[2] = -7.01*u[1]
end

# Jacobian df/dx multiplied by vector vec, i.e. (df/dx)^T*vec or vec^T*(df/dx)
function dfdx_vec( res::Vector{Float64}, x::Vector{Float64}, u::Vector{Float64},
                   p::Vector{Float64}, t::Float64, vec::Vector{Float64}, param::Vector{Float64} )
    res[1] = 0;
    res[2] = vec[1];
end

# Jacobian df/du multiplied by vector vec, i.e. (df/du)^T*vec or vec^T*(df/du)
function dfdu_vec( res::Vector{Float64}, x::Vector{Float64}, u::Vector{Float64},
                   p::Vector{Float64}, t::Float64, vec::Vector{Float64}, param::Vector{Float64} )
    res[1] = -0.04*vec[1] - 7.01*vec[2];
end

####################################################
# Integral cost functions
####################################################
function lfct( t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, xdes::Vector{Float64},
               udes::Vector{Float64}, param::Vector{Float64} )
    cost = (   param[1] * ( x[1] - xdes[1] )^2
             + param[2] * ( x[2] - xdes[2] )^2
             + param[3] * ( u[1] - udes[1] )^2 ) / 2;

    return cost
end

# Gradient dl/dx
function dldx( res::Vector{Float64}, t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64},
               xdes::Vector{Float64}, udes::Vector{Float64}, param::Vector{Float64} )
    res[1] = param[1] * ( x[1] - xdes[1] );
    res[2] = param[2] * ( x[2] - xdes[2] );
end

# Gradient dl/du
function dldu( res::Vector{Float64}, t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64},
               xdes::Vector{Float64}, udes::Vector{Float64}, param::Vector{Float64} )
    res[1] = param[3] * (u[1] - udes[1]);
end

####################################################
# Terminal cost functions
####################################################
function Vfct( T::Float64, x::Vector{Float64}, p::Vector{Float64},xdes::Vector{Float64},
               param::Vector{Float64} )
    cost = (   param[4] * ( x[1] - xdes[1] )^2
             + param[5] * ( x[2] - xdes[2] )^2 ) / 2;

    return cost
end

# Gradient dV/dx
function dVdx( res::Vector{Float64}, T::Float64, x::Vector{Float64}, p::Vector{Float64}, xdes::Vector{Float64},
               param::Vector{Float64} )
    res[1] = param[4] * ( x[1] - xdes[1] );
    res[2] = param[5] * ( x[2] - xdes[2] );
end

####################################################
# Inequality constraints h(t,x(t),u(t),p,uperparam) <= 0
####################################################
function hfct( res::Vector{Float64}, t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64},
               param::Vector{Float64} )
    res[1] =  param[6] - x[1];
    res[2] = -param[7] + x[1];
    res[3] =  param[8] - x[2];
    res[4] = -param[9] + x[2];
end

# Jacobian dh/dx multiplied by vector vec, i.e. (dh/dx)^T*vec or vec^T*(dg/dx)
function dhdx_vec( res::Vector{Float64}, t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64},
                   vec::Vector{Float64}, param::Vector{Float64} )
    res[1] = -vec[1] + vec[2];
    res[2] = -vec[3] + vec[4];
end

# Jacobian dh/du multiplied by vector vec, i.e. (dh/du)^T*vec or vec^T*(dg/du)
function dhdu_vec( res::Vector{Float64}, t::Float64, x::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64},
                   vec::Vector{Float64}, param::Vector{Float64} )
    res[1] = 0;
end

####################################################
# Setup the GRAMPC problem
####################################################
BPmodel = GRAMPC.Model( states = 2, inputs = 1, ineqcon = 4 )

# Add the dynamics to the problem
set_dynamics!( BPmodel, dyn, dx=dfdx_vec, du=dfdu_vec )

# Setup the cost of the problem
set_lagrange_term!( BPmodel, lfct, dx=dldx, du=dldu )
set_mayer_term!( BPmodel, Vfct, dx=dVdx )

# Setup the constraints for the problem
set_inequality_constraints!( BPmodel, hfct, dx=dhdx_vec, du=dhdu_vec )

# Starting state
x0 = [ 0.1;
       0.01 ]

# Starting input
u0 = [ 0.0 ]

# Constraint tolerances
constraintTol = [ 1e-3;
                  1e-3;
                  1e-3;
                  1e-3 ]

# Initialize the problem and setup the options
setup!( BPmodel, Thor=0.3,                        # Length of the prediction horizon
                 Nhor=10,                         # Number of system integration steps
                 dt=0.01,                         # Sampling time
                 t0=0.0,                          # Initial time
                 MaxMultIter=3,                   # Maximum number of augmented lagrangian iterations
                 AugLagUpdateGradientRelTol=1,    # Augmented langrangian gradient tolerance
                 ConstraintsAbsTol=constraintTol, # Constraint tolerances
                 x0=x0,
                 u0=u0 )

# Set the limits after initializing the problem (they are considered an option and not normal constraints)
umax = [  0.0524 ]
umin = [ -0.0524 ]

set_limits!( BPmodel, :input, lb=umin, ub=umax )

# Desired state/input
xdes = [ -0.2;
          0.0 ]
udes = [  0.0 ]

set_setpoint!( BPmodel, input=udes, state=xdes )

# Pass in the user parameters
pCost = [ 100; 10; 180; 100; 10; -0.2; 0.2; -0.1; 0.1 ]
set_user_parameters!( BPmodel, pCost )

estimate_penalty_min!( BPmodel )

# Solve the problem
sol = solve!( BPmodel )

ext: jl


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
using Suppressor
splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

# Writer
# ------

effwriter = @syntax_eff begin
  a = Writer("hello ", 3)
  b = Writer("world!", 5)
  @pure a, b
end
@test effwriter.acc == "hello world!"
@test effwriter.value == (3, 5)

effwriter2 = @syntax_eff begin
  a = Writer("hello ", 3)
  b = collect(a:a+2)
  c = Writer("world!", b*b)
  @pure a, b, c
end

@test effwriter2.acc == "hello world!world!world!"
@test effwriter2.value == [(3, 3, 9), (3, 4, 16), (3, 5, 25)]

effwriter3 = @runhandlers (WriterHandler("PURE"), Vector) @syntax_eff_noautorun begin
  a = Writer("hello ", 3)
  b = collect(a:a+2)
  c = Writer("world!", b*b)
  @pure a, b, c
end
@test effwriter3.acc == "hello world!world!world!PURE"
@test effwriter3.value == [(3, 3, 9), (3, 4, 16), (3, 5, 25)]

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
# AUTO GENERATED FILE - DO NOT EDIT

export iconpickeritem

"""
    iconpickeritem(;kwargs...)

An iconPickerItem component.

Keyword arguments:
- `id` (String; required): The ID used to identify this component in Dash callbacks.
- `icon` (String; required): The name of the icon to render. Example: "FaAdobe"
- `size` (Real; optional): The size of the icon
- `color` (String; optional): The color of the icon. Default: "#000"
- `className` (String; optional): CSS class for div
"""
function iconpickeritem(; kwargs...)
        available_props = Symbol[:id, :icon, :size, :color, :className]
        wild_props = Symbol[]
        return Component("iconpickeritem", "iconPickerItem", "dash_react_font_awesome_icon_picker", available_props, wild_props; kwargs...)
end


text:
function eye(tpe::Type,n::S) where {S <: Integer}

    if tpe <: Number
        return Matrix{tpe}(I,n,n)
    end

end

eye(n::Integer) = eye(Float64,n::Integer)

function tracem(x::Array{T,2}) where {T <: Real}

  # Computes the matrix-trace as defined by Gomme and Klein (2011)

  # We require the number of rows to be greater than the number of columns, so
  # that m is greater than one.

  trans = false

  if size(x,1) < size(x,2)
    x = Matrix(x')
    trans = true
  end

  n = size(x,2)
  m = Int(size(x,1)/n)

  y = zeros(m,1)  # We want this to be a 2-d array for subsequent matrix multiplication

  @inbounds for i = 1:m
    @views y[i,1] = tr(x[(n*(i-1)+1):i*n,1:n])
  end

  if trans == true
    y = y'
  end

  return Matrix(y)

end

function kron_prod_times_vector(a::AbstractArray{T,2},b::AbstractArray{T,2},v::AbstractArray{T,1}) where {T <: Real}

    # This function efficiently computes kron(a,b)*v where a and b are matrices
    # and v is a vector.

    (n1,n2) = size(a)
    (n3,n4) = size(b)
    n5      = length(v)

    if n5 != n2*n4
        error("a, b, and v do not have the correct sizes")
    end

    p = vec(b*reshape(v,n4,n2)*a')

    return p

end

function vector_times_kron_prod(v::AbstractArray{T,1},a::AbstractArray{T,2},b::AbstractArray{T,2}) where {T <: Real}

    product = Matrix(kron_prod_times_vector(a',b',v')')

    return product

end

function kron_prod_times_matrix(a::AbstractArray{T,2},b::AbstractArray{T,2},v::AbstractArray{T,2}) where {T <: Real}

    # This function efficiently computes kron(a,b)*v where a, b, and v are
    # conformable matrices.

    (n1,n2) = size(a)
    (n3,n4) = size(b)
    (n5,n6) = size(v)

    if n5 != n2*n4
        error("a, b, and v do not have the correct sizes")
    end

    p = zeros(n1*n3,n6)
    for i = 1:n6
        @views p[:,i] = kron_prod_times_vector(a,b,v[:,i])
    end

    return p

end

function matrix_times_kron_prod(v::AbstractArray{T,2},a::AbstractArray{T,2},b::AbstractArray{T,2}) where {T <: Real}

    product = Matrix(kron_prod_times_matrix(a',b',v')')

    return product

end

function kron_prod_times_vector(a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2},v::AbstractArray{T,1}) where {T <: Real}

    # This function efficiently computes kron(a,kron(b,c))*v where a, b, and c
    # are matrices and v is a vector.  Included for completeness, but not
    # actually used.

    (n1,n2) = size(a)
    (n3,n4) = size(b)
    (n5,n6) = size(c)
    n7      = length(v)

    if n7 != n2*n4*n6
        error("a, b, c, and v do not have the correct sizes")
    end

    v_tilda = reshape(v,n4*n6,n2)*a'
    p = vec(kron_prod_times_matrix(b,c,v_tilda))

    return p

end

function vector_times_kron_prod(v::AbstractArray{T,1},a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2}) where {T <: Real}

    product = Matrix(kron_prod_times_vector(a',b',c',v')')

    return product

end

function kron_prod_times_matrix(a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2},v::AbstractArray{T,2}) where {T <: Real}

    # This function efficiently computes kron(a,kron(b,c))*v where a, b, c and
    # v are comformable matrices.

    (n1,n2) = size(a)
    (n3,n4) = size(b)
    (n5,n6) = size(c)
    (n7,n8) = size(v)

    if n7 != n2*n4*n6
        error("a, b, c, and v do not have the correct sizes")
    end

    p = zeros(n1*n3*n5,n8)
    for i = 1:n8
        @views v_tilda = reshape(v[:,i],n4*n6,n2)*a'
        p[:,i] .= vec(kron_prod_times_matrix(b,c,v_tilda))
    end

    return p

end

function matrix_times_kron_prod(v::AbstractArray{T,2},a::AbstractArray{T,2},b::AbstractArray{T,2},c::AbstractArray{T,2}) where {T <: Real}

    product = Matrix(kron_prod_times_matrix(a',b',c',v')')

    return product

end

function dsylvester(a::AbstractArray{T,2}, b::AbstractArray{T,2}, c::Union{AbstractArray{T,1},AbstractArray{T,2}}) where {T <: Real}

    #= Uses the Hessenberg-Schur method to find the bounded solution of the
       discrete Sylvester equation:

             X + A*X*B = C

       Based on Golub, Nash, and Van Loan (1979).

    =#

    n = size(a, 1)
    m = size(b, 1)
    x = zeros(size(c))

    (s,u) = schur(Matrix(b'))
    (v,t) = hessenberg(a)

    c = v'*c*u

    j = m
    while j > 0
        j1 = j
        if j == 1
            block_size = 1
        elseif isequal(s[j,j-1],0.0) == false
            block_size = 2
            j -= 1
        else;
            block_size = 1
        end
        @views ajj = kron(s[j:j1,j:j1],t) + I
        @views rhs = vec(c[:,j:j1])
        if j1 < m
            @views rhs2 = t*(x[:,(j+1):m]*s[j:j1,(j+1):m]')
            rhs -= vec(rhs2)
        end
        w = ajj\rhs
        @views x[:,j] = w[1:n]
        if block_size == 2
            @views x[:,j1] = w[(n+1):2*n]
        end
        j -= 1
    end

    x = v*x*u'

    return x

end

function trm(x::AbstractArray{T,2}) where {T <: Real}

    # Computes the matrix trace as defined by Binning (2013).  Used for
    # computing the second-order terms, hss, gss.

    (n1,n2) = size(x)
    k = Int(round(sqrt(n2)))

    y = zeros(n1)
    for i = 1:k
        @views y += x[:,i+(i-1)*k]
    end

    return y

end

function trm2(x::AbstractArray{T,2}) where {T <: Real}

    # Computes the matrix trace as defined by Binning (2013).  Used for
    # computing the third-order terms, hssx, gssx.

    (n1,n2) = size(x)
    k = Int(round(n2^(1//3)))

    y = zeros(n1,k)
    for j = 1:k
        for i = 1:k
            @views y[:,j] += x[:,(j-1)+i+(i-1)*k^2]
        end
    end

    return y

end

function create_omega3(n::S) where {S <: Integer}

    # Creates the combination matrix for a third-order perturbation as defined in Levintal (2017).

    # This function is a simplified version of the create_OMEGA function originally written in
    # Matlab by Oren Levintal for his paper "Fifth Order Perturbation Solution to DSGE Models"
    # published in the Journal of Economic Dynamics and Control, 2017.  Permission to translate
    # this function into Julia and release it within the SolveDSGE module was granted by Oren
    # Levintal on February 5, 2020.

    ind = [1:n^3;]
    M = reshape(ind,1,n,n,n)
    Ix = eye(S,n^3)
    omega3 = (reshape(Ix[:,PermutedDimsArray(M,[1,4,2,3])],n^3,n^3)
           + reshape(Ix[:,PermutedDimsArray(M,[1,2,4,3])],n^3,n^3)
           + reshape(Ix[:,PermutedDimsArray(M,[1,2,3,4])],n^3,n^3))

    return omega3

end

function kron_prod_times_vector(A::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},x::Union{Array{T,1},Array{Complex{T},1}},n::Array{S,1},p::S) where {T <: Real, S <: Integer}

    # Computes y = (A[p] * A[p-1] * ... * A[1] )*x

    N = prod(n[1:p])
    z = copy(x)
    for i = 1:p
        z = (A[i]*reshape(z,n[i],Int(N/n[i])))'
    end
    y = reshape(Matrix(z),N)

    return y

end

function kron_prod_times_matrix(A::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},x::Union{Array{T,2},Array{Complex{T},2}},n::Array{S,1},p::S) where {T <: Real, S <: Integer}

    # Computes y = (A[p] * A[p-1] * ... * A[1] )*x

    N = prod(n[1:p])
    y = Array{Complex{T}}(undef,N,size(x,2))
    for j = 1:size(x,2)
        @views z = x[:,j]
        for i = 1:p
            z = (A[i]*reshape(z,n[i],Int(N/n[i])))'
        end
        y[:,j] .= reshape(Matrix(z),N)
    end

    return y

end

function KPShiftSolve(TT::Union{Array{Array{T,2},1},Array{Array{Complex{T},2},1}},n::Array{S,1},c::Union{Array{T,1},Array{Complex{T},1}},lambda::T,alpha::Union{T,Complex{T}}) where {T <: Real, S <: Integer}

    p = length(n)
    N = prod(n)

    c = copy(c)

    TT[p] = alpha*TT[p]
    if p == 1
        y = (TT[1] + lambda*Matrix{T}(I,n[1],n[1]))\c
    else
        y = Array{Complex{T}}(undef,N)
        mp = Int(N/n[p])
        for i = n[p]:-1:1
            idx = ((i-1)*mp+1):(i*mp)
            y[idx] = KPShiftSolve(TT[1:(p-1)],n[1:(p-1)],c[idx],lambda,TT[p][i,i])
            z = kron_prod_times_vector(TT[1:p],y[idx],n,p-1)
            for j = 1:(i-1)
                jdx = ((j-1)*mp+1):(j*mp)
                c[jdx] = c[jdx] - TT[p][j,i]*z
            end
        end
    end

    return y

end

function martin_van_loan(a::Array{T,2},b::Array{T,2},c::Array{T,2},d::Array{T,2},k::S) where {T <: Real, S <: Integer}

    #= Uses a recursive Schur method to find the bounded solution of the Sylvester
       equation:

       AX + B*X*(Kron^(k)C) = D

       Based on Martin and Van Loan (2006).  This is a simplified implementation
       of their algorithm, but it captures most of the gains over Golub, Nash,
       and van Loan (1979).

    =#

    a = copy(a)
    b = a\copy(b)
    c = copy(c)
    d = a\copy(d)

    (v,s) = hessenberg(b)       # v*s*v' = b
    (t,q) = schur(complex(c'))  # q*t*q' = c'

    v = Matrix(v)

    p = k + 2
    TT = Array{typeof(t)}(undef,p)
    TT[1] = s
    for i = 2:p
        TT[i] = conj(t)
    end

    Q     = fill(q,k+1)
    inv_Q = fill(Matrix(q'),k+1)

    n = fill(size(c,1),p)
    n[1] = size(b,1)
    N = prod(n)

    lambda = 1.0

    e = vec(v'*kron_prod_times_matrix(inv_Q,Matrix(d'),n[2:end],k+1)')

    y = reshape(KPShiftSolve(TT,n,e,lambda,1.0),size(d))
    x = real(Matrix(v*kron_prod_times_matrix(Q,Matrix(y'),n[2:end],k+1)'))

    return x

end

function dlyap(a::Array{T,2}, b::Array{T,2}) where {T <: Real}

    n = size(a,1)
    x = zeros(n,n)
    j = n

    (s,u) = schur(a)
    b = u'b*u

    while j > 0
        j1 = j
        if j == 1
            block = 1
        elseif !isequal(s[j,j-1],0.0)
            block = 2
            j -= 1
        else
            block = 1
        end
        @views lhs = kron(s[j:j1,j:j1],s) - I # I = eye(block*n)
        @views rhs = vec(b[:,j:j1])
        if j1 < n
            @views rhs2 = s*(x[:,(j1+1):n]*s[j:j1,(j1+1):n]')
            rhs += vec(rhs2)
        end
        w = -lhs\rhs
        @views x[:,j] = w[1:n]
        if block == 2
            @views x[:,j1] = w[(n+1):block*n]
        end
        j -= 1
    end

    x = u*x*u'

    return x

end

function ind2sub(i::S,dims::Tuple{S,Vararg{S}}) where {S <: Integer}

    if i < 1 || i > prod(dims)
        error("index is out of bounds.")
    end

    subs = Tuple(CartesianIndices(dims)[i])

    return subs

end

function compute_variances(soln::FirstOrderSolutionStoch)

    hx = soln.hx
    k = soln.k
    gx = soln.gx
    sigma = soln.sigma

    var_states = dlyap(hx,k*sigma*k')
    var_jumps  = gx*var_states*gx'

  return var_states,var_jumps

end

function compute_chebyshev_integrals(eps_nodes,eps_weights,nodes,order,rho,sigma)

  # A simplified implementation of the integration described in Judd et al (2017, section 5.1);
  # it uses the mean integral across the nodes.
  # Using integrals[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
  # leads to something much simplier and equally accurate, but is less recognizably appropriate
  # for the case where the shocks are AR(1) processes and ordinary polynomials are not being used.

  terms_num  = Array{Float64}(undef,length(eps_nodes))
  integrals  = Array{Float64}(undef,order+1,length(nodes))
  integrals2 = Array{Float64}(undef,order+1)
  for i = 1:(order+1)
    integrals2[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
    for j = 1:length(nodes)
      terms_num     .= rho*nodes[j] .+ sqrt(2)*sigma*eps_nodes
      terms_den      = rho*nodes[j]
      terms_num     .= chebyshev_polynomial(i,terms_num)[:,i]
      terms_den      = chebyshev_polynomial(i,terms_den)[i]
      integrals[i,j] = sum((terms_num/terms_den).*eps_weights)*pi^(-1/2)
     end
  end

  nodetoosmall = abs.(nodes) .< sqrt(eps())
  if sum(nodetoosmall) > 0
    if length(nodes) == 1
      integrals[:,1] .= integrals2
    else
      for i = 1:length(nodes)
        if nodetoosmall[i] == 1
          if i == 1
            integrals[:,i] .= integrals[:,i+1]
        elseif i == length(nodes)
            integrals[:,i] .= integrals[:,i-1]
          else
            integrals[:,i] .= (integrals[:,i-1]+integrals[:,i+1])/2
          end
        end
      end
    end
  end

  return reshape(sum(integrals,dims=2)/length(nodes),order+1)

#  integrals = Array{Float64}(undef,order+1)
#  for i = 1:(order+1)
#    integrals[i] = sum(exp.(sqrt(2)*sigma*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
#  end
#  return integrals

end

function compute_smolyak_integrals(eps_nodes,eps_weights,nx,order,grid,RHO,sigma)

  integrals = ones(nx,order+1)
  for j = 1:size(sigma,2)
      nodes = unique(grid[:,j])
      integrals[j,:] .= compute_chebyshev_integrals(eps_nodes,eps_weights,nodes,order,RHO[j,j],sigma[j,j])
  end

#  for j = 1:size(sigma,2)
#      for i = 1:(order+1)
#          integrals[j,i] = sum(exp.(sqrt(2)*sigma[j,j]*(i-1)*eps_nodes).*eps_weights)*pi^(-1/2)
#      end
#  end

  return integrals

end

function weight_scale_factors(eps_nodes,eps_weights,multi_index,nx,grid,RHO,sigma)

  unique_multi_index = sort(unique(multi_index))
  unique_orders = SmolyakApprox.m_i(unique_multi_index).-1

  # Here we construct the base integrals

  base_integrals = Array{Array{Float64,2}}(undef,length(unique_orders))
  for i = 1:length(unique_orders)
    base_integrals[i] = compute_smolyak_integrals(eps_nodes,eps_weights,nx,unique_orders[i],grid,RHO,sigma)
  end

  # Compute the unique polynomial terms from the base polynomials

  unique_base_integrals = Array{Array{Float64,2}}(undef,length(unique_orders))
  for i = length(unique_orders):-1:2
    unique_base_integrals[i] = base_integrals[i][:,size(base_integrals[i-1],2)+1:end]
  end
  unique_base_integrals[1] = base_integrals[1]

  # Construct the first row of the interplation matrix

  new_integrals = unique_base_integrals[multi_index[1,1]][1,:]
  for i = 2:size(multi_index,2)
    new_integrals = kron(new_integrals,unique_base_integrals[multi_index[1,i]][i,:])
  end

  weight_scale_factor = copy(new_integrals)

  # Iterate over nodes, doing the above three steps at each iteration

  for j = 2:size(multi_index,1)

    new_integrals = unique_base_integrals[multi_index[j,1]][1,:]
    for i = 2:size(multi_index,2)
      new_integrals = kron(new_integrals,unique_base_integrals[multi_index[j,i]][i,:])
    end
    weight_scale_factor = [weight_scale_factor; new_integrals]

  end

  return weight_scale_factor

end

function scale_weights(weights,weight_scale_factor)

    scaled_weights = weights.*weight_scale_factor

    return scaled_weights

end

function compute_piecewise_linear_integrals(eps_nodes,eps_weights,sigma)

  integral = 1.0
  #integral = sum(exp.(sqrt(2)*sigma*eps_nodes).*eps_weights)*pi^(-1/2)

  return integral

end


text:
using LinearAlgebra

struct StateSpace
    A::Matrix{Float64}
    B::Vector{Float64}
    C::Matrix{Float64}
    D::Float64

    function StateSpace(a::Matrix, b::Vector, c::Matrix, d::Float64 = 0.)
	return new(a, b, c, d)
    end

    # Creates state space from float values
    function StateSpace(a::Float64, b::Float64, c::Float64, d::Float64 = 0.)
	return new(fill(a, (1,1)), [b], fill(c, (1, 1)), d)
    end
end

function msd(;m::Float64 = 1., b::Float64 = 1., k::Float64 = 1.)
    A = [
	 0 1;
	 (-k/m) (-b/m)
	 ]
    return StateSpace(A, [1., 0.], [1. 0.])
end

# Maybe add low- and high-pass filters in here so that I can easily define them.

function low_pass(wc::Float64)
    return StateSpace(-wc, 1., wc)
end

function high_pass(wc::Float64)
    return StateSpace(-wc, 1., -wc, 1.)
end

mutable struct DiscreteSS
    Ad::Matrix{Float64}
    Bd::Vector{Float64}
    C::Matrix{Float64}
    D::Float64

    function DiscreteSS(ss::StateSpace, dt::Float64)
	this = new()
	this.Ad = exp(ss.A * dt)
	if det(ss.A) == 0
	    this.Bd = dt / 2 * (this.Ad + I) * ss.B
	else
	    this.Bd = inv(ss.A) * (this.Ad - I) * ss.B
	end
	this.C = ss.C
	this.D = ss.D
	return this
    end
end

function get_y(ss::StateSpace, x::Vector{Float64}, u::Float64 = 0.)
    return (ss.C * x)[1, 1] + ss.D * u
end

function get_y(dss::DiscreteSS, x::Vector{Float64}, u::Float64 = 0.)
    return (dss.C * x)[1, 1] + dss.D * u
end

function next_x(dss::DiscreteSS, x::Vector{Float64}, u::Float64 = 0.)
    return dss.Ad * x + dss.Bd * u
end

function rk4(f::Function, x::Vector{Float64}, dt::Float64, params, u::Float64 = 0.)
    k1 = f(x, params, u)
    k2 = f(x + dt * k1 / 2, params, u)
    k3 = f(x + dt * k2 / 2, params, u)
    k4 = f(x + dt * k3, params, u)
    return x + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
end

function rk4(ss::StateSpace, x::Vector{Float64}, dt::Float64, u::Float64 = 0.)
    k1 = ss.A * x + ss.B * u
    k2 = ss.A * (x + dt * k1 / 2) + ss.B * u
    k3 = ss.A * (x + dt * k2 / 2) + ss.B * u
    k4 = ss.A * (x + dt * k3) + ss.B * u
    return x + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
end

function eulers(f::Function, x::Vector{Float64}, dt::Float64, params, u::Float64 = 0.)
    return f(x, params, u) * dt + x
end

function eulers(ss::StateSpace, x::Vector{Float64}, dt::Float64, u::Float64 = 0.)
    return (ss.A * x + ss.B * u) * dt + x
end


text:
using Test

include("definitions/arraypartition.jl")

include("utilities/iterationtools.jl")
include("utilities/lbfgs.jl")
include("utilities/conjugate.jl")

include("problems/test_elasticnet.jl")
include("problems/test_lasso_small.jl")
include("problems/test_lasso_small_v_split.jl")
include("problems/test_lasso_small_h_split.jl")
include("problems/test_linear_programs.jl")
include("problems/test_sparse_logistic_small.jl")
include("problems/test_verbose.jl")
include("problems/test_nonconvex_qp.jl")


text:
module Examples
    using LightGraphs
    using SimpleWeightedGraphs
    medium_pert = SimpleDiGraph(Edge.([
        1 => 2
        1 => 3
        1 => 4
        1 => 5
        4 => 6
        5 => 6
        2 => 7
        3 => 7
        6 => 7
        3 => 8
        6 => 8
        6 => 9
        7 => 10
        8 => 10
        9 => 10
    10 => 11
    ]))

    two_lines = SimpleDiGraph(Edge.([
        1 => 3, 3=>5, 5=>7, 7=>9,
        2 => 4, 4=>6, 6=>8,
    ]))

    loop = SimpleDiGraph(Edge.([
        1 .=> [2, 3];
        2 => 4; 4 => 6;
        3 => 5; 5 => 7;
        [6, 7] .=> 8
    ]))

    cross = SimpleDiGraph(Edge.([
        1 .=> [2, 3];
        2 .=> [3, 4];
        3 => 4;
    ]))

    xcross = SimpleDiGraph(Edge.([
        1 .=> [2, 3, 4];
        2 .=> [3, 4];
        3 => 4;
    ]))

    sankey_3twos = SimpleWeightedDiGraph(6)
    add_edge!(sankey_3twos, 1, 3, 1.0)
    add_edge!(sankey_3twos, 1, 4, 2.0)
    add_edge!(sankey_3twos, 2, 3, 3.0)
    add_edge!(sankey_3twos, 2, 4, 8.0)
    add_edge!(sankey_3twos, 3, 5, 3.0)
    add_edge!(sankey_3twos, 3, 6, 1.0)
    add_edge!(sankey_3twos, 4, 5, 10.0)

    tree = SimpleDiGraph(Edge.([1 => 2, 2 => 3, 4 => 5, 4 => 6, 4 => 7, 4 => 8, 4 => 9, 4 => 10, 5 => 11, 5 => 12, 8 => 15, 8 => 16, 8 => 17, 8 => 18, 8 => 19, 9 => 20, 9 => 21, 10 => 22, 12 => 13, 13 => 14, 23 => 4, 23 => 24, 23 => 25, 23 => 26, 23 => 27, 23 => 28, 23 => 29, 23 => 30, 23 => 31, 28 => 32, 28 => 33, 29 => 35, 30 => 1, 30 => 38, 31 => 40, 33 => 34, 35 => 36, 35 => 37, 38 => 39, 40 => 41, 41 => 42]))

    tiny_depgraph = SimpleDiGraph(Edge.([
        1 => 2
        2 => 4
        2 => 5
        2 => 6
        3 => 2
        3 => 4
        3 => 7
        3 => 5
        3 => 8
        3 => 6
        4 => 9
        4 => 10
        4 => 11
    ]))

    large_depgraph = SimpleDiGraph(Edge.([
        1 => 2
        2 => 3
        3 => 4
        3 => 5
        5 => 6
        2 => 7
        2 => 5
        2 => 8
        8 => 9
        9 => 10
        9 => 11
        11 => 12
        12 => 13
        9 => 12
        8 => 14
        14 => 10
        8 => 15
        8 => 13
        8 => 16
        16 => 17
        17 => 13
        16 => 18
        16 => 19
        16 => 13
        16 => 20
        16 => 21
        8 => 22
        22 => 23
        22 => 24
        8 => 25
        25 => 18
        25 => 16
        25 => 26
        26 => 18
        26 => 19
        26 => 16
        26 => 20
        26 => 21
        25 => 21
        8 => 11
        2 => 27
        27 => 4
        27 => 11
        2 => 28
        28 => 29
        29 => 4
        28 => 8
        28 => 30
        30 => 11
        28 => 27
        28 => 31
        31 => 32
        31 => 33
        33 => 34
        28 => 33
        28 => 35
        35 => 8
        35 => 30
        35 => 33
        35 => 36
        36 => 37
        36 => 38
        35 => 11
        35 => 34
        28 => 4
        28 => 39
        39 => 3
        39 => 5
        39 => 8
        39 => 31
        39 => 33
        39 => 40
        40 => 41
        40 => 26
        28 => 40
        28 => 42
        42 => 8
        42 => 27
        42 => 31
        42 => 33
        42 => 35
        42 => 4
        42 => 43
        43 => 4
        43 => 44
        44 => 14
        44 => 32
        44 => 31
        44 => 22
        44 => 11
        44 => 45
        45 => 31
        45 => 11
        45 => 32
        43 => 31
        43 => 11
        43 => 34
        42 => 11
        28 => 11
        2 => 31
        2 => 33
        2 => 39
        2 => 4
        2 => 46
        46 => 8
        46 => 47
        46 => 31
        46 => 33
        46 => 48
        46 => 4
        46 => 44
        46 => 42
        46 => 11
        2 => 49
        49 => 8
        49 => 32
        49 => 31
        49 => 40
        49 => 44
        49 => 11
        49 => 45
        49 => 38
        2 => 40
        2 => 50
        2 => 44
        2 => 11
        2 => 45
        2 => 26
        2 => 38
    ]))

    extra_large_depgraph = SimpleDiGraph(Edge.([
        1 => 2
        2 => 3
        3 => 4
        3 => 5
        2 => 6
        6 => 7
        2 => 8
        8 => 9
        9 => 10
        10 => 11
        10 => 12
        12 => 13
        13 => 7
        10 => 13
        9 => 14
        14 => 10
        14 => 15
        15 => 11
        14 => 16
        14 => 7
        14 => 17
        17 => 18
        18 => 7
        17 => 19
        17 => 20
        17 => 7
        17 => 21
        17 => 22
        14 => 3
        14 => 23
        23 => 19
        23 => 17
        23 => 24
        24 => 19
        24 => 20
        24 => 17
        24 => 21
        24 => 22
        23 => 22
        14 => 12
        9 => 25
        25 => 26
        25 => 27
        27 => 28
        9 => 17
        9 => 24
        8 => 25
        8 => 14
        8 => 29
        2 => 14
        2 => 30
        30 => 31
        31 => 32
        32 => 33
        31 => 14
        31 => 34
        34 => 12
        31 => 35
        35 => 33
        35 => 12
        31 => 25
        31 => 27
        31 => 36
        36 => 14
        36 => 34
        36 => 27
        36 => 37
        37 => 38
        37 => 39
        36 => 12
        36 => 28
        31 => 33
        31 => 40
        40 => 41
        41 => 33
        41 => 42
        42 => 43
        40 => 42
        40 => 14
        40 => 25
        40 => 27
        40 => 44
        44 => 45
        44 => 24
        31 => 44
        31 => 46
        46 => 14
        46 => 35
        46 => 25
        46 => 27
        46 => 36
        46 => 33
        46 => 47
        47 => 33
        47 => 48
        48 => 15
        48 => 26
        48 => 25
        48 => 3
        48 => 12
        48 => 49
        49 => 25
        49 => 12
        49 => 26
        47 => 25
        47 => 12
        47 => 28
        46 => 12
        31 => 12
        30 => 50
        50 => 41
        50 => 51
        50 => 42
        50 => 14
        50 => 35
        50 => 31
        50 => 25
        50 => 27
        50 => 40
        50 => 33
        50 => 52
        52 => 14
        52 => 53
        52 => 25
        52 => 27
        52 => 54
        52 => 33
        52 => 48
        52 => 46
        52 => 12
        50 => 55
        55 => 14
        55 => 26
        55 => 25
        55 => 44
        55 => 48
        55 => 12
        55 => 49
        55 => 39
        50 => 44
        50 => 56
        50 => 48
        50 => 12
        50 => 49
        50 => 24
        50 => 39
        30 => 57
        57 => 41
        57 => 42
        57 => 14
        57 => 34
        57 => 35
        57 => 25
        57 => 58
        57 => 27
        57 => 36
        57 => 33
        57 => 40
        57 => 47
        57 => 59
        59 => 41
        59 => 60
        59 => 42
        57 => 46
        57 => 12
        30 => 51
        30 => 25
        30 => 40
        30 => 52
        30 => 61
        61 => 51
        61 => 42
        61 => 25
        61 => 62
        62 => 41
        62 => 59
        62 => 63
        62 => 42
        62 => 64
        62 => 65
        65 => 66
        61 => 44
        2 => 67
        2 => 68
        2 => 69
        69 => 25
        69 => 70
        69 => 12
        69 => 49
        2 => 71
        71 => 7
        71 => 72
        72 => 23
        72 => 17
        2 => 44
        2 => 17
        2 => 73
        73 => 17
        73 => 74
        74 => 17
        73 => 25
        2 => 75
        75 => 31
        75 => 50
        75 => 57
        75 => 51
        75 => 25
        75 => 76
        76 => 41
        76 => 42
        76 => 14
        76 => 25
        76 => 27
        76 => 77
        77 => 42
        77 => 78
        77 => 79
        79 => 43
        77 => 80
        80 => 41
        80 => 42
        80 => 81
        81 => 33
        81 => 42
        81 => 65
        80 => 60
        80 => 33
        80 => 65
        80 => 59
        77 => 12
        77 => 43
        76 => 36
        76 => 33
        76 => 40
        76 => 44
        76 => 46
        76 => 12
        75 => 52
        75 => 61
        2 => 39
        2 => 23
        2 => 25
        2 => 48
        2 => 82
        82 => 10
        82 => 50
        82 => 14
        82 => 45
        82 => 25
        82 => 7
        82 => 83
        83 => 84
        83 => 26
        83 => 80
        83 => 44
        83 => 12
        82 => 55
        82 => 80
        82 => 56
        82 => 12
        82 => 49
        82 => 39
        2 => 19
        2 => 50
        2 => 85
        85 => 6
        2 => 86
        86 => 10
        86 => 69
        86 => 25
        86 => 68
        86 => 17
        86 => 23
        86 => 24
        2 => 72
        2 => 7
        2 => 24
    ]))
end  # module


text:
include("serial.jl")
include("../../src/star.jl")

#solves Ax = b in parallel, where A is a symmetric tridiagonal matrix
#upper triangular solve
function usolve(y, z, c, d, l, x)
  #this code performs 5n flops, 4 reads, and 1 write
  n = size(x, 1)
  #y, z are solutions at the previous and next separators
  x[n] = z
  @inbounds for i = n - 1 : -1 : 1
    x[i] = (x[i] - c[i + 1] * x[i + 1] - l[i] * y) / d[i] # 5 flops
  end
end

#lower triangular solve
function lsolve(a, b, c, d, l, x)
  #this code performs 12n flops in total
  #it performs 3 reads, and 3 writes
  n = size(x,1)
  d[1] = a[1]
  x[1] = b[1]
  f = c[1]
  l[1] = f
  r = f / d[1]
  s = f * r
  y = x[1] * r

  @inbounds for i = 2 : n - 1
    l_ = c[i] / d[i - 1]         # 1
    d[i] = a[i] - c[i] * l_      # 2
    x[i] = b[i] - x[i - 1] * l_  # 2
    f *= -l_                     # 2 including negation for fillin
    l[i] = f                     # store fillin
    r = f / d[i]                 # 1
    s += f * r                   # 2
    y += x[i] * r                # 2
  end
  l_ = c[n] / d[n - 1]
  d[n] = a[n] - c[n] * l_
  x[n] = b[n] - x[n - 1] * l_
  f = - r * c[n]  # fillin between previous and next separators

  s, y, f, d[n], x[n]
end

function root(s, y, f, d, x)
  #solve the reduced tridiagonal system
  #d is the diagonal of the reduced system, f the subdiagonal, and
  #x the reduced right hand side
  d[1:end - 1] -= s[2:end]
  x[1:end - 1] -= y[2:end]
  serial_trid(d, x, f[2:end])
  n = length(x)
  x = [0; x]
  return [(x[i], x[i+1]) for i = 1 : n]
end

trid_star(a, b, c, d, l, x) = Star.scan(usolve,
                    root(Star.reduce(lsolve, a, b, c, d, l, x)...), c, d, l, x)

#test routine
function trid_star_test(n)
  P = nworkers() # of processors
  #setup the inputs and outputs
  Da = DArray(I->100 * rand(map(length,I)), (n,), workers(), P, f1)
  Db = DArray(I->rand(map(length,I)), (n,), workers(), P, f1)
  Dc = DArray(I->rand(map(length,I)), (n + 1,), workers(), P, f2)
  Dd = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dl = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dx = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)

  t = @elapsed trid_star(Da, Db, Dc, Dd, Dl, Dx)

  # code to test the results
  x = convert(Array{Float64,1}, Dx)
  a = convert(Array{Float64,1}, Da)
  b = convert(Array{Float64,1}, Db)
  c = convert(Array{Float64,1}, Dc)
  d = convert(Array{Float64,1}, Dd)
  c = sub(c, 2:n)
  println("norm(A * x -b)", norm(matrix_vector_product(a, c, c, x) -b))
  println("serial trid - not storing L")
  x[:] = 0 ;
  t1 = @elapsed serial_trid(a, b, c, d, x)
  println("norm(A * x -b)", norm(matrix_vector_product(a, c, c, x) -b))
  println("serial time ", t1, " parallel time ", t)
end

function trid_star_run(n, num_runs = 1)
  P = nworkers() # of processors
  #setup the inputs and outputs
  Da = DArray(I->100 * rand(map(length,I)), (n,), workers(), P, f1)
  Db = DArray(I->rand(map(length,I)), (n,), workers(), P, f1)
  Dc = DArray(I->rand(map(length,I)), (n + 1,), workers(), P, f2)
  Dd = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dl = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dx = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  T = zeros(num_runs)

  for nruns = 1 : num_runs
    t = @elapsed trid_star(Da, Db, Dc, Dd, Dl, Dx)
    T[nruns] = t
  end
  println(" parallel time ", T)
  println(" min ", minimum(T), " max ", maximum(T), " mean ",  mean(T))
end

function test_serial(n, num_runs = 1)
  P = nworkers() # of processors
  #setup the inputs and outputs
  a = 100 * rand(n)
  b = rand(n)
  c = rand(n)
  d = zeros(n)
  x = zeros(n)
  T = zeros(num_runs)
  for nruns = 1 : num_runs
    t = @elapsed serial_trid(a, b, c, d, x)
    T[nruns] = t
  end
  println("serial time ", T)
  println(" min ", minimum(T), " max ", maximum(T), " mean ",  mean(T))
end

function trid_no_star(n, num_runs = 1)
  P = nworkers() # of processors
  #setup the inputs and outputs
  Da = DArray(I->100 * rand(map(length,I)), (n,), workers(), P, f1)
  Db = DArray(I->rand(map(length,I)), (n,), workers(), P, f1)
  Dc = DArray(I->rand(map(length,I)), (n + 1,), workers(), P, f2)
  Dd = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dl = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)
  Dx = DArray(I->zeros(map(length,I)), (n,), workers(), P, f1)

  ref = Array(Any, P)
  T = zeros(num_runs)
  for nruns = 1 : num_runs
    t = @elapsed begin
      for i = 1 : P
        ref[i] = @spawnat workers()[i] lsolve(localpart(Da), localpart(Db),
                    localpart(Dc), localpart(Dd), localpart(Dl), localpart(Dx))
      end
      s = Array(Float64, P)
      y = Array(Float64, P)
      f = Array(Float64, P)
      d = Array(Float64, P)
      x = Array(Float64, P)
      results = [fetch(ref[i]) for i = 1 : P]
      for i = 1 : P
        (s[i], y[i], f[i], d[i], x[i]) = results[i]
      end
      X = root(s, y, f, d, x)
      @sync begin
        for i = 1 : P
          @spawnat workers()[i] usolve(X[i][1], X[i][2], localpart(Dc),
                                  localpart(Dd), localpart(Dl), localpart(Dx))
        end
      end
    end #@elapsed
    T[nruns] = t
  end #for
  println(" parallel time ", T)
  println(" min ", minimum(T), " max ", maximum(T), " mean ",  mean(T))
#=
  # code to test the results
  x = convert(Array{Float64,1}, Dx)
  a = convert(Array{Float64,1}, Da)
  b = convert(Array{Float64,1}, Db)
  c = convert(Array{Float64,1}, Dc)
  d = convert(Array{Float64,1}, Dd)
  c = sub(c, 2:n)
  println("norm(A * x -b)", norm(matrix_vector_product(a, c, c, x) -b))
  println("serial trid - not storing L")
  x[:] = 0 ;
  t1 = @elapsed serial_trid(a, b, c, d, x)
  println("norm(A * x -b)", norm(matrix_vector_product(a, c, c, x) -b))
  println("serial time ", t1, " parallel time ", t)
=#
end

#auxiliary functions to partition the darrays among processors.
function f1(dims, procs) # auxiliary function to partition all darrays but the
                         # subdiagonal
  n = dims[1]
  P = size(procs, 1)
  leaf_size::Int64 = max(floor(n / P), 2) #size of a leaf
  cuts = ones(Int, P + 1) * leaf_size
  cuts[1] = 1
  cuts = cumsum(cuts)
  cuts[end] = n + 1
  indexes = Array(Any, P)
  for i = 1:P
    indexes[i] = (cuts[i] : cuts[i+1] - 1,)
  end
  indexes, Any[cuts]
end

function f2(dims, procs) # auxiliary function to partition the subdiagonal
  n = dims[1]
  P = size(procs, 1)
  leaf_size::Int64 = max(floor((n - 1) / P), 2) #size of a leaf
  cuts = ones(Int, P + 1) * leaf_size
  cuts[1] = 1
  cuts = cumsum(cuts)
  cuts[end] = n + 1
  indexes = Array(Any, P)
  for i = 1:P
    indexes[i] = (cuts[i] : cuts[i+1] - 1,)
  end
  indexes, Any[cuts]
end

#overload the DArray constructor to provide user defined partitions.
#=
function DArray(init, dims, procs, dist, distfunc::Function)
    np = prod(dist)
    procs = procs[1:np]
    idxs, cuts = distfunc([dims...], procs)

    chunks = Array(RemoteRef, dist...)
    for i = 1:np
        chunks[i] = remotecall(procs[i], init, idxs[i])
    end
    return DistributedArrays.construct_darray(dims, chunks, procs, idxs, cuts)
end
=#

#function DArray(init, dims, procs, dist)
function DArray(init, dims, procs, dist, distfunc::Function)
    np = prod(dist)
    procs = reshape(procs[1:np], ntuple(i->dist[i], length(dist)))
    idxs, cuts = distfunc([dims...], procs)
    #idxs, cuts = chunk_idxs([dims...], dist)
    identity = DistributedArrays.next_did()

    return DistributedArrays.construct_darray(identity, init, dims, procs, idxs, cuts)
end


text:
using RecurrenceAnalysis
using DynamicalSystemsBase, Random, Statistics, SparseArrays
using Graphs, LinearAlgebra
using Test
using DataStructures
RA = RecurrenceAnalysis

rng = Random.seed!(194)
# Trajectories of 200 points
# Examples of the Hénon map based on:
#     Webber CL & Zbilut JP, "Recurrence Quantification Analysis of Nonlinear
#     Dynamical Systems", in: Riley MA & Van Orden GC, Tutorials in Contemporary
#     Nonlinear Methods for the Behavioral Sciences, 2005, 26-94.
#     https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf
#
trajectories = Dict(
    "Sine wave" => RA.Dataset(map(x->[sin.(x) cos.(x)], StepRangeLen(0.0,0.2,200))),
    "White noise" => RA.Dataset(randn!(zeros(200,2))),
    "Hénon (chaotic)" => trajectory(Systems.henon(a=1.4, b=0.3), 199, Ttr=1000),
    "Hénon (periodic)" => trajectory(Systems.henon(a=1.054, b=0.3), 199, Ttr=1000)
)
embed_params = Dict(       #(m, τ)
    "Sine wave"   => (9, 7),
    "White noise" => (1, 1),
    "Hénon (chaotic)" => (3, 1),
    "Hénon (periodic)" => (3, 1)
)
rqa_threshold = Dict(
    "Sine wave"   => 0.15,
    "White noise" => 0.15,
    "Hénon (chaotic)" => 0.15,
    "Hénon (periodic)" => 0.15
)

dict_keys = ["Sine wave","White noise","Hénon (chaotic)","Hénon (periodic)"]
@testset "$k" for k in dict_keys
    data = trajectories[k]
    x = data[:,1].*0.00000001.*randn(length(data[:,1]))
    y = data[1:100,2]
    if k ≠ "White noise"
        xe = embed(x, embed_params[k]...)
        ye = embed(y, embed_params[k]...)
    else
        xe = x
        ye = copy(x)
    end

    # Distance and recurrence matrices
    ε = rqa_threshold[k]
    dmat = distancematrix(xe, ye)
    crmat = CrossRecurrenceMatrix(xe, ye, ε)
    @test Matrix(crmat) == (dmat .≤ ε)
    rmat = RecurrenceMatrix(xe, ε; parallel = false)
    rmat_p = RecurrenceMatrix(xe, ε; parallel = true)
    jrmat = JointRecurrenceMatrix(xe, ye, ε; parallel = false)
    sz = size(jrmat)
    @test (jrmat.data .& rmat.data[1:sz[1],1:sz[2]]) == jrmat.data
    @test (jrmat.data .& rmat_p.data[1:sz[1],1:sz[2]]) == jrmat.data
    # Compare metrics
    m_euc = rmat
    m_euc_p = rmat_p
    m_max = RecurrenceMatrix(xe, ε, metric="max")
    m_min = RecurrenceMatrix(xe, ε, metric="manhattan")
    @test (m_max.data .& m_euc.data) == m_euc.data
    @test (m_euc.data .& m_min.data) == m_min.data
    @test (m_max.data .& m_euc_p.data) == m_euc.data
    @test (m_euc_p.data .& m_min.data) == m_min.data
    # Compare scales
    m_scalemax = CrossRecurrenceMatrix(xe, ye, 0.1, scale=maximum)
    m_scalemean = CrossRecurrenceMatrix(xe, ye, 0.1, scale=mean)
    m_scalemax_p = CrossRecurrenceMatrix(xe, ye, 0.1, scale=maximum, parallel = true)
    m_scalemean_p = CrossRecurrenceMatrix(xe, ye, 0.1, scale=mean, parallel = true)
    @test (m_scalemax.data .& m_scalemean.data) == m_scalemean.data
    @test (m_scalemax_p.data .& m_scalemean_p.data) == m_scalemean.data
    # Fixed rate for recurrence matrix
    crmat_fixed = CrossRecurrenceMatrix(xe, ye, 0.05; fixedrate=true)
    crmat_fixed_p = CrossRecurrenceMatrix(xe, ye, 0.05; fixedrate=true, parallel = true)
    @test .04 < recurrencerate(crmat_fixed) < .06
    @test .04 < recurrencerate(crmat_fixed_p) < .06
    @test recurrencerate(crmat_fixed) ≈ recurrencerate(crmat_fixed_p)
    # fan method for recurrence threshold
    cr_fan = CrossRecurrenceMatrix{FAN}(xe, ye, 0.05; fixedrate=true, parallel = false)
    cr_fan_p = CrossRecurrenceMatrix{FAN}(xe, ye, 0.05; parallel = true)
    rp_fan = RecurrenceMatrix{FAN}(xe, 0.05; fixedrate=true, parallel = false)
    rp_fan_p = RecurrenceMatrix{FAN}(Dataset(xe), 0.05; parallel = true)
    n = length(xe)
    @test all(.04 < nnz(cr_fan[:,i])/n < .06 for i in axes(cr_fan, 2))
    @test all(.04 < (nnz(rp_fan[:,i])-1)/n < .06 for i in axes(rp_fan, 2))
    @test .04 < recurrencerate(cr_fan) < .06
    @test .04 < recurrencerate(cr_fan_p) < .06
    @test .04 < recurrencerate(rp_fan) < .06
    @test .04 < recurrencerate(rp_fan_p) < .06
    @test cr_fan == cr_fan_p
    @test rp_fan == rp_fan_p
    @test recurrencerate(cr_fan) ≈ recurrencerate(cr_fan_p)
    @test recurrencerate(rp_fan_p; theiler=0) > recurrencerate(rp_fan_p)

    # Recurrence plot
    crp = grayscale(crmat, width=125)
    szplot = size(crp)
    szmat  = size(crmat)
    @test szplot[1] ≈ szplot[2]*szmat[2]/szmat[1] atol = 1

    # RQA
    rqapar = rqa(rmat, theiler=1, lmin=3, border=20)
    rqadiag = rqa(rmat, theiler=1, lmin=3, border=20, onlydiagonal=true)
    for p in keys(rqadiag)
        @test rqapar[p]==rqadiag[p]
    end

    # Windowed RQA
    rmatw = @windowed RecurrenceMatrix(xe, ε, metric=RecurrenceAnalysis.Chebyshev()) 50
    @windowed RecurrenceMatrix{FAN}(xe, ε) 50 # not meaningful, only to check that it does not error
    crmatw = @windowed(CrossRecurrenceMatrix(xe, ye, ε),30)
    @windowed jrmatw = JointRecurrenceMatrix(xe, ye, ε) 30
    @test jrmatw[3 .+ (1:30), 3 .+ (1:30)] == jrmat[3 .+ (1:30), 3 .+ (1:30)]
    @windowed(rrw = recurrencerate(rmatw), width=50, step=40)
    @windowed rqaw = rqa(rmatw) width=50 step=40
    @test rqaw[:RR] == rrw

    # Recurrence networks
    graph = SimpleGraph(rmat)
    amat = adjacency_matrix(graph)
    @test amat == Matrix(rmat) - I
    rna_dict = rna(rmat)
    rna_dict[:density] == recurrencerate(rmat)

    #test sparse matrices
    rqapar_sparse = rqa(rmat.data, theiler=1, lmin=3, border=20)
    rqadiag_sparse = rqa(rmat.data, theiler=1, lmin=3, border=20, onlydiagonal=true)
    for p in keys(rqadiag_sparse)
        @test rqapar_sparse[p]==rqadiag_sparse[p]
    end

    #check for erroring here with sparse matrices
    @windowed(rrw2 = recurrencerate(rmat.data,theiler=1), width=50, step=40)
    @windowed(rqaw2 = rqa(rmat.data,theiler=1), width=50, step =40)
    @test rqaw2[:RR]==rrw2
    rc1 = coordinates(rmat)
    rc2 = coordinates(rmat.data)
    @test isequal(rc1[1],rc2[1])
    @test isequal(rc2[2],rc2[2])

    @test isa(rqa(rmat.data).data,Dict)
    @test isa(rqa(Dict,rmat.data),Dict)
    @test isa(rqa(OrderedDict,rmat.data),OrderedDict)
end


text:
using SuperSub
using Base.Test

@test superscript(52) == "⁵²"
@test subscript(23852943//29387423) == "₂₃₈₅₂₉₄₃.₂₉₃₈₇₄₂₃"


text:
# Mapping from Julia (as given by MOI) to SCIP operators
const OPMAP = Dict{Symbol, SCIP_ExprOp}(
    :+ => SCIP_EXPR_SUM,       # n-ary
    :* => SCIP_EXPR_PRODUCT,   # n-ary
    :- => SCIP_EXPR_MINUS,     # unary, binary
    :/ => SCIP_EXPR_DIV,       # unary
    :^ => SCIP_EXPR_REALPOWER, # binary (or INTPOWER)
    :sqrt => SCIP_EXPR_SQRT,   # unary
    :exp => SCIP_EXPR_EXP,     # unary
    :log => SCIP_EXPR_LOG,     # unary
    :abs => SCIP_EXPR_ABS,     # unary
    :min => SCIP_EXPR_MIN,     # binary
    :max => SCIP_EXPR_MAX,     # binary
)

"""Subexpressions and variables referenced in an expression tree.

Used to convert Julia expression to SCIP expression using recursive calls to the
mutating push_expr!.
"""
mutable struct NonlinExpr
    vars::Vector{Ptr{SCIP_VAR}}
end

"Extract operators from Julia expr recursively and convert to SCIP expressions."
function push_expr!(nonlin::NonlinExpr, mscip::ManagedSCIP, expr::Expr)
    # Storage for SCIP_EXPR*
    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)
    num_children = length(expr.args) - 1

    if Meta.isexpr(expr, :comparison) # range constraint
        # args: [lhs, <=, mid, <=, rhs], lhs and rhs constant
        @assert length(expr.args) == 5
        @assert expr.args[2][1] == expr.args[4][1] == :<=

        # just call on middle expression, bounds are handled outside
        return push_expr!(nonlin, mscip, expr.args[3])

    elseif Meta.isexpr(expr, :call) # operator
        op = expr.args[1]
        if op in [:(==), :<=, :>=]
            # args: [op, lhs, rhs], rhs constant
            @assert length(expr.args) == 3

            # just call on lhs expression, bounds are handled outside
            return push_expr!(nonlin, mscip, expr.args[2])

        elseif op == :^
            # Special case: power with constant exponent.
            # The Julia expression considers the base and exponent to be
            # subexpressions. SCIP does in principle support constant
            # expressions, but in the case of SCIP_EXPR_REALPOWER, the exponent
            # value is stored directly, as the second child.
            @assert num_children == 2

            # Base (first child) is proper sub-expression.
            base = push_expr!(nonlin, mscip, expr.args[2])

            # Exponent (second child) is stored as value.
            @assert isa(expr.args[3], Number)
            if isa(expr.args[3], Integer)
                exponent = Cint(expr.args[3])
                @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, SCIP_EXPR_INTPOWER, base, exponent)
            else
                exponent = Cdouble(expr.args[3])
                @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, SCIP_EXPR_REALPOWER, base, exponent)
            end

        elseif op == :- && num_children == 1
            # Special case: unary version of minus. SCIP only supports binary
            # minus, so we will represent it as :(0 - child).

            # First, insert constant 0 subexpression:
            left = push_expr!(nonlin, mscip, 0.0)

            # Then, insert the actual subexpression:
            right = push_expr!(nonlin, mscip, expr.args[2])

            # Finally, add the (binary) minus:
            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], left, right)

        elseif op in [:sqrt, :exp, :log, :abs]
            # Unary operators
            @assert num_children == 1

            # Insert child expression:
            child = push_expr!(nonlin, mscip, expr.args[2])

            # Add this operator on top
            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], child)

        elseif op in [:-, :/, :min, :max]
            # Binary operators
            @assert num_children == 2

            # Create left and right subexpression.
            left = push_expr!(nonlin, mscip, expr.args[2])
            right = push_expr!(nonlin, mscip, expr.args[3])

            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], left, right)

        elseif op in [:+, :*]
            # N-ary operators
            @assert num_children >= 1

            # Create all children
            children = Ptr{SCIP_EXPR}[]
            for subexpr in expr.args[2:end]
                child = push_expr!(nonlin, mscip, subexpr)
                push!(children, child)
            end
            @assert length(children) == num_children

            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], Cint(num_children), children)

        else
            error("Operator $op (in $expr) not supported by SCIP.jl!")
        end

    elseif Meta.isexpr(expr, :ref) # variable
        # It should look like this:
        # :(x[MathOptInterface.VariableIndex(1)])
        @assert expr.args[1] == :x
        @assert num_children == 1
        vi = expr.args[2] # MOI.VariableIndex
        vr = VarRef(vi.value)
        v = var(mscip, vr)

        # 0-based indexing for SCIP
        op = SCIP_EXPR_VARIDX
        @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, op, Cint(length(nonlin.vars)))
        push!(nonlin.vars, v)

    else
        error("Expression $expr not supported by SCIP.jl!")
    end

    # Return this expression to be referenced by parent expressions
    return expr__[]
end

function push_expr!(nonlin::NonlinExpr, mscip::ManagedSCIP, expr::Number)
    # Storage for SCIP_EXPR*
    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)

    op = SCIP_EXPR_CONST
    value = Cdouble(expr)

    @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, op, value)

    # double check whether value was saved correctly
    value_stored = SCIPexprGetOpReal(expr__[])
    if value != value_stored
        error("Failed to create SCIP_EXPR / $(op):\n" *
              "passed in constant value $(value),\n" *
              "retrieved $(value_stored) instead,\n" *
              "in expression stored at $(expr__[]).")
    end

    return expr__[]
end

"""
Add nonlinear constraint to problem, return cons ref.

    lhs ≤ expression ≤ rhs

# Arguments
- `expr::Expr`: Julia expression of nonlinear function, as given by MOI.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

"""
function add_nonlinear_constraint(mscip::ManagedSCIP, expr::Expr, lhs::Float64, rhs::Float64)
    nonlin = NonlinExpr([])

    # convert expression recursively, extract root and variable pointers
    root = push_expr!(nonlin, mscip, expr)
    vars = nonlin.vars

    # create expression graph object
    tree__ = Ref{Ptr{SCIP_EXPRTREE}}(C_NULL)
    @SC SCIPexprtreeCreate(SCIPblkmem(mscip), tree__, root, length(vars), 0, C_NULL)
    @SC SCIPexprtreeSetVars(tree__[], length(vars), vars)

    # create and add cons_nonlinear
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicNonlinear(mscip, cons__, "", 0, C_NULL, C_NULL,
                                     1, tree__, C_NULL, lhs, rhs)
    @SC SCIPaddCons(mscip, cons__[])

    # free memory
    @SC SCIPexprtreeFree(tree__)

    # register and return cons ref
    return store_cons!(mscip, cons__)
end


text:
# Autogenerated wrapper script for libfyaml_jll for x86_64-unknown-freebsd
export fy_compose, fy_dump, fy_filter, fy_join, fy_testsuite, fy_tool, fy_ypath, libfyaml_c

JLLWrappers.@generate_wrapper_header("libfyaml")
JLLWrappers.@declare_library_product(libfyaml_c, "libfyaml.so.0")
JLLWrappers.@declare_executable_product(fy_compose)
JLLWrappers.@declare_executable_product(fy_dump)
JLLWrappers.@declare_executable_product(fy_filter)
JLLWrappers.@declare_executable_product(fy_join)
JLLWrappers.@declare_executable_product(fy_testsuite)
JLLWrappers.@declare_executable_product(fy_tool)
JLLWrappers.@declare_executable_product(fy_ypath)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfyaml_c,
        "lib/libfyaml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        fy_compose,
        "bin/fy-compose",
    )

    JLLWrappers.@init_executable_product(
        fy_dump,
        "bin/fy-dump",
    )

    JLLWrappers.@init_executable_product(
        fy_filter,
        "bin/fy-filter",
    )

    JLLWrappers.@init_executable_product(
        fy_join,
        "bin/fy-join",
    )

    JLLWrappers.@init_executable_product(
        fy_testsuite,
        "bin/fy-testsuite",
    )

    JLLWrappers.@init_executable_product(
        fy_tool,
        "bin/fy-tool",
    )

    JLLWrappers.@init_executable_product(
        fy_ypath,
        "bin/fy-ypath",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Random
using LinearAlgebra
Random.seed!(843)

include("feast.jl")
using ..feastLinear
using IterativeSolvers

n=100

alpha=0.99 # #determine how nonsymmetric our matrix is; alpha=1.0 makes a symmetric matrix
#generate B-matrix eigenvectors and eigenvalues
#vectors:
R=rand(n,n)
(Xb,R)=qr(R)
Xb=Xb+(1-alpha)*ones(n,n)
Yb=inv(Xb)'
#values:
Lb=rand(n)*0.5 .+ 0.5 #random positive numbers

#generate A-matrix eigenvectors:
R=rand(n,n)
(X,R)=qr(R)
X=alpha*X+(1-alpha)*ones(n,n)
Y=inv(X)'

#generate A-matrix eigenvalues
L1=zeros(n)
for i in 1:n
    L1[i]=i
end

#generate eigenvalue problem matrices
A=Xb*diagm(0 => sqrt.(Lb))*X*diagm(0 => L1)*Y'*diagm(0 => sqrt.(Lb))*Yb'
B=Xb*diagm(0 => Lb)*Yb'

#FEAST parameters
emid=1.0+0.0*im #contour center
ra=0.5 #contour radius 1
rb=0.5 #contour radius 2
nc=8 #number of contour points
m0=2 #subspace dimension
eps=1e-9 #residual convergence tolerance
maxit=100 #maximum FEAST iterations
x0=rand(n,m0) #eigenvector initial guess
y0=copy(x0)

alphals=1e-16 #linear system solve accuracy
isMaxit=30 #maximum number of linear system iterations
#(ye,le,xe)=nseig(A,B)
println("Standard FEAST:")
(ye,le,xe,data)=feastNS_linear(A,B,x0,y0,nc,emid,ra,rb,eps,maxit; log=true)
println("Iterative FEAST:")
(ye,le,xe)=ifeastNS_linear(A,B,x0,y0,alphals,isMaxit,nc,emid,ra,rb,eps,maxit)

p=sortperm(real(le))
println(le[p])

max_stars_repo_name: brendanedwardgavin/feastjl

max_issues_repo_name: brendanedwardgavin/feastjl

max_forks_repo_name: brendanedwardgavin/feastjl


text:
using Mimi
using GraphPlot
using Compose

function _open_file(filename)
    if Sys.isapple()
        run(`open $(filename)`)
    elseif Sys.islinux()
        run(`xdg-open $(filename)`)
    elseif Sys.iswindows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    else
        @warn "Showing plots is not supported on $(Sys.KERNEL)"
    end
end

"""
    plot_comp_graph(m::Model, filename::Union{Nothing, Symbol} = nothing)

Plot the DAG of component connections within model `m` and save to `filename`. If
no `filename` is given, plot will simply display.
"""
function plot_comp_graph(m::Model, filename::Union{Nothing, String} = nothing)

    graph = comp_graph(m.md)
    paths = map(i -> get_prop(graph, i, :path), vertices(graph))
    names = map(path -> path.names[end], paths)

    plot = gplot(graph, nodelabel=names, nodesize=6, nodelabelsize=6)
    if filename !== nothing
        draw(PDF(filename, 16cm, 16cm), plot)
        return _open_file(filename)
    else
        return plot
    end
end


text:
#!/usr/bin/env julia
#
# Simple DMRG tutorial.  This code contains a basic implementation of the
# infinite system algorithm
#
# Copyright 2013-2015 James R. Garrison and Ryan V. Mishmash.
# Open source under the MIT license.  Source code at
# <https://github.com/simple-dmrg/simple-dmrg/>

using Arpack
using LinearAlgebra
using SparseArrays

# Data structures to represent the block and enlarged block objects.
struct Block
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol,AbstractMatrix{Float64}}
end

struct EnlargedBlock
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol,AbstractMatrix{Float64}}
end

# For these objects to be valid, the basis size must match the dimension of
# each operator matrix.
isvalid(block::Union{Block,EnlargedBlock}) =
    all(op -> size(op) == (block.basis_size, block.basis_size), values(block.operator_dict))

# Model-specific code for the Heisenberg XXZ chain
model_d = 2  # single-site basis size

Sz1 = [0.5 0.0; 0.0 -0.5]  # single-site S^z
Sp1 = [0.0 1.0; 0.0 0.0]  # single-site S^+

H1 = [0.0 0.0; 0.0 0.0]  # single-site portion of H is zero

function H2(Sz1, Sp1, Sz2, Sp2)  # two-site part of H
    # Given the operators S^z and S^+ on two sites in different Hilbert spaces
    # (e.g. two blocks), returns a Kronecker product representing the
    # corresponding two-site term in the Hamiltonian that joins the two sites.
    J = 1.0
    Jz = 1.0
    return (J / 2) * (kron(Sp1, Sp2') + kron(Sp1', Sp2)) + Jz * kron(Sz1, Sz2)
end

# conn refers to the connection operator, that is, the operator on the edge of
# the block, on the interior of the chain.  We need to be able to represent S^z
# and S^+ on that site in the current basis in order to grow the chain.
initial_block = Block(1, model_d, Dict{Symbol,AbstractMatrix{Float64}}(
    :H => H1,
    :conn_Sz => Sz1,
    :conn_Sp => Sp1,
))

function enlarge_block(block::Block)
    # This function enlarges the provided Block by a single site, returning an
    # EnlargedBlock.
    mblock = block.basis_size
    o = block.operator_dict

    # Create the new operators for the enlarged block.  Our basis becomes a
    # Kronecker product of the Block basis and the single-site basis.  NOTE:
    # `kron` uses the tensor product convention making blocks of the second
    # array scaled by the first.  As such, we adopt this convention for
    # Kronecker products throughout the code.
    I1 = sparse(1.0I, model_d, model_d)
    I_block = sparse(1.0I, mblock, mblock)
    enlarged_operator_dict = Dict{Symbol,AbstractMatrix{Float64}}(
        :H => kron(o[:H], I1) + kron(I_block, H1) + H2(o[:conn_Sz], o[:conn_Sp], Sz1, Sp1),
        :conn_Sz => kron(I_block, Sz1),
        :conn_Sp => kron(I_block, Sp1),
    )

    return EnlargedBlock(block.length + 1,
                         block.basis_size * model_d,
                         enlarged_operator_dict)
end

function rotate_and_truncate(operator, transformation_matrix)
    # Transforms the operator to the new (possibly truncated) basis given by
    # `transformation_matrix`.
    return transformation_matrix' * (operator * transformation_matrix)
end

function single_dmrg_step(sys::Block, env::Block, m::Int)
    # Performs a single DMRG step using `sys` as the system and `env` as the
    # environment, keeping a maximum of `m` states in the new basis.

    @assert isvalid(sys)
    @assert isvalid(env)

    # Enlarge each block by a single site.
    sys_enl = enlarge_block(sys)
    if sys === env  # no need to recalculate a second time
        env_enl = sys_enl
    else
        env_enl = enlarge_block(env)
    end

    @assert isvalid(sys_enl)
    @assert isvalid(env_enl)

    # Construct the full superblock Hamiltonian.
    m_sys_enl = sys_enl.basis_size
    m_env_enl = env_enl.basis_size
    sys_enl_op = sys_enl.operator_dict
    env_enl_op = env_enl.operator_dict
    I_sys_enl = sparse(1.0I, m_sys_enl, m_sys_enl)
    I_env_enl = sparse(1.0I, m_env_enl, m_env_enl)
    superblock_hamiltonian = kron(sys_enl_op[:H], I_env_enl) + kron(I_sys_enl, env_enl_op[:H]) +
                             H2(sys_enl_op[:conn_Sz], sys_enl_op[:conn_Sp], env_enl_op[:conn_Sz], env_enl_op[:conn_Sp])

    # Call ARPACK to find the superblock ground state.  (:SR means find the
    # eigenvalue with the "smallest real" value.)
    #
    # But first, we explicitly modify the matrix so that it will be detected as
    # Hermitian by `eigs`.  (Without this step, the matrix is effectively
    # Hermitian but won't be detected as such due to small roundoff error.)
    superblock_hamiltonian = (superblock_hamiltonian + superblock_hamiltonian') / 2
    (energy,), psi0 = eigs(superblock_hamiltonian, nev=1, which=:SR)

    # Construct the reduced density matrix of the system by tracing out the
    # environment
    #
    # We want to make the (sys, env) indices correspond to (row, column) of a
    # matrix, respectively.  Since the environment (column) index updates most
    # quickly in our Kronecker product structure, psi0 is thus row-major.
    # However, Julia stores matrices in column-major format, so we first
    # construct our matrix in (env, sys) form and then take the transpose.
    psi0 = transpose(reshape(psi0, (env_enl.basis_size, sys_enl.basis_size)))
    rho = Hermitian(psi0 * psi0')

    # Diagonalize the reduced density matrix and sort the eigenvectors by
    # eigenvalue.
    fact = eigen(rho)
    evals, evecs = fact.values, fact.vectors
    permutation = sortperm(evals, rev=true)

    # Build the transformation matrix from the `m` overall most significant
    # eigenvectors.
    my_m = min(length(evals), m)
    indices = permutation[1:my_m]
    transformation_matrix = evecs[:, indices]

    truncation_error = 1 - sum(evals[indices])
    println("truncation error: ", truncation_error)

    # Rotate and truncate each operator.
    new_operator_dict = Dict{Symbol,AbstractMatrix{Float64}}()
    for (name, op) in sys_enl.operator_dict
        new_operator_dict[name] = rotate_and_truncate(op, transformation_matrix)
    end

    newblock = Block(sys_enl.length, my_m, new_operator_dict)

    return newblock, energy
end

function infinite_system_algorithm(L::Int, m::Int)
    block = initial_block
    # Repeatedly enlarge the system by performing a single DMRG step, using a
    # reflection of the current block as the environment.
    while 2 * block.length < L
        println("L = ", block.length * 2 + 2)
        block, energy = single_dmrg_step(block, block, m)
        println("E/L = ", energy / (block.length * 2))
    end
end

infinite_system_algorithm(100, 20)

max_stars_repo_name: simple-dmrg/simple-dmrg


text:
module Plasm

#export centroid, cuboidGrid, mkpol, view, hpc_exploded, lar2hpc

using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using DataStructures
using SparseArrays
using PyCall

p = pyimport("pyplasm")

import Base.view

"""
	Points = Matrix

Alias declation of LAR-specific data structure.
Dense `Matrix` ``M times N`` to store the position of *vertices* (0-cells)
of a *cellular complex*. The number of rows ``M`` is the dimension
of the embedding space. The number of columns ``N`` is the number of vertices.
"""
const Points = Matrix

"""
	Cells = Array{Array{Int,1},1}

Alias declation of LAR-specific data structure.
Dense `Array` to store the indices of vertices of `P-cells`
of a cellular complex.
The linear space of `P-chains` is generated by `Cells` as a basis.
Simplicial `P-chains` have `P+1` vertex indices for `cell` element in `Cells` array.
Cuboidal `P-chains` have ``2^P`` vertex indices for `cell` element in `Cells` array.
Other types of chain spaces may have different numbers of vertex indices for `cell`
element in `Cells` array.
"""
const Cells = Array{Array{Int,1},1}

"""
	Chain = SparseVector{Int8,Int}

Alias declation of LAR-specific data structure.
Binary `SparseVector` to store the coordinates of a `chain` of `N-cells`. It is
`nnz=1` with `value=1` for the coordinates of an *elementary N-chain*, constituted by
a single *N-chain*.
"""
const Chain = SparseVector{Int8,Int}

"""
	ChainOp = SparseMatrixCSC{Int8,Int}

Alias declation of LAR-specific data structure.
`SparseMatrix` in *Compressed Sparse Column* format, contains the coordinate
representation of an operator between linear spaces of `P-chains`.
Operators ``P-Boundary : P-Chain -> (P-1)-Chain``
and ``P-Coboundary : P-Chain -> (P+1)-Chain`` are typically stored as
`ChainOp` with elements in ``{-1,0,1}`` or in ``{0,1}``, for
*signed* and *unsigned* operators, respectively.
"""
const ChainOp = SparseMatrixCSC{Int8,Int}

"""
	ChainComplex = Array{ChainOp,1}

Alias declation of LAR-specific data structure. It is a
1-dimensional `Array` of `ChainOp` that provides storage for either the
*chain of boundaries* (from `D` to `0`) or the transposed *chain of coboundaries*
(from `0` to `D`), with `D` the dimension of the embedding space, which may be either
``R^2`` or ``R^3``.
"""
const ChainComplex = Array{ChainOp,1}

"""
	LARmodel = Tuple{Points,Array{Cells,1}}

Alias declation of LAR-specific data structure.
`LARmodel` is a pair (*Geometry*, *Topology*), where *Geometry* is stored as
`Points`, and *Topology* is stored as `Array` of `Cells`. The number of `Cells`
values may vary from `1` to `N+1`.
"""
const LARmodel = Tuple{Points,Array{Cells,1}}

"""
	LAR = Tuple{Points,Cells}

Alias declation of LAR-specific data structure.
`LAR` is a pair (*Geometry*, *Topology*), where *Geometry* is stored as
`Points`, and *Topology* is stored as `Cells`.
"""
const LAR = Tuple{Points,Cells}

"""
	Hpc = PyCall.PyObject

Alias declation of LAR-specific data structure.
`Hpc` stands for *Hierarchical Polyhedral Complex* and is the geometric data structure
used by `PLaSM` (Programming LAnguage for Solid Modeling). See the Wiley's book
[*Geometric Programming for Computer-Aided Design*]
(https://onlinelibrary.wiley.com/doi/book/10.1002/0470013885) and its
current `Python` library [*https://github.com/plasm-language/pyplasm*]
(https://github.com/plasm-language/pyplasm).
"""
const Hpc = PyCall.PyObject

"""
	const color::OrderedDict(String,PyObject)

Dictionary of predefined `color::PyObject` functions, with `string` keys and
`PyObject` value, usable as Julia functions of type `Hpc -> Hpc`.
Notice that in order to be visualized by `PyPlasm` viewer, color functions must
be applied to `Hpc` values, and return `Hpc` values.

# Example

```julia
julia> julia> @show keys(colors);
keys(colors) = ["orange", "red", "green", "blue", "cyan", "magenta",
				"yellow", "white", "black", "purple", "gray", "brown"]

julia> green = colors["green"]

julia> V,CV = Lar.cuboid([1,4,9])
([0.0 0.0 … 1.0 1.0; 0.0 0.0 … 4.0 4.0; 0.0 9.0 … 0.0 9.0], Array{Int64,1}[[1, 2, 3, 4, 5, 6, 7, 8]])

julia> hpc = Plasm.lar2hpc(V,CV)
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type 'std::shared_ptr< Hpc > *' at 0x12943c690> >

julia> green(hpc)
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type 'std::shared_ptr< Hpc > *' at 0x12943c6c0> >

julia> Plasm.view(Plasm.green(hpc))
```
"""
function color(key::String)::PyObject
	p = PyCall.pyimport("pyplasm")
	colors = OrderedDict([
	"white" => p["COLOR"](p["WHITE"]),
	"red" => p["COLOR"](p["RED"]),
	"green" => p["COLOR"](p["GREEN"]),
	"blue" => p["COLOR"](p["BLUE"]),
	"cyan" => p["COLOR"](p["CYAN"]),
	"magenta" => p["COLOR"](p["MAGENTA"]),
	"yellow" => p["COLOR"](p["YELLOW"]),
	"black" => p["COLOR"](p["BLACK"]),
	"orange" => p["COLOR"](p["ORANGE"]),
	"purple" => p["COLOR"](p["PURPLE"]),
	"gray" => p["COLOR"](p["GRAY"]),
	"brown" => p["COLOR"](p["BROWN"]) ])
	return  colors[key]
end

const colorkey = ["white", "red", "green", "blue", "cyan", "magenta", "yellow",
	"black", "orange", "purple", "gray", "brown"]

"""
	Color4f

`PyObject` <function COLOR at 0x126d6f730>

# Example

```julia
julia> black = Color4f(0,0,0,1)
PyObject Color4f(0.000000e+00,0.000000e+00,0.000000e+00,1.000000e+00)

julia> b = color(black)
PyObject <function COLOR.<locals>.COLOR0 at 0x135dfcb70>

julia> white = Color4f(1,1,1,1)
PyObject Color4f(1.000000e+00,1.000000e+00,1.000000e+00,1.000000e+00)

julia> w = color(white)
PyObject <function COLOR.<locals>.COLOR0 at 0x135dfcc80>
```
"""
function Color4f(red::Number, green::Number, blue::Number, alpha::Number)::PyCall.PyObject
	p = PyCall.pyimport("pyplasm")
	return p["Color4f"](red, green, blue, alpha)
end

"""
	Color(color4f::PyCall.PyObject)::PyCall.PyObject

Generate an object of type `PyObject Color4f()` to be used with `pyplasm` package.

# Example

```julia
julia> black = Plasm.Color4f(0,0,0,1)
PyObject Color4f(0.000000e+00,0.000000e+00,0.000000e+00,1.000000e+00)

julia> b = color(black)
PyObject <function COLOR.<locals>.COLOR0 at 0x135dfcb70>

julia> white = Plasm.Color4f(1,1,1,1)
PyObject Color4f(1.000000e+00,1.000000e+00,1.000000e+00,1.000000e+00)

julia> w = Plasm.color(white)
PyObject <function COLOR.<locals>.COLOR0 at 0x135dfcc80>
```
"""
function Color(color4f::PyCall.PyObject)::PyCall.PyObject
	p = PyCall.pyimport("pyplasm")
	return p["COLOR"](color4f)
end

#function color(mycolor::"String", model::Lar.LAR)::PyCall.PyObject
#	p = PyCall.pyimport("pyplasm")
#	hpc = lar2hpc(model)
#	return color(mycolor)(hpc)
#end

include("./fenvs.jl")
include("./views.jl")
include("./graphic_text.jl")
include("./utilities.jl")

#include("./viewer/Viewer.jl")

end # module


text:
if VERSION < v"1.1.0-DEV.472"
    export isnothing
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end


text:
### fnt.jl
#
# Copyright (C) 2020 Andrey Oskin.
# Copyright (C) 2020 Jakub Wronowski.
#
# Maintainer: Jakub Wronowski <jakubwro@users.noreply.github.com>
# Keywords: number theoretic transform, fermat number transform
#
# This file is a part of NumberTheoreticTransforms.jl.
#
# License is MIT.
#
### Commentary:
#
# This file contains implementation of Fermat Number Transform.
#
### Code:

export fnt, fnt!, ifnt, ifnt!, modfermat, isfermat

"""
    isfermat(n)

Checks if a given number is a Fermat number \$ 2^{2^t}+1 \$.
"""
function isfermat(number::T) where {T<:Integer}
    if !ispow2(number - one(T))
        return false
    end

    return ispow2(trailing_zeros(number - one(T)))
end

"""
    modfermat(n, q)

Equivalent of mod(n, q) but uses faster algorithm.
Constraints:
-   `q` must be a Fermat number \$ 2^{2^t}+1 \$
-   `n` must be smaller or equal to \$ (q-1)^2 \$
-   `n` must be grater or equal to \$ 0 \$
If above constraints are not met, the result is undefined.
"""
function modfermat(n::T, q::T) where T <: Integer
    x = n & (q - T(2)) - n >>> trailing_zeros(q - T(1)) + q
    x = x >= q ? x - q : x
end

"""
Order input to perform radix-2 structured calculation.
It sorts array by bit-reversed 0-based sample index.
"""
function radix2sort!(data::Array{T, 1}) where {T<:Integer}
    N = length(data)
    @assert ispow2(N)

    l = 1
    for k in 1:N
        if l > k
            data[l], data[k] = data[k], data[l]
        end

        l = l - 1
        m = N
        while l & (m >>= 1) != 0
            l &= ~m
        end
        l = (l | m) + 1
    end

    return data
end

function fnt!(x::Array{T, 1}, g::T, q::T) where {T<:Integer}
    N = length(x)
    @assert ispow2(N)
    @assert isfermat(q)
    @assert all(v -> 0 <= v < q, x)

    radix2sort!(x)

    logN = trailing_zeros(N)

    for l in 1:logN
        M = 1 << (l-1) # [1,2,4,8,...,N/2]
        interval = 1 << l # [2,4,8,...,N]
        p = 1 << (logN - l) # [N/2,...,4,2,1]
        gp = powermod(g, p, q)
        W = 1
        for m in 1:M
            for i in m:interval:N
               j = i + M
               xi, xj = x[i], x[j]
               Wxj = modfermat(W * xj, q)
               xi, xj = xi + Wxj, xi - Wxj + q
               xi = xi >= q ? xi - q : xi # mod q
               xj = xj >= q ? xj - q : xj # mod q
               x[i], x[j] = xi, xj
            end
            W = modfermat(W * gp, q)
        end
    end

    return x
end

"""
    fnt!(x, g, q)

In-place version of `fnt`. That means it will store result in the `x` array.
"""
function fnt!(x::Array{T,2}, g::T, q::T) where {T<:Integer}
    N, M = size(x)
    @assert N == M #TODO: make it work for N != M (need different g for each dim)

    for n in 1:N
        x[n, :] = fnt!(x[n, :], g, q)
    end

    for m in 1:M
        x[:, m] = fnt!(x[:, m], g, q)
    end

    return x
end

"""
    fnt(x, g, q)

The Fermat Number Transform returns the same result as `ntt` function using
more performant algorithm. When `q` has \$ 2^{2^t}+1 \$ form the calculation
can be performed with O(N*log(N)) operation instead of O(N^2) for `ntt`.
"""
function fnt(x::Array{T}, g::T, q::T) where {T<:Integer}
    return fnt!(copy(x), g, q)
end

"""
    ifnt!(y, g, q)

In-place version of `ifnt`. That means it will store result in the `y` array.
"""
function ifnt!(y::Array{T,1}, g::T, q::T) where {T<:Integer}
    N = length(y)
    inv_N = invmod(N, q)
    inv_g = invmod(g, q)

    x = fnt!(y, inv_g, q)
    for i in eachindex(x)
        x[i] = mod(inv_N * x[i], q)
    end

    return x
end

function ifnt!(y::Array{T,2}, g::T, q::T) where {T<:Integer}
    N, M = size(y)

    for m in 1:M
        y[:, m] = ifnt!(y[:, m], g, q)
    end

    for n in 1:N
        y[n, :] = ifnt!(y[n, :], g, q)
    end

    return y
end

"""
    ifnt(y, g, q)

Calculates inverse of Fermat Number Transform for array `y` using
mod \$ 2^{2^t}+1 \$ arithmetic.

The input must be array of integers calculated by `fnt` function with the same
`g` and `q` params.
"""
function ifnt(y::Array{T}, g::T, q::T) where {T<:Integer}
    return ifnt!(copy(y), g, q)
end


text:
typealias XT Array{Float64,1}
typealias YT ASCIIString

const labels = YT[ "sunny", "rainy", "foggy" ]

function features(yt::YT, x::AbstractArray{XT}, t::Int)
    features("", yt, x, t)
end

function features(yp::YT, yt::YT, x::AbstractArray{XT}, t::Int)
    res = Features(6+9+36+33+84)

    for yval in labels
        @append! res (yp == yval)
        @append! res (yt == yval)
    end

    for ypval in labels, ytval in labels
        @append! res ((yp == ypval) & (yt == ytval))
    end

    for ypval in labels, ytval in labels
        @append! res ((yp == ypval) & (yt == ytval) & (x[t][1] > 0.0))
        @append! res ((yp == ypval) & (yt == ytval) & (x[t][2] > 0.0))
        @append! res ((yp == ypval) & (yt == ytval) & (x[t][1] < 0.0))
        @append! res ((yp == ypval) & (yt == ytval) & (x[t][2] < 0.0))
    end

    for ytval in labels, diff = 0.0:10.0
        @append! res ((yt == ytval) & (diff <= abs(x[t][1] - x[t][2]) < (diff + 1.0)))
    end

    for yval in labels, tres = 0.0:13.0
        @append! res ((yt == yval) & (tres <= abs(x[t][1]) < (tres + 1.0)))
        @append! res ((yt == yval) & (tres <= abs(x[t][2]) < (tres + 1.0)))
    end

    return convert(Array{Float64,1}, res)
end


text:
abstract type LesserElement <: Element end

mutable struct BabelCall <: LesserElement
    name::SubString{String}
end
BabelCall(name::String) = BabelCall(SubString(name))

abstract type Block <: LesserElement end

mutable struct CommentBlock <: Block
    contents::Vector{SubString{String}}
end
mutable struct ExampleBlock <: Block
    contents::Vector{SubString{String}}
end
mutable struct ExportBlock <: Block
    backend::SubString{String}
    contents::Vector{SubString{String}}
end
mutable struct SourceBlock <: Block
    lang::Union{SubString{String}, Nothing}
    arguments::Union{SubString{String}, Nothing}
    contents::Vector{SubString{String}}
end
mutable struct VerseBlock <: Block
    contents::Vector{LesserElement}
end

mutable struct Clock{T <: Union{TimestampInactive, TimestampInactiveRange}} <: LesserElement
    timestamp::T
    duration::Union{Nothing, Tuple{Integer, Integer}}
end

mutable struct DiarySexp <: LesserElement
    sexp::SubString{String}
end

mutable struct Planning <: LesserElement
    deadline::Union{Timestamp, Nothing}
    scheduled::Union{Timestamp, Nothing}
    closed::Union{Timestamp, Nothing}
end

mutable struct Comment <: LesserElement
    contents::Vector{SubString{String}}
end

mutable struct FixedWidth <: LesserElement
    contents::Vector{SubString{String}}
end

struct HorizontalRule <: LesserElement end

mutable struct Keyword{V <: Union{SubString{String}, Vector{Object}, Nothing}} <: LesserElement
    key::SubString{String}
    value::V
end
Keyword(key::String, value::Union{Vector{<:Object}, Nothing}=nothing) =
    Keyword(SubString(key), value)
Keyword(key::String, value::String) = Keyword(SubString(key), SubString(value))
Keyword((key, value)::Pair{String, <:Union{String, Vector{<:Object}, Nothing}}) =
    Keyword(key, value)

mutable struct AffiliatedKeyword{V <: Union{SubString{String}, Vector{Object}, Nothing}} <: LesserElement
    key::SubString{String}
    optval::Union{SubString{String}, Vector{Object}, Nothing}
    value::V
end
AffiliatedKeyword(key::AbstractString, value) =
    AffiliatedKeyword(SubString(key), nothing,
                      if value isa String; SubString(value) else value end)
AffiliatedKeyword((key, value)::Pair{<:AbstractString, <:Any}) =
    AffiliatedKeyword(key, value)

mutable struct LaTeXEnvironment <: LesserElement
    name::SubString{String}
    contents::Vector{SubString{String}}
end
LaTeXEnvironment(name::String, contents::Vector{String}) =
    LaTeXEnvironment(SubString(name), SubString.(contents))
LaTeXEnvironment(name::String, contents::String) =
    LaTeXEnvironment(SubString(name), split(contents, '\n'))

mutable struct NodeProperty <: LesserElement
    name::SubString{String}
    additive::Bool
    value::SubString{String}
end
NodeProperty(name::String, value::String, additive::Bool=false) =
    NodeProperty(SubString(name), additive, SubString(value))
NodeProperty((name, value)::Pair{String, String}, additive::Bool=false) =
    NodeProperty(name, value, additive)

mutable struct Paragraph <: LesserElement
    contents::Vector{Object}
end
Paragraph(obj::Object) = Paragraph(Object[obj])
Paragraph(text::String) = Paragraph(TextPlain(text))
Paragraph() = Paragraph(Object[])

mutable struct TableRow <: LesserElement
    cells::Vector{TableCell}
end
TableRow(cellvals::Vector=[]) =
    TableRow(TableCell.(Vector{Union{String, Vector{Object}}}(cellvals)))

struct TableHrule <: LesserElement end


text:
function myf(x)
    return [x[1]+x[2], 5*x[1]^2+3*x[2]-5]
end

function myj(x)
    return [1  1;
            10*x[1]  3]
end

x = newtonraphson(myf, [1,1], myj)


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function normpatterns(patlist::Vector{Vector{String}}, config::TextConfig)
    A = Vector{String}[]
    for patterns in patlist
        B = String[]
        for pat in patterns
            T = tokenize(pat, config)
            length(T) > 0 && push!(B, join(T, ' '))
        end
        length(B) > 0 && push!(A, B)
    end

    Set(A)
end

function normpatternsset(patlist::Vector{Vector{String}}, config::TextConfig)
    A = Vector{Set{String}}[]

    for patterns in patlist
        B = Set{String}[]
        for pat in patterns
            T = tokenize(pat, config)
            length(T) > 0 && push!(B, Set(T))
        end
        length(B) > 0 && push!(A, B)
    end

    Set(A)
end

function intersectmatch(tokens, pat)
    for w in pat
        if !(w in tokens)
            return false
        end
    end
    return true
end

function intersectmatch(fun, wrap, input, config::TextConfig, key, patlist::Vector{Vector{String}})
    patlist = normpatternsset(patlist, config)

    itertweets(input) do tweet
        tokens = tokenize(tweet[key], config) |> Set

        matches = 0
        for (i, patterns) in enumerate(patlist)
            for pat in patterns
                if intersectmatch(tokens, pat)
                    matches +=1
                    break
                end
            end
            matches != i && break
        end

        wrap(matches == length(patlist)) && fun(tweet)
    end
end

function kerrorsmatch{T1 <: Any, T2 <: Any}(a::T1, b::T2, errors::Int)::Bool
    # if length(a) < length(b)
    #     a, b = b, a
    # end

    alen::Int = length(a)
    blen::Int = length(b)

    alen == 0 && return alen == blen
    blen == 0 && return true

    C::Vector{Int} = Vector{Int}(0:blen)

    for i in 1:alen
	prevA::Int = 0
	prevC::Int = C[1]
	j::Int = 1

	while j <= blen
	    cost::Int = 1
	    @inbounds if a[i] == b[j]
		cost = 0
	    end
	    @inbounds C[j] = prevA
	    j += 1
	    @inbounds prevA = min(C[j]+1, prevA+1, prevC+cost)
	    @inbounds prevC = C[j]
	end
	@inbounds C[j] = prevA
        if prevA <= errors
            return true
        end
    end

    return false
end

function kerrorsmatch(fun, wrap, input, config::TextConfig, key, patlist::Vector{Vector{String}}, maxerrors::Int=1)
    # config.stem = false
    patlist = normpatterns(patlist, config)

    itertweets(input) do tweet
        text = normtext(tweet[key], config)

        matches = 0
        for (i, patterns) in enumerate(patlist)
            for pat in patterns
                if kerrorsmatch(text, pat, maxerrors)
                    matches += 1
                    break
                end
            end
            matches != i && break
        end

        wrap(matches == length(patlist)) && fun(tweet)
    end
end

function substringmatch(fun, wrap, input, config::TextConfig, key, patlist::Vector{Vector{String}})
    # patlist = [[normtext(pat, config) for pat in patterns] for patterns in patlist]
    patlist = normpatterns(patlist, config)

    itertweets(input) do tweet
        text = normtext(tweet[key], config)

        matches = 0
        for (i, patterns) in enumerate(patlist)
            for pat in patterns
                if contains(text, pat)
                    matches += 1
                    break
                end
            end
            matches != i && break
        end

        wrap(matches == length(patlist)) && fun(tweet)
    end
end

function main()
    config = TextConfig()
    patlist = Vector{String}[]
    stopwords = ""
    kind = ""
    key = "text"
    errors = 1
    wrapfun = identity

    for arg in ARGS
        arg = strip(arg)
        if arg[1] == '{'
            _config = JSON.parse(arg)
            for (k, v) in JSON.parse(arg)
                setfield!(config, Symbol(k), v)
            end
        elseif arg[1] == '['
            push!(patlist, JSON.parse(arg))
        elseif startswith(arg, "file:")
            push!(patlist, readlines(arg[6:end]))
        elseif startswith(arg, "stopwords:")
            stopwords = arg[length("stopwords:")+1:end]
        else
            if startswith(arg, "substring")
                kind = "substring"
            elseif startswith(arg, "intersect")
                kind = "intersect"
            elseif startswith(arg, "kerror")
                kind = "kerrors"
                if ':' in kind
                    errors = parse(Int, split(arg, ':')[2])
                end
            elseif startswith(arg, "voc")
                kind = "voc"
            elseif startswith(arg, "key")
                key = split(arg, ':')[2]
            elseif startswith(arg, "neg")
                wrapfun = (x) -> !x
            else
                error("Unknown command $arg")
            end
        end
    end

    load(config, stopwords)

    if kind == "substring"
        substringmatch(wrapfun, STDIN, config, key, patlist) do tweet
            write(STDOUT, JSON.json(tweet), '\n')
        end
    elseif kind == "intersect"
        intersectmatch(wrapfun, STDIN, config, key, patlist) do tweet
            write(STDOUT, JSON.json(tweet), '\n')
        end
    elseif kind == "kerrors"
        kerrorsmatch(wrapfun, STDIN, config, key, patlist, errors) do tweet
            write(STDOUT, JSON.json(tweet), '\n')
        end
    else
        info("""
Usage: julia PatternSearch/matches.jl KIND KEY CONFIG STOPWORDS PAT1 PAT2 ...
    Where

    KIND is substring|intersect|kerrors|voc|radius
    - in the case of kerrors you can specify the maximum number of errors using the syntax kerrors:numerrors
    - the default kind is not specified, it is a mandatory argument

    KEY indicates the keywords containing the text
    - defaults to "text"

    CONFIG is a json-dictionary representing a TextConfig object, it controls how the text is preprocessed for the match
    - please check PatternSearch/textmodel.jl to check the fields and default values
    - the default value is created with the default constructor of TextConfig

    STOPWORDS indicates the path of the stopwords file
    - the syntax is stopwords:filename

    PAT is a json-list of strings, each string is a pattern in the list
    - PatternSearch/matches.jl basically evaluates a intersection-union tree
    - each PAT defines a list of patterns, the text needs to match at least one to be a candidate (union)
    - a real match is a candidate matching all given PATs (intersection)

    The input is always read from STDIN

    Note: This tool is intended to be used with unix pipes and command line tools like awk, perl, grep, PatternSearch/matches.jl, etc.;
    so a lot of possible features will never be implemented.

    """)
        error("please read the help")
    end
end

main()

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
mutable struct LinearRegressor
    _w::AbstractArray{Float64,1}
    _var::Float64
end

function fitting(
    regressor::LinearRegressor,
    Phi::AbstractArray{Float64,2},
    t::AbstractArray{Float64,1},
)
    # Phi is the design matrix of size (n_samples, n_feature)
    # t is the target vector of size (sample_size, 1)
    # `pinv(Phi)` = (Phi^T * Phi)^(-1) * Phi^T
    w = pinv(Phi) * t
    error = Phi * w - t
    var = sum(error .^ 2) / size(t)[1]

    regressor._w = w
    regressor._var = var
end

function predict(
    regressor::LinearRegressor,
    phi::AbstractArray{Float64,1},
    return_std::Bool,
)
    # phi is the transformed vector of size n_features
    y = transpose(phi) * regressor._w
    std = regressor._var

    if return_std == true
        return y, std
    else
        return y
    end
end

function predict(
    regressor::LinearRegressor,
    Phi::AbstractArray{Float64,2},
    return_std::Bool,
)
    # Phi is the design matrix [phi(x_1), phi(x_2), ,,, phi(x_N)]^T of size (n_samples, n_features)
    y = Phi * regressor._w
    std = regressor._var

    if return_std == true
        return y, std
    else
        return y
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# this file is only used to plot the results computed at luanrc

using PyPlot
using DataFrames

# scaled up problem
alg_prop_da = Matrix(readtable("DWPSDE model/analyses/alg prop/alg_prop_da_simdatadt.csv"))
alg_prop_ada = Matrix(readtable("DWPSDE model/analyses/alg prop/alg_prop_ada_simdatadt.csv"))

# small problem
alg_prop_da = Matrix(readtable("DWPSDE model/analyses/alg prop/alg_prop_da_smallsimdatadt.csv"))
alg_prop_ada = Matrix(readtable("DWPSDE model/analyses/alg prop/alg_prop_ada_smallsimdatadt.csv"))

function print_stats(x::Vector)
  @printf "--------------------------\n"
  @printf "Mean:           %.4f\n"  mean(x)
  @printf "Std:            %.4f\n"  std(x)
  describe(x)
end

# text and lable size
text_size = 25
label_size = 20

# run time

println("Runtime:")
print_stats(alg_prop_da[:,1])
print_stats(alg_prop_ada[:,1])

# box plot
run_times = zeros(size(alg_prop_da,1), 2)
run_times[:,1] = alg_prop_da[:,1]
run_times[:,2] = alg_prop_ada[:,1]

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.boxplot(run_times)
PyPlot.xticks([1, 2], ["DA", "ADA"])
ax[:tick_params]("both",labelsize = label_size)

# histogram
PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](alg_prop_da[:,1],10, alpha = 0.6)
PyPlot.plt[:hist](alg_prop_ada[:,1],10, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

speed_up = alg_prop_da[:,1]./alg_prop_ada[:,1]
print_stats(speed_up)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](speed_up,20, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

# nbr pf eval
println("Nbr pf eval:")
print_stats(alg_prop_da[:,2])
print_stats(alg_prop_ada[:,2])

# box plot
pf_eval = zeros(size(alg_prop_da,1), 2)
pf_eval[:,1] = alg_prop_da[:,2]
pf_eval[:,2] = alg_prop_ada[:,2]

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.boxplot(pf_eval)
PyPlot.xticks([1, 2], ["DA", "ADA"])
ax[:tick_params]("both",labelsize = label_size)

# histogram
PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](alg_prop_da[:,2],10, alpha = 0.6)
PyPlot.plt[:hist](alg_prop_ada[:,2],10, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

diff_nbr_pf = alg_prop_da[:,2] - alg_prop_ada[:,2]

print_stats(diff_nbr_pf)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](diff_nbr_pf,10, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

# nbr pf eval in secound stage

println("Nbr pf eval in secound stage:")

print_stats(alg_prop_da[:,3])
print_stats(sum(alg_prop_ada[:,end-3:end],2)[:])

pf_eval = zeros(size(alg_prop_da,1), 2)
pf_eval[:,1] = alg_prop_da[:,3]
pf_eval[:,2] = sum(alg_prop_ada[:,end-3:end],2)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.boxplot(pf_eval)
ax[:tick_params]("both",labelsize = label_size)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](alg_prop_da[:,3],10, alpha = 0.6)
PyPlot.plt[:hist](sum(alg_prop_ada[:,end-3:end],2)[:],10, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

diff_nbr_pf_secound_stage = alg_prop_da[:,3]./sum(alg_prop_ada[:,end-3:end],2)[:]

print_stats(diff_nbr_pf_secound_stage)

PyPlot.figure(figsize=(10,10))
ax = axes()
PyPlot.plt[:hist](diff_nbr_pf_secound_stage,20, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

# nbr ord. mh.

print_stats(alg_prop_da[:,end])
print_stats(alg_prop_ada[:,3])

PyPlot.figure()
ax = axes()
PyPlot.plt[:hist](alg_prop_da[:,end],5, alpha = 0.6)
PyPlot.plt[:hist](alg_prop_ada[:,3],5, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

# nbr times in secound stage

print_stats(alg_prop_da[:,3])
print_stats(sum(alg_prop_ada[:,4:5],2)[:])

PyPlot.figure()
ax = axes()
PyPlot.plt[:hist](alg_prop_da[:,3],5, alpha = 0.6)
PyPlot.plt[:hist](sum(alg_prop_ada[:,4:5],2),5, alpha = 0.6)
ax[:tick_params]("both",labelsize = label_size)

# nbr cases1_3 and cases2_4

print_stats(alg_prop_ada[:,4])
print_stats(alg_prop_ada[:,5])

PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,4],5, alpha = 0.6)

PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,5],5, alpha = 0.6)

# analyses of cases

# nbr cases
print_stats(alg_prop_ada[:,6])
print_stats(alg_prop_ada[:,7])
print_stats(alg_prop_ada[:,8])
print_stats(alg_prop_ada[:,9])

PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,6],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,7],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,8],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,9],5, alpha = 0.6)

# pf eval per case
print_stats(alg_prop_ada[:,10])
print_stats(alg_prop_ada[:,11])
print_stats(alg_prop_ada[:,12])
print_stats(alg_prop_ada[:,13])

PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,10],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,11],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,12],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](alg_prop_ada[:,13],5, alpha = 0.6)

# % pf eval per case

percentage_pf_eval = alg_prop_ada[:,10:13]./alg_prop_ada[:,6:9]

print_stats(percentage_pf_eval[:,1])
print_stats(percentage_pf_eval[:,2])
print_stats(percentage_pf_eval[:,3])
print_stats(percentage_pf_eval[:,4])

PyPlot.figure()
PyPlot.plt[:hist](percentage_pf_eval[:,1],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](percentage_pf_eval[:,2],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](percentage_pf_eval[:,3],5, alpha = 0.6)
PyPlot.figure()
PyPlot.plt[:hist](percentage_pf_eval[:,4],5, alpha = 0.6)

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using Spring

sp=SpringParam(ls=0.1,ms=1,ks=2.5,ns=1000)
ip=InitParam(λ=2.0,shift=0.1,pos=sp.ls*sp.ns/2)
ap=AnimParam(δt=0.2,nδt=5000,nδtperframe=100)
# sp=SpringParam(ls=0.1,ms=1,ks=1,ns=1000)
# ip=InitParam(λ=1.0,shift=1.5,pos=sp.ls*sp.ns/2)
# ap=AnimParam(δt=1,nδt=1000,nδtperframe=10)

animate_spring(sp,ip,ap)