text:
#Tensors in Kruskal (CP) format + functions

export ktensor, randktensor, arrange, arrange!, cp_als, display, extract, fixsigns, fixsigns!, full, innerprod, isequal, krank, kten2TT, minus, mtimes, mttkrp, ncomponents, ndims
export normalize, normalize!, nvecs, permutedims, plus, redistribute, redistribute!, size, tenmat, tocell, ttensor, ttm, ttv, uminus, norm

"""
    ktensor(fmat)
    ktensor(lambda,fmat)

Tensor in Kruskal format defined by its vector of weights lambda and factor matrices. Default lambda: vector of ones.
For ktensor X, X.isorth=true if factor matrices are othonormal.
"""
mutable struct ktensor
	lambda::AbstractVector{<:Number}
	fmat::MatrixCell
	isorth::Bool
	function ktensor(lambda::AbstractVector{<:Number},fmat::MatrixCell,isorth::Bool)
    @assert(length(unique([size(fmat[n],2) for n=1:length(fmat)]))==1,"Factor matrices must have same number of columns.")
    @assert(length(lambda)==size(fmat[1],2),"Dimension mismatch.")
		for A in fmat
			if norm(A'*A-eye(size(A,2)))>(size(A,1)^2)*eps()
				isorth=false
			end
		end
		new(lambda,fmat,isorth)
	end
end
#ktensor(lambda::AbstractArray{T1,1},fmat::MatrixCell,isorth::Bool) where {T1<:Number,T2<:Number}=ktensor{T1,T2}(lambda,fmat,isorth)
ktensor(lambda::AbstractVector{<:Number},fmat::MatrixCell)=ktensor(lambda,fmat,true)
ktensor(lambda::AbstractVector{<:Number},mat::AbstractMatrix{<:Number})=ktensor(lambda,collect(mat),true)
ktensor(fmat::MatrixCell)=ktensor(ones(size(fmat[1],2)),fmat)
# function ktensor(lambda::AbstractVector{<:Number},mat...)
#   fmat=MatrixCell(undef,0)
#   for M in mat
# 			push!(fmat,M)
# 	end
# 	ktensor(lambda,fmat,true)
# end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ktensor(lambda::AbstractVector{<:Number},fmat::Array{Matrix{T},1},isorth::Bool) where {T<:Number}=ktensor(lambda,MatrixCell(fmat),isorth)
ktensor(lambda::AbstractVector{<:Number},fmat::Array{Matrix{T},1}) where {T<:Number}=ktensor(lambda,MatrixCell(fmat),true)
ktensor(fmat::Array{Matrix{T},1}) where {T<:Number}=ktensor(MatrixCell(fmat))
"""
    randktensor(I::Vector,R::Integer)
    randktensor(I::Integer,R::Integer,N::Integer)

Create random Kruskal tensor of size I with R components, or of order N and size I × ⋯ × I.
"""
function randktensor(sz::AbstractVector{<:Integer},R::Integer)
  fmat=[randn(sz[n],R) for n=1:length(sz)] #create random factor matrices
  ktensor(fmat)
end
randktensor(sz::Number,R::Number,N::Integer)=randktensor(repeat([sz],N),R)
#For input defined as tuples or nx1 matrices - ranttensor(([I,I,I],[R,R,R]))

"""
    arrange(X[,n::Integer])
    arrange(X[,P::Vector])

Arrange the rank-1 components of a ktensor: normalize the columns of the factor matrices and then sort the ktensor components by magnitude, greatest to least.

## Arguments:
- `n`: Absorb the weights into the nth factor matrix instead of lambda.
- `P`: Rearrange the components of X according to the permutation P. P should be a permutation of 1 to ncomponents(X).
"""
function arrange(X::ktensor,mode=-1)
    N=ndims(X)
    Xnorm=normalize(X) #Ensure that matrices are normalized
    p=sortperm(Xnorm.lambda,rev=true) #Sort
    lambda=Xnorm.lambda[p]
    fmat=MatrixCell(undef,N)
    [fmat[n]=Xnorm.fmat[n][:,p] for n=1:N]
    #Absorb the weight into one factor, if requested
    if mode>0
        fmat[mode]=fmat[mode].*repeat(lambda',size(X,mode),1) #fmat[mode]*diagm(lambda)
        lambda = fill(1,size(lambda))#ones(lambda)
    end
    ktensor(lambda,fmat)
end
function arrange(X::ktensor,perm::AbstractVector{<:Integer})
    N=ndims(X)
    @assert(collect(1:ncomponents(X))==sort(perm),"Invalid permutation.")
    lambda = X.lambda[perm]
    fmat=MatrixCell(undef,N)
    [fmat[n]=X.fmat[n][:,perm] for n=1:N]
    ktensor(lambda,fmat)
end

"""
    arrange!(X[,n::Integer])
    arrange!(X[,P::Vector])

Arrange the rank-1 components of a ktensor: normalize the columns of the factor matrices and then sort the ktensor components by magnitude, greatest to least.
Rewrite ktensor.

## Arguments:
- `n`: Absorb the weights into the nth factor matrix instead of lambda.
- `P`: Rearrange the components of X according to the permutation P. P should be a permutation of 1 to ncomponents(X).
"""
function arrange!(X::ktensor,mode=-1)
    N=ndims(X)
    normalize!(X) #Ensure that matrices are normalized
    p=sortperm(X.lambda,rev=true) #Sort
    X.lambda=X.lambda[p]
    [X.fmat[n]=X.fmat[n][:,p] for n=1:N]
    #Absorb the weight into one factor, if requested
    if mode>0
        X.fmat[mode]=X.fmat[mode].*repeat(X.lambda',size(X,mode),1) #X.fmat[mode]*diagm(X.lambda)
        X.lambda = fill(1,size(X.lambda))#ones(X.lambda)
    end
    X
end
function arrange!(X::ktensor,perm::AbstractVector{<:Integer})
    N=ndims(X)
    @assert(collect(1:ncomponents(X))==sort(perm),"Invalid permutation.")
    X.lambda = X.lambda[perm]
    [X.fmat[n]=X.fmat[n][:,perm] for n=1:N]
    X
end

#Compute a CP decomposition with R components of a tensor X. **Documentation in tensor.jl.
function cp_als(X::ktensor,R::Integer;init="rand",tol=1e-4,maxit=1000,dimorder=[])
    N=ndims(X)
    nr=norm(X)
    K=ktensor
    if length(dimorder) == 0
        dimorder=collect(1:N)
    end
    fmat=MatrixCell(undef,N)
    if isa(init,Vector) || isa(init,MatrixCell)
        @assert(length(init)==N,"Wrong number of initial matrices.")
        for n in dimorder[2:end]
            @assert(size(init[n])==(size(X,n),R),"$(n)-th initial matrix is of wrong size.")
            fmat[n]=init[n]
        end
    elseif init=="rand"
        [fmat[n]=rand(size(X,n),R) for n in dimorder[2:end]]
    elseif init=="eigs" || init=="nvecs"
        [fmat[n]=nvecs(X,n,R) for n in dimorder[2:end]]
    else
        error("Initialization method wrong.")
    end
    G = zeros(R,R,N); #initalize gramians
    [G[:,:,n]=fmat[n]'*fmat[n] for n in dimorder[2:end]]
    fit=0
    for k=1:maxit
        fitold=fit
        lambda=[]
        for n in dimorder
            fmat[n]=mttkrp(X,fmat,n)
            W=reshape(prod(G[:,:,setdiff(collect(1:N),n)],dims=3),Val(2))
            fmat[n]=fmat[n]/W
            if k == 1
                lambda = sqrt.(sum(fmat[n].^2,dims=1))[:] #2-norm
            else
                lambda = maximum(maximum(abs.(fmat[n]),dims=1),dims=1)[:] #max-norm
            end
            fmat[n] = fmat[n]./lambda'
            G[:,:,n] = fmat[n]'*fmat[n]
        end
        K=ktensor(lambda,fmat)
        if nr==0
            fit=norm(K)^2-2*innerprod(X,K)
        else
            nr_res=sqrt.(nr^2+norm(K)^2-2*innerprod(X,K))
            fir=1-nr_res/nr
        end
        fitchange=abs.(fitold-fit)
        if k>1 && fitchange<tol
            break
        end
    end
    arrange!(K)
    fixsigns!(K)
    K
end

#Display a ktensor. **Documentation in ttensor.jl.
function display(X::ktensor,name="ktensor")
    println("Kruskal tensor of size ",size(X),":\n")
    println("$name.lambda: ")
	flush(stdout)
	display(X.lambda)
    for n=1:ndims(X)
        println("\n\n$name.fmat[$n]:")
		flush(stdout)
		display(X.fmat[n])
    end
end

"""
    extract(X,factors)

Create a new ktensor with only the specified factors.
"""
function extract(X::ktensor,factors::AbstractVector{<:Integer})
    N=ndims(X)
    lambda = X.lambda[factors]
    fmat = MatrixCell(undef,N)
    [fmat[n]=X.fmat[n][:,factors] for n=1:N]
    ktensor(lambda, fmat)
end
extract(X::ktensor,factors::Integer)=extract(X,[factors])

"""
    fixsign(X)

Fix sign ambiguity of a ktensor.
"""
function fixsigns(X::ktensor)
    N=ndims(X)
    sgn=zeros(Int,N)
    fmat=deepcopy(X.fmat)
     for r=1:ncomponents(X)
        for n=1:N
            maxind=argmax(abs.(X.fmat[n][:,r]))
            sgn[n]=sign(X.fmat[n][maxind,r])
        end
        negind=findall(sgn.==-1)
        nflip=2*floor(Int,length(negind)/2)
        for i=1:nflip
            fmat[negind[i]][:,r]=-X.fmat[negind[i]][:,r]
        end
    end
    ktensor(X.lambda,fmat)
end
"""
    fixsign(X)

Fix sign ambiguity of a ktensor. Rewrite ktensor.
"""
function fixsigns!(X::ktensor)
    N=ndims(X)
    sgn=zeros(Int,N)
     for r=1:ncomponents(X)
        for n=1:N
            maxind=argmax(abs.(X.fmat[n][:,r]))
            sgn[n]=sign(X.fmat[n][maxind,r])
        end
        negind=findall(sgn.==-1)
        nflip=2*floor(Int,length(negind)/2)
        for i=1:nflip
            X.fmat[negind[i]][:,r]=-X.fmat[negind[i]][:,r]
        end
    end
end

#Makes full tensor out of a ktensor. **Documentation in ttensor.jl.
function full(X::ktensor)
  reshape(Array(X.lambda'*khatrirao(X.fmat[end:-1:1])'),size(X)...)
end

#Inner product of two ktensors. **Documentation in tensor.jl.
function innerprod(X1::ktensor,X2::ktensor)
  @assert(size(X1)==size(X2),"Dimension mismatch.")
  inpr=X1.lambda*X2.lambda'
  for n=1:ndims(X1)
    inpr=inpr.*(X1.fmat[n]'*X2.fmat[n])
  end
  sum(inpr[:])
end
function innerprod(X1::ktensor,X2)
  @assert(isa(X2,ttensor) || isa(X2,Array), "Inner product not available for type $(typeof(X2)).")
  @assert(size(X1)==size(X2),"Dimension mismatch.")
  inpr=0
  for r=1:ncomponents(X1)
    inpr+=X1.lambda[r]*ttv(X2,[X1.fmat[n][:,r] for n=1:ndims(X1)])[1]
  end
  inpr
end
innerprod(X1,X2::ktensor)=innerprod(X2,X1)

#Two ktensors are equal if they have equal components. **Documentation in ttensor.jl.
function isequal(X1::ktensor,X2::ktensor)
  if (X1.lambda == X2.lambda) && (X1.fmat == X2.fmat)
    true
  else
    false
  end
end
==(X1::ktensor,X2::ktensor)=isequal(X1,X2)

"""
    krank(X::ktensor)

Represenation k-rank. Equal to number of columns of factor matrices if ktensor X.
"""
function krank(X::ktensor)
    size(X.fmat[1],2)
end

"""
    kten2TT(X::ktensor)

Transform ktensor to TTtensor.
"""
function kten2TT(X::ktensor)
    N=ndims(X)
    Isz=size(X)
    r=krank(X)
    G=CoreCell(undef,N)
    G[1]=reshape(X.fmat[1],(1,Isz[1],r))
    G[N]=reshape(X.fmat[N]',(r,Isz[N],1))
    for n=2:N-1
        G[n]=zeros(r,Isz[n],r)
        for i=1:Isz[n]
            G[n][:,i,:]=Diagonal(X.fmat[n][i,:])
        end
    end
    TTtensor(G)
end

#Subtraction of two ktensors. **Documentation in ttensor.jl.
function minus(X1::ktensor,X2::ktensor)
  1*X1+(-1)*X2
end
-(X1::ktensor,X2::ktensor)= minus(X1,X2)

#Scalar times tensor. **Documentation in ttensor.jl.
function mtimes(α::Number,X::ktensor)
	ktensor(α*X.lambda,X.fmat);
end
*(α::Number,X::ktensor)=mtimes(α,X)
*(X::ktensor,α::Number)=*(α,X)

#Matricized ktensor times Khatri-Rao product. **Documentation in tensor.jl.
function mttkrp(X::ktensor,M::MatrixCell,n::Integer)
  N=ndims(X)
  @assert(length(M) == N,"Wrong number of matrices")
  n==1 ? R=size(M[2],2) : R=size(M[1],2)
  modes=setdiff(1:N,n)
  W=repeat(X.lambda,1,R) #matrix of weights
  for m in modes
    W=W.*(X.fmat[m]'*M[m])
  end
  X.fmat[n]*W
end
mttkrp(X::ktensor,M::Array{Matrix{T},1},n::Integer) where {T<:Number}=mttkrp(X,MatrixCell(M),n)

"""
    ncomponents(X)

Number of components of a ktensor.
"""
function ncomponents(X::ktensor)
  length(X.lambda)
end

#Number of modes of a ttensor. **Documentation in Base.
function ndims(X::ktensor)
  length(X.fmat)
end

"""
    normalize(X[,n;normtype,factor])

Normalize columns of factor matrices of a ktensor. Also ensures that lambda is positive.

## Arguments:
- `n`: Absorbe the excess weight into nth factor matrix.
   If n=0, equally divide the weights across the factor matrices and set lambda to vector of ones.
   If n="sort", sort the components by lambda value, from greatest to least.
- `normtype`: Default: 2.
- `factor`: Just normalize specified factor.
"""
function normalize(X::ktensor,mode=-1;normtype=2,factor=-1)
  if mode=="sort"
    mode=-2
  end
  N=ndims(X)
  lambda=deepcopy(float(X.lambda))
  fmat=deepcopy(X.fmat)
  if factor>0
    @assert(0<factor<=N,"Factor exceeds dimension.")
    for r = 1:ncomponents(X)
      nr = norm(fmat[factor][:,r],normtype);
      lambda[r]=lambda[r]*nr
      if nr > 0
         fmat[factor][:,r]=fmat[factor][:,r]/nr
      end
    end
  end
  for r=1:ncomponents(X)
    for n=1:N
      nr=norm(fmat[n][:,r],normtype)
      lambda[r]=lambda[r]*nr
      if nr>0
        fmat[n][:,r]=fmat[n][:,r]/nr
      end
    end
  end
  #Check that all the lambda values are positive
  negind = findall(lambda .< 0)
  for i in negind
    lambda[i] = -1 * lambda[i]
    fmat[1][:,i] = -1 * fmat[1][:,i]
  end
  if mode == 0
    d=(lambda').^(1/N)
   [fmat[n]=fmat[n].*repeat(d,size(X,n),1) for n=1:N] #fmat[n]*diagm(d)
    lambda = fill(1,size(lambda))#ones(lambda)
  elseif mode > 0
    fmat[mode] = fmat[mode].*repeat(lambda',size(X,mode),1) #fmat[mode]*diagm(lambda)
    lambda = fill(1,size(lambda))#ones(lambda)
  elseif mode==-2
    p = sortperm(lambda,rev=true)
    return arrange(ktensor(lambda,fmat),p)
  end
  ktensor(lambda,fmat)
end

"""
    normalize!(X[,n,normtype,factor])

Normalize columns of factor matrices of a ktensor. Rewrite ktensor. Also ensures that lambda is positive.

## Arguments:
- `n`: Absorbe the excess weight into nth factor matrix.
   If n=0, equally divide the weights across the factor matrices and set lambda to vector of ones.
   If n="sort", sort the components by lambda value, from greatest to least.
- `normtype`: Default: 2.
- `factor`: Just normalize specified factor.
"""
function normalize!(X::ktensor,mode=-1;normtype=2,factor=-1)
  if mode=="sort"
    mode=-2
  end
  N=ndims(X)
  X.lambda=float(X.lambda)
  if factor>0
    @assert(0<factor<=N,"Factor exceeds dimension.")
    for r = 1:ncomponents(X)
      nr = norm(X.fmat[factor][:,r],normtype);
      X.lambda[r]=X.lambda[r]*nr
      if nr > 0
         X.fmat[factor][:,r]=X.fmat[factor][:,r]/nr
      end
    end
  end
  for r=1:ncomponents(X)
    for n=1:N
      nr=norm(X.fmat[n][:,r],normtype)
      X.lambda[r]=X.lambda[r]*nr
      if nr>0
        X.fmat[n][:,r]=X.fmat[n][:,r]/nr
      end
    end
  end
  #Check that all the lambda values are positive
  negind = findall(X.lambda .< 0)
  for i in negind
    X.fmat[1][:,i] = -1 * X.fmat[1][:,i]
    X.lambda[i] = -1 * X.lambda[i]
  end
  if mode == 0
     d=(X.lambda').^(1/N)
     [X.fmat[n]=X.fmat[n].*repeat(d,size(X,n),1) for n=1:N] #X.fmat[n]*diagm(d)
     X.lambda = fill(1,size(X.lambda))#ones(X.lambda)
  elseif mode > 0
    X.fmat[mode] = X.fmat[mode].*repeat(X.lambda',size(X,mode),1) #X.fmat[n]*diagm(X.lambda)
    X.lambda = fill(1,size(X.lambda))#ones(X.lambda)
  elseif mode==-2
    p = sortperm(X.lambda,rev=true)
    X = arrange(X,p)
  end
   X
end

#Computes the r leading singular vectors of mode-n matricization of a tensor X.
# **Documentation in tensor.jl.
function nvecs(X::ktensor,n::Integer,r=0;flipsign=false)
    M = X.lambda * X.lambda'
    for m in setdiff(1:ndims(X),n)
        M = M .* (X.fmat[m]' * X.fmat[m])
    end
    #V=eigs(Symmetric(X.fmat[n] * M * X.fmat[n]'),nev=r,which=:LM)[2] #has bugs
    V=eigen(Symmetric(X.fmat[n] * M * X.fmat[n]')).vectors[:,end:-1:end-r+1]
    if flipsign
        maxind = findmax(abs.(V),1)[2]
        for i = 1:r
            ind=ind2sub(size(V),maxind[i])
            if V[ind...] < 0
                V[:,ind[2]] = V[:,ind[2]] * -1
            end
        end
    end
    V
end

#Permute dimensions of a ktensor. **Documentation in Base.
function permutedims(X::ktensor,perm::AbstractVector{<:Integer})
    N=ndims(X)
    @assert(collect(1:N)==sort(perm),"Invalid permutation.")
    ktensor(X.lambda,X.fmat[perm])
end

#Addition of two ktensors. **Documentation in ttensor.jl.
function plus(X1::ktensor,X2::ktensor)
	@assert(size(X1) == size(X2),"Dimension mismatch.")
  	lambda=[X1.lambda; X2.lambda]
	fmat=[[X1.fmat[n] X2.fmat[n]] for n=1:ndims(X1)] #concatenate factor matrices
	ktensor(lambda,fmat)
end
+(X1::ktensor,X2::ktensor)=plus(X1,X2)

"""
    redistribute(X,n)

Distribute lambda values to a mode n.
"""
function redistribute(X::ktensor,mode::Integer)
    lambda=deepcopy(float(X.lambda))
    fmat=deepcopy(X.fmat)
    for r = 1:ncomponents(X)
        fmat[mode][:,r] = fmat[mode][:,r] * lambda[r];
        lambda[r] = 1;
    end
    ktensor(lambda,fmat)
end

"""
    redistribute!(X,n)

Distribute lambda values to a mode n. Rewrite ktensor X.
"""
function redistribute!(X::ktensor,mode::Integer)
    for r = 1:ncomponents(X)
        X.fmat[mode][:,r] = X.fmat[mode][:,r] * X.lambda[r];
        X.lambda[r] = 1;
    end
    X
end

#Size of a ktensor. **Documentation in Base.
function size(X::ktensor)
  tuple([size(X.fmat[n],1) for n=1:length(X.fmat)]...)
end
#Size of n-th mode of a ktensor.
function size(X::ktensor,n::Integer)
  size(X.fmat[n],1)
end

function Base.show(io::IO,X::ktensor)
    display(X)
end

#Mode-n matricization of a ktensor. **Documentation in tensor.jl.
tenmat(X::ktensor,n::Integer)=tenmat(full(X),n)

"""
tocell(X[,n])

Converts ktensor X into MatrixCell. If n specified, absorb the weights into the nth factor matrix.
"""
function tocell(X::ktensor,mode::Integer)
    Y = normalize(X,mode)
    Y.fmat
end
function tocell(X::ktensor)
    if isequal(X.lambda,ones(size(X.lambda)))
        return X.fmat;
    end
    N=ndims(X)
    d = (X.lambda').^(1/N)
    fmat=MatrixCell(undef,N)
    for n = 1:N
        fmat[n] = X.fmat[n].*repeat(d,size(X.fmat[n],1),1)
    end
    fmat
end

"""
---
    ttensor(X::ktensor)

Create ttensor out of a ktensor.
"""
function ttensor(X::ktensor)
    Xred=redistribute(X,1)
    ttensor(neye(ncomponents(X),order=ndims(X)),Xred.fmat)
end

#Mode-n multiplication of ktensor and matrix. **Documentation in tensor.jl.
#t='t' transposes matrices
function ttm(X::ktensor,M::MatrixCell,modes::AbstractVector{<:Integer},t='n')
  if t=='t'
	 [M[n]=M[n]' for n=1:length(M)]
	end
	@assert(length(modes)<=length(M),"Too few matrices")
	@assert(length(M)<=ndims(X),"Too many matrices")
  if length(modes)<length(M)
    M=M[modes]
  end
  fmat=deepcopy(X.fmat)
  for n=1:length(modes)
    fmat[modes[n]]=M[n]*fmat[modes[n]]
  end
    ktensor(X.lambda,fmat)
end
ttm(X::ktensor,M::MatrixCell,modes::AbstractRange{<:Integer},t='n')=ttm(X,M,collect(modes),t)
ttm(X::ktensor,M::Matrix{<:Number},n::Integer,t='n')=ttm(X,MatrixCell([M]),[n],t)
ttm(X::ktensor,M::MatrixCell,t::Char)=ttm(X,M,1:length(M),t)
ttm(X::ktensor,M::MatrixCell)=ttm(X,M,1:length(M))
function ttm(X::ktensor,M::MatrixCell,n::Integer,t='n')
	if n>0
		ttm(X,M[n],n,t)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,M,modes,t)
	end
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ttm(X::ktensor,M::Array{Matrix{T},1},modes::AbstractVector{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::ktensor,M::Array{Matrix{T},1},modes::AbstractRange{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::ktensor,M::Array{Matrix{T},1},t::Char) where {T<:Number}=ttm(X,MatrixCell(M),t)
ttm(X::ktensor,M::Array{Matrix{T},1}) where {T<:Number}=ttm(X,MatrixCell(M))
ttm(X::ktensor,M::Array{Matrix{T},1},n::Integer,t='n') where {T<:Number}=ttm(X,MatrixCell(M),n,t)

#Mode-n multiplication ktensor and vector. **Documentation in tensor.jl.
#t='t' transposes matrices
function ttv(X::ktensor,V::VectorCell,modes::AbstractVector{<:Integer})
  N=ndims(X)
  remmodes=setdiff(1:N,modes)
  lambda=deepcopy(float(X.lambda))
  if length(modes) < length(V)
    V=V[modes]
  end
  for n=1:length(modes)
      lambda=lambda.*(X.fmat[modes[n]]'*V[n])
  end
  if isempty(remmodes)
    sum(lambda)
  else
    ktensor(lambda,X.fmat[remmodes])
  end
end
ttv(X::ktensor,v::AbstractVector{<:Number},n::Integer)=ttv(X,VectorCell([v]),[n])
ttv(X::ktensor,V::VectorCell,modes::AbstractRange{<:Integer})=ttv(X,V,collect(modes))
ttv(X::ktensor,V::VectorCell)=ttv(X,V,1:length(V))
function ttv(X::ktensor,V::VectorCell,n::Integer)
	if n>0
		ttm(X,V[n],n)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,A,modes)
	end
end
#If array of vectors isn't defined as VectorCell, but as V=[v1,v2,...,vn]:
ttv(X::ktensor,V::Array{Vector{T},1},modes::AbstractVector{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::ktensor,V::Array{Vector{T},1},modes::AbstractRange{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::ktensor,V::Array{Vector{T},1}) where {T<:Number}=ttv(X,VectorCell(V))
ttv(X::ktensor,V::Array{Vector{T},1},n::Integer) where {T<:Number}=ttv(X,VectorCell(V),n)

#**Documentation in ttensor.jl.
uminus(X::ktensor)=mtimes(-1,X)
-(X::ktensor)=uminus(X)

#Frobenius norm of a ktensor. **Documentation in Base.
function norm(X::ktensor)
  nrm=X.lambda*X.lambda'
  for n=1:ndims(X)
    nrm=nrm.*(X.fmat[n]'*X.fmat[n])
  end
  sqrt.(abs.(sum(nrm[:])))
end


text:
"""
    hessian!{T <: AbstractFloat}(
        output::AbstractArray,
        f::Function,
        x::AbstractArray{T},
        buffer::AbstractArray{T},
    )

# Description

Evaluate the Hessian of `f` at `x` using finite differences. Store the results
into `output`. Work with a user-provided `buffer` to ensure that we can work
without copies, but also without mutating `x`. In Rust jargon, we take
ownership of both `output` and `buffer`, but do not require ownership of `x`.
We just need read-only access to `x`.

# Arguments

* `output::AbstractArray`: An array that will be mutated to contain the
    gradient.
* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T}`: The value of `x` at which to evaluate the Hessian of
    `f`.
* `buffer::AbstractArray{T}`: A buffer that is equivalent to `similar(x)`. Used
    as a temporary copy of `x` that can be mutated.

# Returns

* `nothing::Void`: This function is called for its side effects.

# Examples

```julia
import FiniteDiff: hessian!
x = [0.0, 0.0]
output, buffer = similar(x, 2, 2), similar(x, 2)
hessian!(output, x -> sin(x[1]) + 2 * sin(x[2]), x, buffer)
```
"""
function hessian!{T <: AbstractFloat}(
    output::AbstractArray,
    f::Function,
    x::AbstractArray{T},
    buffer::AbstractArray{T},
)::Void
    # Validate that all inputs have the expected number of dimensions.
    n = length(x)
    if size(output) != (n, n) || length(buffer) != n
        throw(DomainError())
    end

    # Cache the value of f(x) to avoid multiple calls to f, which is pure.
    f_x = f(x)

    # Copy x into the user-supplied buffer.
    copy!(buffer, x)

    # Iterate over the rows of the input.
    for i = 1:n
        # Make a copy of the true value of x[i].
        x_i = x[i]

        # Temporarily change buffer[i] to a new value and evaluate f at the
        # new x stored in the buffer.
        ϵ = step_size(HessianMode(), x_i)
        buffer[i] = x_i + ϵ
        f_xpp = f(buffer)

        # Temporarily change buffer[i] to a new value and evaluate f at the
        # new x stored in the buffer.
        buffer[i] = x_i - ϵ
        f_xmm = f(buffer)

        # Store the diagonal component of the Hessian into the output array.
        output[i, i] = (f_xpp - 2 * f_x + f_xmm) / (ϵ * ϵ)

        # Determine the step-size to use for this row of the Hessian on the
        # off-diagonal elements.
        ϵ_i = step_size(CentralMode(), x_i)

        # Iterate over the columns of the input for the current row.
        for j = (i + 1):n
            # Make a copy of the true value of x[i].
            x_j = x[j]

            # Determine the step-size to use for this entry of the Hessian.
            ϵ_j = step_size(CentralMode(), x_j)

            # Temporarily change buffer[i] and buffer[j] to new values and
            # evaluate f at the new values of x stored in the buffer.
            buffer[i] = x_i + ϵ_i
            buffer[j] = x_j + ϵ_j
            f_xpp = f(buffer)
            buffer[j] = x_j - ϵ_j
            f_xpm = f(buffer)

            # Temporarily change buffer[i] and buffer[j] to new values and
            # evaluate f at the new values of x stored in the buffer.
            buffer[i] = x_i - ϵ_i
            buffer[j] = x_j + ϵ_j
            f_xmp = f(buffer)
            buffer[j] = x_j - ϵ_j
            f_xmm = f(buffer)

            # Store the i, j-th component of the Hessian into the output array.
            t4 = convert(T, 4)
            output[i, j] = (f_xpp - f_xpm - f_xmp + f_xmm) / (t4 * ϵ_i * ϵ_j)

            # Restore the true value of x[j] in the buffer.
            buffer[j] = x_j
        end

        # Restore the true value of x[i] in the buffer.
        buffer[i] = x_i
    end

    # Copy the upper triangular values into the lower triangular values to
    # symmetrize the results.
    Base.LinAlg.copytri!(output, 'U')

    # Return nothing. The result is already stored in output.
    return
end

"""
    hessian{T <: AbstractFloat}(f::Function, x::AbstractArray{T})

# Description

Evaluate the Hessian of `f` at `x` using finite differences. See
`hessian!(f, x, buffer)` for more details.

# Arguments

* `f::Function`: The function to be differentiated.
* `x::AbstractArray{T <: AbstractArray}`: The value of `x` at which to evaluate
    the Hessian of `f`.

# Returns

* `output::Array{T <: AbstractArray}`: The Hessian of `f` at `x`.

# Examples

```julia
import FiniteDiff: hessian
x = [0.0, 0.0]
H = hessian(x -> sin(x[1]) + 2 * sin(x[2]), x)
```
"""
function hessian{T <: AbstractFloat}(
    f::Function,
    x::AbstractArray{T},
)::AbstractArray{T}
    # Determine the number of dimensions of the input.
    n = length(x)

    # Allocate memory for both the output and the temporary buffer.
    output, buffer = similar(x, n, n), similar(x, n)

    # Compute the Hessian using the mutating variant of this function.
    hessian!(output, f, x, buffer)

    # Return the output, but discard the temporary buffer.
    return output
end

"""
    hessian(f::Function; mutates::Bool = false)

# Description

Construct a new function that will evaluate the Hessian of `f` at any value
of `x`. The user can whether to return a mutating or non-mutating function
using a keyword argument.

# Arguments

* `f::Function`: The function to be differentiated.

# Keyword Arguments

* `mutates::Bool = false`: Determine whether the resulting function will mutate
    its inputs or will be a pure function.

# Returns

* `h::Function`: A function that will evaluate the Hessian of `f` at any `x`.

# Examples

```julia
import FiniteDiff: hessian
h = hessian(x -> sin(x[1]) + 2 * sin(x[2]), mutates = false)
x = [0.0, 0.0]
h(x)
```
"""
function hessian(f::Function; mutates::Bool = false)::Function
    if mutates
        return (output, x, buffer) -> hessian!(output, f, x, buffer)
    else
        return (x, ) -> hessian(f, x)
    end
end

ext: jl


text:
function numberstotex(d)
    join(map(n -> begin
                 if divrem(n, 1)[2] == 0
                     string(Int(n))
                 else
                     string(n)
                 end
             end, d'[:]), " & ")
end
function datatotextable(raw, rearranged)
    s = "\\begin{tabular}{l|c|cccccc||cccccc}\n"
    s *= "\\multicolumn{2}{c|}{} & \\multicolumn{6}{c||}{\\cellcolor{gray!20} Raw Data} & \\multicolumn{6}{c}{\\cellcolor{gray!20}Sorted Data}\\\\\n"
    s *= " \\multicolumn{2}{c|}{} & \\multicolumn{2}{c}{A} &
    \\multicolumn{2}{c}{B} & \\multicolumn{2}{c||}{C} &
    \\multicolumn{2}{c}{K\$_1\$}  & \\multicolumn{2}{c}{K\$_2\$} &
    \\multicolumn{2}{c}{K\$_3\$} \\\\"
    s *= "Bird & Sequence & P & M & P & M & P & M & f\$_1\$ & f\$_2\$  & f\$_1\$  & f\$_2\$  & f\$_1\$  & f\$_2\$
          \\\\\\midrule\n"
    for (k, (c, f, d)) in raw
        s *= "$k & $(join(map(string, f)))-$(join(map(string, c))) & $(numberstotex(d)) & $(numberstotex(rearranged[k].cached)) \\\\\n"
     end
     s *= "\\end{tabular}"
end
function datatotextable2(raw, rearranged)
    s = "\\begin{tabular}{l|c|ccc||ccc}\n"
    s *= "\\multicolumn{2}{c|}{} & \\multicolumn{3}{c||}{\\cellcolor{gray!20} Raw Data} & \\multicolumn{3}{c}{\\cellcolor{gray!20}Sorted Data}\\\\\n"
    s *= "Bird & Sequence & A & B & C & K\$_1\$  & K\$_2\$ & K\$_3\$ \\\\\\midrule\n"
    for (k, (c, f, d)) in raw
        s *= "$k & $(join(map(string, f)))-$(join(map(string, c))) & $(numberstotex(d)) & $(numberstotex(rearranged[k].cached)) \\\\\n"
     end
     s *= "\\end{tabular}"
end

function plotindividualdata(data, title;
                            isaxis = true, isticks = true,
                            y_axis_line_style = "{<->}",
                            ylabel = "f\$_2\$\\hspace{1.2cm} f\$_1\$",
                            xticklabels = "")
    if length(data) == 3
        append!(data, zeros(3)) # this is a hack to deal with exp 2
    end
    @pgf p = Plot(Coordinates([1:3; 1:3], [data[1:3]; -data[4:end]]))
    if isaxis
        a1 = @pgf Axis({ybar, bar_width = "10pt", title = title,
                        width = "3.5cm", ymin = - maximum(data[4:end]) - 1,
                        ymax = maximum(data[1:3]) + 1,
                        axis_x_line = "middle", axis_y_line = "left",
                        xmin = .5, xmax = 3.5, xtick = [1, 2, 3],
                        xticklabels = xticklabels,
                        xticklabel_shift = "2mm",
                        y_axis_line_style = y_axis_line_style,
                        x_axis_line_style = "{-}",
                        ytick = "\\empty",
                        ylabel = ylabel,
                        ylabel_style = {yshift = "0pt"},
                        yticklabels = ["f\$_2\$", "f\$_1\$"],
                        yticklabel = "{\\scriptsize
                                   \\ifdim\\tick pt < 0pt
                                   \\pgfmathparse{-1*\\tick}
                                   \\pgfmathprintnumber{\\pgfmathresult}
                                   \\else
                                   \\pgfmathprintnumber{\\tick}
                                   \\fi}" }, p)
        if !isticks
            a2["xtick"] = "\\empty"
        end
        TikzPicture(a1)
    else
        p
    end
end

function plotallindividualdata(data; exp = 1, dir, xlabelat = "",
                               titles = Dict(), kwargs...)
    for (k, (c, f, d)) in data
        td = TikzDocument(plotindividualdata(d, haskey(titles, k) ? titles[k] : string(k); xticklabels =
                                             (string(k) == xlabelat ? ["K\$_1\$", "K\$_2\$", "K\$_3\$"] : ""),
                                             kwargs...))
        PGFPlotsX.savetex(joinpath(dir, "data$(replace(string(k), " " => ""))-$exp.tex"), td,
                          include_preamble = false)
    end
end

function plotcomparison(y, models; title = "", ymax = 1, legend = true,
                            legendentries = ["bird independent", "bird dependent",
                                            "food independent", "food dependent"],
                            legendtoname = "modelcomparisonlegend",
                            altpattern = "crosshatch",
                            colors = ["red!70", "blue!70"],
                            exp2 = false
                           )
    plots = Any[]
    coords = [0, 1, 1.6, 2.2, 2.8]
    idx1 = findall(x -> x.bird_indep && x.food_indep, models)
    push!(plots, @pgf Plot({fill = "gray", draw = "gray", bar_shift = exp2 ? "16pt" : "0pt"},
                           Coordinates(coords[1:1], y[1:1])))
    push!(plots, @pgf Plot({fill = colors[2], draw = colors[2]},
                           Coordinates(coords[(2:length(idx1) + 1)], y[idx1])))
    idx2 = findall(x -> !x.bird_indep && x.food_indep, models)[2:end]
    push!(plots, @pgf Plot({fill = colors[1], draw = colors[1]},
                           Coordinates(coords[(2:length(idx2) + 1)], y[idx2])))
    if !exp2
        idx3 = findall(x -> x.bird_indep && !x.food_indep, models)
        push!(plots, @pgf Plot({pattern = altpattern, draw = colors[2],
                                pattern_color = colors[2]},
                               Coordinates(coords[1:length(idx3)], y[idx3])))
        idx4 = findall(x -> !x.bird_indep && !x.food_indep, models)
        push!(plots, @pgf Plot({pattern = altpattern, draw = colors[1],
                                pattern_color = colors[1]},
                               Coordinates(coords[1:length(idx3)], y[idx4])))
    end
    a = @pgf Axis({ybar, bar_width = "14pt", xmin = -.3, xmax = exp2 ? 3.1 : 3.5, width = "\\textwidth",
                   height = "7cm", line_width = "0pt", xtick = coords, ymax = ymax,
                   major_tick_length = 0,
                   xticklabels = [(exp2 ? "\\hspace{32pt}" : "\\hspace{32pt}") *"\\parbox{2cm}{\\centering compartement\\\\ independent}",
                                  (exp2 ? "\\hspace{16pt}" : "\\hspace{16pt}") * "\\parbox{2cm}{\\centering compartement\\\\dependent}",
                                  (exp2 ? "\\hspace{16pt}" : "\\hspace{48pt}") * " CCH",
                                  (exp2 ? "\\hspace{16pt}" : "\\hspace{48pt}") * " FPH 1",
                                  (exp2 ? "\\hspace{16pt}" : "\\hspace{48pt}") * " FPH 2"],
                          ymode = "log", log_origin = "infty",
                          title = title}, plots...)
    if legend
        prepend!(a.contents, ["\\addlegendimage{color = $(colors[2])}
                              \\addlegendentry[color = black]{$(legendentries[1])}
                         \\addlegendimage{color = $(colors[1])}
                         \\addlegendentry[color = black]{$(legendentries[2]) }" *
                         (exp2 ? "" : "\\addlegendimage{fill = black, line width = 0pt}
                         \\addlegendentry[color = black]{$(legendentries[3]) },
                         \\addlegendimage{pattern = $altpattern, line width = 0pt}
                         \\addlegendentry[color = black]{$(legendentries[4]) }")])
        a["legend_to_name"] = legendtoname
        a["legend_columns"] = "-1"
        a["legend_style"] = "draw = none, nodes = {inner sep = 8pt}"
    end
    a
end


text:
"""
    BilateralStructure

An anatomical structure that occurs bilaterally.
"""
struct BilateralStructure{S<:AnatomicalStructure} <: AbstractAnatomicalStructure
    structure::S
end

BilateralStructure(s::Symbol) = BilateralStructure(AnatomicalStructure(s))

BilateralStructure(s::String) = BilateralStructure(Symbol(s))

function BilateralStructure(; name, children::Tuple=(), traits::Tuple=())
    BilateralStructure(AnatomicalStructure(name, children, traits))
end

structure(s::BilateralStructure) = getproperty(s, :structure)
Base.names(s::BilateralStructure) = names(structure(s))

isbilateral(s::AnatomicalStructure) = false
isbilateral(s::BilateralStructure) = true
isbilateral(s::AbstractAnatomicalStructure) = isbilateral(structure(s))

"""
    LeftStructure

Represents the left side of a `BilateralStructure`.
"""
struct LeftStructure{S<:BilateralStructure} <: AbstractAnatomicalStructure
    structure::S

    LeftStructure(s::BilateralStructure) = new{typeof(s)}(s)
end

left(s::BilateralStructure) = LeftStructure(s)

structure(s::LeftStructure) = getproperty(s, :structure)

substructuretype(::S) where {S<:LeftStructure} = substructuretype(S)
substructuretype(::Type{<:LeftStructure{S}}) where S = S

Base.names(s::LeftStructure) = names(structure(s))

Base.show(io::IO, v::LeftStructure) = print(io, Symbol("left $(names(v))"))

"""
    RightStructure

Represents the right side of a `BilateralStructure`.
"""
struct RightStructure{S<:BilateralStructure} <: AbstractAnatomicalStructure
    structure::S

    RightStructure(s::BilateralStructure) = new{typeof(s)}(s)
end

right(s::BilateralStructure) = RightStructure(s)

structure(s::RightStructure) = getproperty(s, :structure)

Base.names(s::RightStructure) = names(structure(s))

Base.show(io::IO, v::RightStructure) = print(io, Symbol("right $(names(v))"))


text:
macro dowhile(condition, block)
    quote
        while true
            $(esc(block))
            if !$(esc(condition))
                break
            end
        end
    end
end

macro dountil(condition, block)
    quote
        while true
            $(esc(block))
            if $(esc(condition))
                break
            end
        end
    end
end

using Primes

arr = [7, 31]

@dowhile (!isprime(arr[1]) && !isprime(arr[2])) begin
    println(arr)
    arr .+= 1
end
println("Done.")

@dountil (isprime(arr[1]) || isprime(arr[2])) begin
    println(arr)
    arr .+= 1
end
println("Done.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module MultipleDispatchKoans

#=
Links:
    - https://docs.julialang.org/en/v1/manual/methods/
=#

struct Dog
    name::String
end

struct Cat
    name::String
end

function hello()
    "hello stranger"
end

function hello(name)
    "hello $name"
end

function pet_talk(pet::Dog)
    "$(pet.name): bark"
end

function pet_talk(pet::Cat)
    "$(pet.name): meow"
end

end


text:
using JetPack, Jets, Test

@testset "jop_permute" begin
    A = JopPermute(JetSpace(Float64,4,2), (1,), ([3;1;2;4],))
    m = [11. 12. ; 21. 22. ; 31. 32. ; 41. 42.]
    d = A*m
    @test d ≈ [31. 32. ; 11. 12. ; 21. 22. ; 41. 42.]

    A = JopPermute(JetSpace(Float64,3,2), (1,2), ([3;2;1],[2;1]))
    m = [11. 12. ; 21. 22. ; 31. 32.]
    d = A*m
    @test d ≈ [32. 31. ; 22. 21. ; 12. 11.]

    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A,m,d)
    @test lhs ≈ rhs

    A = JopPermute(JetSpace(Float64,4), (1,), ([3;2;1;4],))
    m = [1. ; 2. ; 3. ; 4.]
    d = A*m
    @test d ≈ [3. ; 2. ; 1. ; 4.]

    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A,m,d)
    @test lhs ≈ rhs
end


text:
macro func_loader(class, methods...)
    blocks = Expr(:block)
    for method in methods
        block = quote
            function $(esc(method))(args...; kwargs...)
                $(class).$(method)(args...; kwargs...)
            end
        end
        push!(blocks.args, block)
    end
    blocks
end


text:
include("datafitting.jl")
function f(x)
    t = collect(0:3)
    y = [2.0; 0.7; 0.3; 0.1]
    r = y - x[1].*exp.(x[2].*t)
    f =  sum(r.^2)/2
    jac = [-exp.(x[2].*t) -x[1].*t.*exp.(x[2].*t)]
    g = jac'*r
    return f, jac, g, r
end

# good initial guess
x0 = [1;0]
x = gaussnewton(x0, f)
println(x)

# bad initial guess
x0 = [-3;3]
x = gaussnewton(x0, f)
println(x)

# Levenberg-Marquardt with good initial guess
x0 = [1;0]
x = levenbergmarquardt(x0, f)
println(x)

# Levenberg-Marquardt with bad initial guess
x0 = [-3;3]
x = levenbergmarquardt(x0, f)
println(x)


text:
using Documenter, DiffEqParamEstim

include("pages.jl")

makedocs(
    sitename="DiffEqParamEstim.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqParamEstim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://diffeqparamestim.sciml.ai/stable/"),

    pages=pages
)

deploydocs(
    repo="github.com/SciML/DiffEqParamEstim.jl";
    push_preview=true
)


text:
# Smith Normal Form

function divisable(y::R, x::R ) where {R}
  x == zero(R) && return y == zero(R)
  return div(y,x)*x == y
end

function divide(y::R, x::R) where {R}
    if x != -one(R)
        return div(y,x)
    else
        return y * x
    end
end

function rcountnz(X, i)
    n = size(X, 1)
    c = 0
    for j in 1:n
        if X[j, i] != 0
           c += 1
        end
    end
    return c
end

function ccountnz(X, j)
    m = size(X, 2)
    c = 0
    for i in 1:m
        if X[j, i] != 0
           c += 1
        end
    end
    return c
end

function rswap!(M::AbstractArray{R,2}, r1::Int, r2::Int) where {R}
    r1 == r2 && return M
    m = size(M, 2)
    @inbounds for i in 1:m
        tmp = M[r1, i]
        M[r1, i] = M[r2, i]
        M[r2, i] = tmp
    end
    return M
end

function cswap!(M::AbstractArray{R,2}, c1::Int, c2::Int) where {R}
    c1 == c2 && return M
    n = size(M, 1)
    @inbounds for j in 1:n
        tmp = M[j, c1]
        M[j, c1] = M[j, c2]
        M[j, c2] = tmp
    end
    return M
end

function rowelimination(D::AbstractArray{R,2}, a::R, b::R, c::R, d::R, i::Int, j::Int) where {R}
    m = size(D, 2)
    @inbounds for k in 1:m
        t = D[i,k]
        s = D[j,k]
        D[i,k] = a*t + b*s
        D[j,k] = c*t + d*s
    end
    return D
end

function colelimination(D::AbstractArray{R,2}, a::R, b::R, c::R, d::R, i::Int, j::Int) where {R}
    n = size(D, 1)
    @inbounds for k in 1:n
        t = D[k,i]
        s = D[k,j]
        D[k,i] = a*t + b*s
        D[k,j] = c*t + d*s
    end
    return D
end

function rowpivot(U::AbstractArray{R,2},
                  Uinv::AbstractArray{R,2},
                  D::AbstractArray{R,2},
                  i, j; inverse=true) where {R}
    for k in reverse!(findall(!iszero, view(D, :, j)))
        a = D[i,j]
        b = D[k,j]

        i == k && continue

        s, t, g = bezout(a, b)
        x = divide(a, g)
        y = divide(b, g)

        rowelimination(D, s, t, -y, x, i, k)
        inverse && rowelimination(Uinv, s, t, -y, x, i, k)
        colelimination(U, x, y, -t, s, i, k)
    end
end

function colpivot(V::AbstractArray{R,2},
                  Vinv::AbstractArray{R,2},
                  D::AbstractArray{R,2},
                  i, j; inverse=true) where {R}
    for k in reverse!(findall(!iszero, view(D, i, :)))
        a = D[i,j]
        b = D[i,k]

        j == k && continue

        s, t, g = bezout(a, b)
        x = divide(a, g)
        y = divide(b, g)

        colelimination(D, s, t, -y, x, j, k)
        inverse && colelimination(Vinv, s, t, -y, x, j, k)
        rowelimination(V, x, y, -t, s, j, k)
    end
end

function smithpivot(U::AbstractArray{R,2},
                    Uinv::AbstractArray{R,2},
                    V::AbstractArray{R,2},
                    Vinv::AbstractArray{R,2},
                    D::AbstractArray{R,2},
                    i, j; inverse=true) where {R}

    pivot = D[i,j]
    @assert pivot != zero(R) "Pivot cannot be zero"
    while ccountnz(D,i) > 1 || rcountnz(D,j) > 1
        colpivot(V, Vinv, D, i, j, inverse=inverse)
        rowpivot(U, Uinv, D, i, j, inverse=inverse)
    end
end

function init(M::AbstractSparseMatrix{R,Ti}; inverse=true) where {R, Ti}
    D = copy(M)
    rows, cols = size(M)

    U = spzeros(R, rows, rows)
    for i in 1:rows
        U[i,i] = one(R)
    end
    Uinv = inverse ? copy(U) : spzeros(R, 0, 0)

    V = spzeros(R, cols, cols)
    for i in 1:cols
        V[i,i] = one(R)
    end
    Vinv = inverse ? copy(V) : spzeros(R, 0, 0)

    return U, V, D, Uinv, Vinv
end

function init(M::AbstractMatrix{R}; inverse=true) where {R}
    D = copy(M)
    rows, cols = size(M)

    U = zeros(R, rows, rows)
    for i in 1:rows
        U[i,i] = one(R)
    end
    Uinv = inverse ? copy(U) : zeros(R, 0, 0)

    V = zeros(R, cols, cols)
    for i in 1:cols
        V[i,i] = one(R)
    end
    Vinv = inverse ? copy(V) : zeros(R, 0, 0)

    return U, V, D, Uinv, Vinv
end

function snf(M::AbstractMatrix{R}; inverse=true, debug=false) where {R}
    rows, cols = size(M)
    U, V, D, Uinv, Vinv = init(M, inverse=inverse)

    t = 1
    for j in 1:cols
        debug && println("Working on column $j out of ", size(D,2))
        debug && display(collect(D))

        rcountnz(D,j) == 0 && continue

        prow = 1
        if D[t,t] != zero(R)
            prow = t
        else
            # Good pivot row for j-th column is the one
            # that have a smallest number of elements
            rsize = typemax(R)
            for i in 1:rows
                if D[i,j] != zero(R)
                    c = count(!iszero, view(D, i, :))
                    if c < rsize
                        rsize = c
                        prow = i
                    end
                end
            end
        end

        debug && "Pivot Row selected: t = $t, pivot_row = $prow" |> println
        debug && display(collect(D))
        debug && "Swapping rows. ( $t, $prow )" |> println

        rswap!(D, t, prow)
        inverse && rswap!(Uinv, t, prow)
        cswap!(U, t, prow)

        debug && println("Performing the pivot step at (t=$t, j=$j)")
        debug && display(collect(D))

        smithpivot(U, Uinv, V, Vinv, D, t, j, inverse=inverse)

        cswap!(D, t, j)
        inverse && cswap!(Vinv, t, j)
        rswap!(V, t, j)

        t += 1

        debug && println("D: $(size(D))")
        debug && display(collect(D))
        debug && println("U: $(size(U))")
        debug && display(collect(U))
        debug && println("V: $(size(V))")
        debug && display(collect(V))
        debug && println("Vinv:")
        debug && display(collect(Vinv))
    end

    # Make sure that d_i is divisible be d_{i+1}.
    r = minimum(size(D))
    pass = true
    while pass
        pass = false
        for i in 1:r-1
            divisable(D[i+1,i+1], D[i,i]) && continue
            pass = true
            D[i+1,i] = D[i+1,i+1]

            colelimination(Vinv, one(R), one(R), zero(R), one(R), i, i+1)
            rowelimination(V, one(R), zero(R), -one(R), one(R), i, i+1)

            smithpivot(U, Uinv, V, Vinv, D, i, i, inverse=inverse)
        end
    end

    if issparse(D)
        return dropzeros!(U), dropzeros!(V), dropzeros!(D), dropzeros!(Uinv), dropzeros!(Vinv)
    else
        return U, V, D, Uinv, Vinv
    end
end


text:
using Pkg
Pkg.activate(".")

using Genie, HTTP
using Genie.Router

route("/options", method = OPTIONS) do
  push!(@params(:RESPONSE).headers, "X-Foo-Bar"=>"Baz")
end
Genie.AppServer.startup(verbose = true)

response = HTTP.request("OPTIONS", "http://localhost:8000") # unhandled, should get default response
@show response

response = HTTP.request("OPTIONS", "http://localhost:8000/options") # handled
@show response

exit(0)


text:
function oem_fit(x::AbstractMatrix{<:Real},
                 y::AbstractVector{<:Real};
                 intercept::Bool = true,
                 standardize::Bool = true,
                 nlambda::Int = 100,
                 maxit::Int = 1000,
                 tol::BLAS.BlasFloat = 1e-5,
                 lambda_min_ratio::BLAS.BlasFloat = 1e-3,
                 penalty_factor::AbstractVector{<:Real} = ones(Float64, size(x)[2]))

    n, p = size(x)

    if intercept
        p += 1
        penalty_factor = vcat(0.0, penalty_factor)
    end

    xtx, xty, colsquare_inv, intscale = compute_xtx_y(x, y, intercept, standardize)

    # get eigenvalues, stepsize, A
    steps, A, d = make_steps_A(xtx)

    # initialize beta as vector of zeros
    beta  = zeros(p)

    # compute maximum lambda value
    lambda_max = compute_maxlam(xty, penalty_factor)

    # compute sequence of values for lambda
    lambda = exp.(range(log(lambda_max),
                        stop   = log(lambda_min_ratio * lambda_max),
                        length = nlambda))

    beta_mat  = zeros(p, nlambda)
    iters_vec = repeat([maxit], outer = nlambda)

    # loop over each lambda value
    for l in 1:nlambda
        # set up penalty amounts
        penalty_val   = repeat([lambda[l]], outer = p)
        # adjust penalty by user-specified factors
        penalty_val .*= penalty_factor

        iters = maxit

        # start oem iterations
        for i in 1:maxit
            beta_prev = beta

            u = xty + A * beta_prev

            # M-step
            beta = softthresh(u, penalty_val) ./ steps

            #for j in 1:p
            #    if beta[j] != beta_prev[j]
            #        resid .+= x[:,j] * (beta[j] - beta_prev[j])
            #    end
            #end

            if all(abs.(beta - beta_prev) .<= tol)
                iters = i
                break
            end
        end # end loop over iterations

        ## un-standardize coefficients
        if standardize
            if intercept
                beta[2:p] .*= colsquare_inv[1,:]
            else
                beta .*= colsquare_inv[1,:]
            end
        end

        beta_mat[:,l] = beta
        iters_vec[l] = iters
    end # end loop over lambda

    if intercept
        beta_mat[1,:] ./= intscale
    end

    beta_mat, iters_vec, lambda, d
end


text:
struct Arg4Closure{B,T}
    f::B
    arg4::T
end
(f::Arg4Closure)(arg1, arg2, arg3) = f.f(arg1, arg2, arg3, f.arg4)

function dP!(t, p, out, P)
    B!(t, p, out, P)
    out .= out .+ out' - a(t, P)
    out
end

function partialbridgeode!(S::R3!, t, L, Σnoise, v, νt, Σt, P, ϵ)
    m, d = size(L)
    # print(typeof(Σt))
    # print(typeof(νt))
    # print(typeof(inv(L' * inv(Σ) * L + ϵ * I)   ))
    Γnoise = iszero(Σnoise) ? inv(eps()*one(Σnoise)) : inv(Σnoise)
    Σt[end] .=  inv(L' * Γnoise * L + ϵ * I)
    νt[end] .=  Σt[end] * L' * Γnoise * v
    wsν = workspace(S, νt[end])
    wsΣ = workspace(S, Σt[end])

    for i in length(t)-1:-1:1
        dt = t[i] - t[i+1]
        kernelr3!(Arg4Closure(b!, P), t[i+1], νt[i+1], wsν, νt[i], dt)
        kernelr3!(Arg4Closure(dP!, P), t[i+1], Σt[i+1], wsΣ, Σt[i], dt)
    end
    νt, Σt
 end

struct PartialBridge!{T,TP,TPt,Tv,Tν,TH} <: ContinuousTimeProcess{T}
    Target::TP
    Pt::TPt
    tt::Vector{Float64}
    v::Tv
    ν::Vector{Tν}
    H::Vector{TH}

    function PartialBridge!(tt_, P, Pt, L, v::Tv, ϵ, Σnoise = Bridge.outer(zero(v))) where {Tv}
        tt = collect(tt_)
        N = length(tt)
        m, d = size(L)
        Σ = L'*L
        ν = diag(Σ)
        Σt = [zero(outer(ν)) for i in 1:N]
        νt = [zero(ν) for i in 1:N]

        TH = typeof(Σ)
        Tν = typeof(ν)

        partialbridgeode!(R3!(), tt, L, Σnoise, v, νt, Σt, Pt,ϵ)
        map!(inv, Σt, Σt)
        new{Bridge.valtype(P),typeof(P),typeof(Pt),Tv,Tν,TH}(P, Pt, tt, v, νt, Σt)
    end
end

function _b!((i,t), x, out, P::PartialBridge!)
    b!(t, x, out, P.Target)
    out .+= a(t, x, P.Target)*(P.H[i]*(P.ν[i] - x))
    out
end
σ!(t, x, dw, out, P::PartialBridge!) = σ!(t, x, dw, out, P.Target)

function rti!((t,i), x, out, P::PartialBridge!)
    out .= (P.H[i]*(P.ν[i] - x))
    out
end

btitilde!((t,i), x, out, P) = _b!((t,i), x, out, P.Pt)

constdiff(P::PartialBridge!) = constdiff(P.Target) && constdiff(P.Pt)

function llikelihood(::LeftRule, Xcirc::SamplePath, Po::PartialBridge!; skip = 0)
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    rout = copy(xx[1])
    bout = copy(rout)
    btout = copy(rout)
    for i in 1:length(tt)-1-skip #skip last value, summing over n-1 elements
        s = tt[i]
        x = xx[i]
        rti!((s,i), x, rout, Po)
        b!(s, x, bout, target(Po))
        btitilde!((s,i), x, btout, Po)

        som += dot(bout, rout) * (tt[i+1]-tt[i])
        som -= dot(btout, rout) * (tt[i+1]-tt[i])

        if !constdiff(Po)
            error("not implemented")
        end
    end
    som
end


text:
using CSV
using DataFrames
using MagneticLocSuchowiak
using MLJ
using MLJModels
using NearestNeighbors
using GLM
using StableRNGs
using Plots

DATA_PATH = "C:/Users/Sebastian/PracaInz/data"
CURVES = DATA_PATH * "/curves"
LINES = DATA_PATH * "/lines"
TESTS = DATA_PATH * "/tests"
train_dirs = vcat(readdir(LINES, join=true), readdir(CURVES, join=true))
test_files = readdir(TESTS, join=true)
train_files = collect(Iterators.flatten(readdir.(train_dirs, join=true)))

test = CSV.read("C:/Users/Sebastian/PracaInz/train_data.csv", DataFrame)
a = findall(test[:, :path_id] .== "l1n")
df = test[a, :]

b = findall(test[:, :path_id] .== "l3n")
df2 = test[b, :]

createfullgraph(df, measures=["magnetometr", "acceleromter", "orientation"])
createcomparisongraphs(df, df2, measures=["magnetometr", "acceleromter", "orientation"])

ALL_COLUMNS = [:timestep, :magnetometr_x, :magnetometr_y, :magnetometr_z,
        :acceleromter_x, :acceleromter_y, :acceleromter_z,
        :orientation_x, :orientation_y, :orientation_z,
        :lon, :lat, :path_id, :path_sample]
TRAIN_COLUMNS = [:timestep, :magnetometr_x, :magnetometr_y, :magnetometr_z,
        :lon, :lat, :path_id, :path_sample]
TEST_COLUMNS = [:timestep, :magnetometr_x, :magnetometr_y, :magnetometr_z, :lon, :lat]
traindf = readdata(train_files, columns_to_get=TRAIN_COLUMNS)
testdf = readdata(test_files, columns_to_get=TEST_COLUMNS, add_tag=false)

folds = MagneticLocSuchowiak.getfolds(traindf)

using DecisionTree
knn = @load RandomForestRegressor pkg=DecisionTree
r1 = range(knn, :n_trees; upper=20, lower=1, scale=:log);
knn_x = TunedModel(model=knn, resampling=folds, tuning=Grid(resolution=10), range=[r1], measure=l2);
knn_y = TunedModel(model=knn, resampling=folds, tuning=Grid(resolution=10), range=[r1], measure=l2);

y, X = MagneticLocSuchowiak.preparedata(traindf)
y1 = coerce(y[:,1], Continuous)

m_x = machine(knn_x, X, y1);
m_y = machine(knn_x, X, y[:, 2]);

println(scitype(X))
println(scitype(y1))
println(typeof(y1))
println(typeof(X))
fit!(m_x)
fit!(m_y, verbosity=0)

fitandmean!(knn_x, knn_y, traindf, testdf)

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
# as C's expression shares some similarities with Julia's expression, we just do some tweaks
# on the tokens and then try to parse them using Julia's `Meta.parse`

const C_KEYWORDS_DATATYPE = [
    "char", "double", "float", "int", "long", "short", "signed", "unsigned", "void",
    "_Bool", "_Complex", "_Noreturn"
    ]
const C_KEYWORDS_CVR = ["const", "volatile", "restrict"]
const C_KEYWORDS_UNSUPPORTED = [
    "auto", "break", "case",  "continue", "default", "do",
    "else", "enum", "extern", "for", "goto", "if", "inline", "register",
    "return", "sizeof", "static", "struct", "switch", "typedef", "union",
    "while", "_Alignas", "_Alignof", "_Atomic", "_Decimal128", "_Decimal32",
    "_Decimal64", "_Imaginary", "_Static_assert", "_Thread_local"
]
const C_KEYWORDS = [
    C_KEYWORDS_DATATYPE..., C_KEYWORDS_UNSUPPORTED...,
    "_Generic", "_Decimal128", "_Decimal32", "_Decimal64"
]

const C_DATATYPE_TO_JULIA_DATATYPE = Dict(
    "char"                      => :Cchar,
    "double"                    => :Float64,
    "float"                     => :Float32,
    "int"                       => :Cint,
    "long"                      => :Clong,
    "short"                     => :Cshort,
    "void"                      => :Cvoid,
    "_Bool"                     => :Cuchar,
    "_Complex float"            => :ComplexF32,
    "_Complex double"           => :ComplexF64,
    "_Noreturn"                 => :(Union{}),
    "signed char"               => :Int8,
    "signed int"                => :Cint,
    "signed long"               => :Clong,
    "signed short"              => :Cshort,
    "unsigned char"             => :Cuchar,
    "unsigned int"              => :Cuint,
    "unsigned long"             => :Culong,
    "unsigned short"            => :Cushort,
    "long long"                 => :Clonglong,
    "long long int"             => :Clonglong,
    "signed long long int"      => :Clonglong,
    "unsigned long long"        => :Culonglong,
    "unsigned long long int"    => :Culonglong
)

const LITERAL_LONG = ["L", "l"]
const LITERAL_ULONG = ["UL", "Ul", "uL", "ul", "LU", "Lu", "lU", "lu"]
const LITERAL_ULONGLONG = ["ULL", "Ull", "uLL", "ull", "LLU", "LLu", "llU", "llu"]
const LITERAL_LONGLONG = ["LL", "ll"]
const LITERAL_NON_FLOAT_SUFFIXES = [
    LITERAL_ULONGLONG..., LITERAL_LONGLONG..., LITERAL_ULONG..., LITERAL_LONG..., "U", "u",
]
const LITERAL_FLOAT_SUFFIXES = ["F", "f"]
const LITERAL_SUFFIXES = vcat(LITERAL_NON_FLOAT_SUFFIXES, LITERAL_FLOAT_SUFFIXES)

function c_non_float_literal_to_julia(sfx)
    sfx = lowercase(sfx)
    # integers following http://en.cppreference.com/w/cpp/language/integer_literal
    unsigned = occursin("u", sfx)
    if unsigned && (endswith(sfx, "llu") || endswith(sfx, "ull"))
        return "Culonglong"
    elseif !unsigned && endswith(sfx, "ll")
        return "Clonglong"
    elseif occursin("ul", sfx) || occursin("lu", sfx)
        return "Culong"
    elseif !unsigned && endswith(sfx, "l")
        return "Clong"
    else
        return unsigned ? "Cuint" : "Cint"
    end
end

function c_float_literal_to_julia(txt, sfx)
    sfx = lowercase(sfx)
    # floats following http://en.cppreference.com/w/cpp/language/floating_literal
    float64 = occursin(".", txt) && occursin("l", sfx)  # long double
    float32 = occursin("f", sfx)

    if float64 || float32
        float64 && return "Float64"
        float32 && return "Float32"
    end
end

function normalize_literal_type(text::AbstractString)
    txt = strip(text)
    for sfx in LITERAL_NON_FLOAT_SUFFIXES
        if endswith(txt, sfx)
            txt_no_sfx = replace(txt, Regex(sfx*"\$")=>"")
            type = c_non_float_literal_to_julia(sfx)
            return "$(type)($txt_no_sfx)"
        end
    end
    for sfx in LITERAL_FLOAT_SUFFIXES
        if !occursin("0x", lowercase(txt)) && endswith(txt, sfx)
            txt_no_sfx = replace(txt, Regex(sfx*"\$")=>"")
            type = c_float_literal_to_julia(txt_no_sfx, sfx)
            return "$(type)($txt_no_sfx)"
        end
    end
    # Char to Cchar
    if match(r"^'.*'$", txt) !== nothing
        return "Cchar($txt)"
    end
    return txt
end

function normalize_literal(text)
    m = match(r"^0[0-9]*\d$", text)
    if m !== nothing && m.match !== "0"
        strs = "0o"*normalize_literal_type(text)
    else
        strs = normalize_literal_type(text)
    end
    if occursin('\$', strs)
        return "($(replace(strs, "\$"=>"\\\$")))"
    else
        return "($strs)"
    end
end

normalize_punctuation(text) = text == "/" ? "÷" : text == "^" ? "⊻" : text

is_literal(tok::CLToken) = tok.kind == CXToken_Literal
is_identifier(tok::CLToken) = tok.kind == CXToken_Identifier
is_punctuation(tok::CLToken) = tok.kind == CXToken_Punctuation
is_keyword(tok::CLToken) = tok.kind == CXToken_Keyword
is_comment(tok::CLToken) = tok.kind == CXToken_Comment

const C_UNARY_OPS = ["-", "+", "--", "++", #="&", "*",=# "~", "!"]

const DUMMY_TOKEN = CXToken((0,0,0,0), C_NULL)

function tweak_exprs(dag::ExprDAG, toks::Vector)
    new_toks = []
    i = 1
    while i ≤ length(toks)
        tok = toks[i]
        if is_literal(tok)
            push!(new_toks, Literal(DUMMY_TOKEN, CXToken_Literal, normalize_literal(tok.text)))
            i += 1
        elseif is_punctuation(tok)
            if !isempty(tok.text)
                push!(new_toks, Punctuation(DUMMY_TOKEN, CXToken_Punctuation, normalize_punctuation(tok.text)))
            end
            i += 1
        elseif is_keyword(tok)
            j = i
            while j ≤ length(toks) && is_keyword(toks[j])
                j += 1
            end
            txt = ""
            for k = i:j-1
                txt = txt * " " * toks[k].text
            end
            txt = strip(txt)
            if haskey(C_DATATYPE_TO_JULIA_DATATYPE, txt)
                txt = string(C_DATATYPE_TO_JULIA_DATATYPE[txt])
                push!(new_toks, Keyword(DUMMY_TOKEN, CXToken_Keyword, txt))
                i = j
            else
                push!(new_toks, tok)
                i += 1
            end
        elseif i != 1 && is_identifier(tok)
            # whether this identifier is a type-cast
            if i+1 ≤ length(toks) && is_punctuation(toks[i-1]) && is_punctuation(toks[i+1])
                j = findnext(x->is_punctuation(x) && x.text ∈ C_UNARY_OPS, toks, i)
                if j != nothing && j+1 ≤ length(toks) && is_literal(toks[j+1])
                    s = Symbol(tok.text)
                    if haskey(dag.tags, s) ||
                        (haskey(dag.ids, s) && dag.nodes[dag.ids[s]].type isa AbstractTypedefNodeType)
                        # only support tags and typedefs for now
                        # FIXME: support macros
                        op = toks[j].text
                        toks[j] = Punctuation(DUMMY_TOKEN, CXToken_Punctuation, "")
                        toks[j+1] = Literal(DUMMY_TOKEN, CXToken_Literal, op*toks[j+1].text)
                    end
                end
            end
            push!(new_toks, Identifier(DUMMY_TOKEN, CXToken_Identifier, "($(tok.text))"))
            i += 1
        else
            push!(new_toks, tok)
            i += 1
        end
    end
    return new_toks
end

# `#define SOMETHING`
is_macro_definition_only(toks::Vector) = length(toks) == 1 && is_identifier(toks[1])
is_macro_definition_only(toks::TokenList) = toks.size == 1 && is_identifier(toks[1])
is_macro_definition_only(cursor::CLCursor) = is_macro_definition_only(tokenize(cursor))

# identifier and keyword denylist
const MACRO_IDK_DENYLIST = [
    C_KEYWORDS_UNSUPPORTED...,
    "_Pragma",
    "__attribute__",
    "__typeof__",
    "__inline__",
    "__func__",
    "noexcept",
    "##",
    "#",
]

"""
    is_macro_unsupported(cursor::CLCursor) -> Bool
Return true if the macro should be ignored.
"""
function is_macro_unsupported(cursor::CLCursor)
    for tok in tokenize(cursor)
        is_identifier(tok) && tok.text ∈ MACRO_IDK_DENYLIST && return true
        is_keyword(tok) && tok.text ∈ MACRO_IDK_DENYLIST && return true
        is_punctuation(tok) && tok.text ∈ MACRO_IDK_DENYLIST && return true
    end
    return false
end

function add_spaces_for_macros(lhs, rhs)
    if startswith(rhs, "(")  # handle function call
        if endswith(lhs, "?") || endswith(lhs, ":") # handle trinary operator
            return lhs * " " * rhs
        else
            return lhs * rhs
        end
    else
        return lhs * " " * rhs
    end
end

function get_comment_expr(tokens)
    code = join([tok.text for tok in tokens], " ")
    return Expr(:block, "# Skipping MacroDefinition: " * replace(code, "\n" => "\n#"))
end

"""
    macro_emit!
Emit Julia expression for macros.
"""
function macro_emit! end

macro_emit!(dag::ExprDAG, node::ExprNode, options::Dict) = dag

function macro_emit!(dag::ExprDAG, node::ExprNode{MacroDefinitionOnly}, options::Dict)
    ignore_header_guards = get(options, "ignore_header_guards", true)
    suffixes = get(options, "ignore_header_guards_with_suffixes", [])
    if ignore_header_guards
        endswith(string(node.id), "_H") && return dag
        for suffix in suffixes
            endswith(string(node.id), suffix) && return dag
        end
    end
    if !get(options, "ignore_pure_definition", true)
        tokens = tokenize(node.cursor)
        sym = make_symbol_safe(tokens[1].text)
        push!(node.exprs, Expr(:const, Expr(:(=), sym, :nothing)))
    end
    return dag
end

function macro_emit!(dag::ExprDAG, node::ExprNode{MacroUnsupported}, options::Dict)
    print_comment = get(options, "add_comment_for_skipped_macro", true)
    toks = collect(tokenize(node.cursor))
    print_comment && push!(node.exprs, get_comment_expr(toks))
    return dag
end

"Try to parse the string, return `nothing` if the string is incomplete or has syntax error."
function try_parse(str)
    try
        ex = Meta.parse(str)
        Meta.isexpr(ex, :incomplete) && throw(Meta.ParseError("failed to parse:  $str"))
        return ex
    catch err
        return nothing
    end
end

"Parse a string potentially as a pointer type (`type *``). Issue #311"
function parse_as_pointer(str)
    m = match(r"^(.+?)((?:\s(:?\*|const|volatile))+)$", str)

    base_type = isnothing(m) ? str : m[1]
    ptr_postfix = isnothing(m) ? "" : m[2]
    ex = try_parse(base_type)
    isnothing(ex) && return nothing

    ptr_count = count("*", ptr_postfix)
    for _ in 1:ptr_count
        ex = :(Ptr{$ex})
    end
    return ex
end

function macro_emit!(dag::ExprDAG, node::ExprNode{MacroDefault}, options::Dict)
    mode = get(options, "macro_mode", "basic")
    print_comment = get(options, "add_comment_for_skipped_macro", true)
    ignore_header_guards = get(options, "ignore_header_guards", true)
    suffixes = get(options, "ignore_header_guards_with_suffixes", [])
    if ignore_header_guards
        endswith(string(node.id), "_H") && return dag
        for suffix in suffixes
            endswith(string(node.id), suffix) && return dag
        end
    end

    cursor = node.cursor
    toks = collect(tokenize(cursor))
    tokens = tweak_exprs(dag, toks)

    @assert length(tokens) > 1

    sym = make_symbol_safe(tokens[1].text)
    txts = Vector{String}(undef, length(tokens)-1)
    for (i, tok) in enumerate(collect(tokens)[2:end])
        txts[i] = tok.text
    end
    str = reduce(add_spaces_for_macros, txts)
    ex = parse_as_pointer(str)
    if isnothing(ex)
        print_comment && push!(node.exprs, get_comment_expr(toks))
    else
        push!(node.exprs, Expr(:const, Expr(:(=), sym, ex)))
    end

    return dag
end

function macro_emit!(dag::ExprDAG, node::ExprNode{MacroFunctionLike}, options::Dict)
    print_comment = get(options, "add_comment_for_skipped_macro", true)
    mode = get(options, "macro_mode", "basic")
    allowlist = get(options, "functionlike_macro_allowlist", String[])

    cursor = node.cursor
    toks = collect(tokenize(cursor))
    tokens = tweak_exprs(dag, toks)
    id = tokens[1].text
    mode != "aggressive" && id ∉ allowlist && return dag

    lhs_sym = make_symbol_safe(id)

    i = findfirst(x->x.text == ")", tokens)
    sig_ex = Meta.parse(mapreduce(x->x.text, *, tokens[1:i]))
    sig_ex.args[1] = lhs_sym
    if i == length(tokens)
        push!(node.exprs, Expr(:(=), sig_ex, nothing))
    else
        body_toks = tokens[1+i:end]
        txts = [tok.text for tok in body_toks]
        str = reduce(add_spaces_for_macros, txts)
        try
            push!(node.exprs, Expr(:(=), sig_ex, Meta.parse(str)))
        catch err
            print_comment && push!(node.exprs, get_comment_expr(toks))
        end
    end

    return dag
end


text:
# # How to Make Forecast Plots

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ```julia
# using PowerGraphics
# ```

# # Make a plot of the reserves

# ```julia
# plots = plot_reserves(results)
# ```

# ![this one](plots-4/Stack_Up_Reserves.png)
# ![this one](plots-4/Stack_Down_Reserves.png)

# # Make demand plots from results

# ```julia
# plots = plot_demand(results)
# ```

# ![this one](plots-4/Example_Demand_Plot.png)

# ### Make demand plots from a system

# ```julia
# plots = plot_demand(system)
# ```

# ![this one](plots-4/Example_Demand_Plot_From_System.png)

# ### Make demand plots from a Subsection of Time

# ```julia
# initial_time = Dates.Date(2024, 01, 01, 02, 0, 0)
# horizon = 6
# plots = plot_demand(system; horizon = horizon, initial_time = initial_time)
# ```

# ![this one](plots-4/Example_Demand_Plot_Subsection.png)

# ## Make Demand Plots by Power Systems DataType

# ### sort by PowerLoad, Bus, System

# ```julia
# plots = plot_demand(system; aggregate = System)
# ```

# ![this one](plots-4/Example_Demand_Plot_by_Type.png)

# ### Make demand plots without interpolation between points (stair plot)

# ```julia
# plots = plot_demand(system; stair = true)
# ```

# ![this one](plots-4/Example_Stair_Demand_Plot.png)

# ### Set different colors for the plots

# ```julia
# colors = [:orange :pink :blue :red :grey]
# plot_demand(system; seriescolor = colors)
# ```

# ![this one](plots-4/Example_Demand_Plot_with_Different_Colors.png)

# ### Set a title

# ```julia
# title = "Example Demand Plot with Title"
# plot_demand(system; title = title)
# ```

# ![this one](plots-4/Example_Demand_Plot_with_Title.png)

# ### For saving the plot with the PlotlyJS backend, you can set a different format for saving
# ```julia
# plot_demand(system; save = path, format = "html")
# ```
# Default format for saving is png.
# Optional formats for saving include png, html, and svg.


text:
##########################################################
"""
$(TYPEDEF)

Abstract type for physics.
"""
abstract type AbstractPhysics end

##########################################################
"""
$(TYPEDEF)

Abstract type for user data.
!!! compat
    Will be removed in future versions
"""
abstract type AbstractData end

isdata(::Nothing)=false
isdata(::Any)=true

#
# Dummy callbacks
#
function nofunc(f,u,node,data)
end

function nosrc(f,node,data)
end

function default_storage(f,u,node,data)
    f.=u
end

function nofunc2(f,u,node)
end

function nosrc2(f,node)
end

function default_storage2(f,u,node)
    f.=u
end

function nofunc_generic(f,u,sys)
end

function nofunc_generic_sparsity(sys)
end

##########################################################
"""
````
struct Physics
````

Physics data record with the following fields:

$(TYPEDFIELDS)
"""
struct Physics{Flux<:Function,
               Reaction<:Function,
               Storage<:Function,
               Source<:Function,
               BFlux<:Function,
               BReaction<:Function,
               BSource<:Function,
               BStorage<:Function,
               GenericOperator<:Function,
               GenericOperatorSparsity<:Function,
               Data} <: AbstractPhysics
    """
    Flux between neigboring control volumes: `flux(f,_u,edge)` or `flux(f,_u,edge,data)`
    should return in `f[i]` the flux of species i along the edge joining circumcenters
    of neigboring control volumes. `u=unknowns(_u)` returns a 2D array such that for species i,
    `u[i,1]` and `u[i,2]` contain the unknown values at the corresponding ends of the edge.
    """
    flux::Flux

    """
    Storage term (term under time derivative): `storage(f,u,node)` or `storage(f,u,node,data)`

    It should return in `f[i]` the storage term for the i-th equation. `u[i]` contains the value of
    the i-th unknown.
    """
    storage::Storage

    """
    Reaction term:  `reaction(f,u,node)` or `reaction(f,u,node,data)`

    It should return in `f[i]` the reaction term for the i-th equation. `u[i]` contains the value of
    the i-th unknown.
    """
    reaction::Reaction

    """
    Source term: `source(f,node)` or `source(f,node,data)`.

    It should return the in `f[i]` the value of the source term for the i-th equation.
    """
    source::Source

    """
    Flux between neighboring control volumes on the boundary
    """
    bflux::BFlux

    """
    Boundary reaction term:  `breaction(f,u,node)` or `breaction(f,u,node,data)`
    Similar to reaction, but restricted to the inner or outer boundaries.
    """
    breaction::BReaction

    """
    Boundary source term: `bsource(f,node)` or `bsource(f,node,data)`.

    It should return in `f[i]` the value of the source term for the i-th equation.
    """
    bsource::BSource
    """
    Boundary storage term: `bstorage(f,u,node)` or `bstorage(f,u,node,data)`
    Similar to storage, but restricted to the inner or outer boundaries.
    """
    bstorage::BStorage

    """
    Generic operator  `generic_operator(f,u,sys)`.
    This operator acts on the full solution `u` of a system. Sparsity
    is detected automatically  unless `generic_operator_sparsity` is given.
    """
    generic_operator::GenericOperator

    """
    Function defining the sparsity structure of the generic operator.
    This should return the sparsity pattern of the `generic_operator`.
    """
    generic_operator_sparsity::GenericOperatorSparsity

    """
    User data (parameters).
    This allows to pass various parameters to the callback functions.
    """
    data::Data

    """
    Number of species including boundary species.
    """
    num_species::Int8

end

##########################################################
"""
````
Physics(;num_species=0,
         data=nothing,
         flux,
         reaction,
         storage,
         source,
         breaction,
         bstorage,
         generic,
         generic_sparsity
)
````

Constructor for physics data. For the meaning of the optional function
callbacks, see [`Physics`](@ref).

There are two variants of this constructor. It `data` is given, all callback functions
should accept a last `data` argument. Otherwise, no data are passed explicitely, and it
is assumed that constitutive callbacks take parameters from the closure where the function
is defined.

"""
function Physics(;num_species=0,
                 data=nothing,
                 flux::Function=nofunc,
                 reaction::Function=nofunc,
                 storage::Function=default_storage,
                 source::Function=nosrc,
                 bflux::Function=nofunc,
                 breaction::Function=nofunc,
                 bsource::Function=nosrc,
                 bstorage::Function=nofunc,
                 generic::Function=nofunc_generic,
                 generic_sparsity::Function=nofunc_generic_sparsity
                 )
    if !isdata(data)
        flux==nofunc ? flux=nofunc2 : true
        reaction==nofunc ? reaction=nofunc2 : true
        storage==default_storage ? storage=default_storage2 : true
        source==nosrc ? source=nosrc2 : true
        bflux==nofunc ? bflux=nofunc2 : true
        breaction==nofunc ? breaction=nofunc2 : true
        bsource==nosrc ? bsource=nosrc2 : true
        bstorage==nofunc ? bstorage=nofunc2 : true
    end

    return Physics(flux,
                   storage,
                   reaction,
                   source,
                   bflux,
                   breaction,
                   bsource,
                   bstorage,
                   generic,
                   generic_sparsity,
                   data,
                   Int8(num_species)
                   )
end

"""
$(SIGNATURES)

Check if physics object has data
"""
hasdata(physics::Physics)=isdata(physics.data)

"""
$(SIGNATURES)

Show physics object
"""
function Base.show(io::IO,physics::AbstractPhysics)
    str=@sprintf("VoronoiFVM.Physics(num_species=%d",physics.num_species)
    if isdata(physics.data)
        str=str*", data=$(typeof(physics.data))"
    end
    function addfunc(func,name)
        if func!=nofunc
            str=str*", $(name)=$(nameof(func))"
        end
    end

    for name in fieldnames(typeof(physics))
        if (name!=:num_species)  && (name!=:data) && getfield(physics,name)!=nofunc
             str=str*", $(name)=$(nameof(getfield(physics,name)))"
        end
    end
    str=str*")"
    println(io,str)
end

"""
```
@create_physics_wrappers(physics,node,bnode,edge,bedge)
```

Create wrapper functions around physics callbacks which
fit the API of `ForwardDiff.jacobian!` and pass the `data` parameter
if necessary. These are meant to  be defined
before performing assembly loops. The macro creates the follwing variables:

- wrapper functions: `fluxwrap`,`storagewrap`,`reactionwrap`,`bstoragewrap`, `breactionwrap`
- flag variables: `issource`, `isreaction`,`isbreaction`,`isbstorage`

"""
macro create_physics_wrappers(physics,node,bnode,edge,bedge)
    return quote
        data=$(esc(physics)).data
        if isdata(data)
            global issource=($(esc(physics)).source!=nofunc)
            global isreaction=($(esc(physics)).reaction!=nofunc)
            global isbreaction=($(esc(physics)).breaction!=nofunc)
            global isbstorage=($(esc(physics)).bstorage!=nofunc)
            global isbsource=($(esc(physics)).bsource!=nofunc)
            global isbflux=($(esc(physics)).bflux!=nofunc)

            global fluxwrap=function(y, u)
                y.=0
                $(esc(physics)).flux( rhs($(esc(edge)),y),
                                      unknowns($(esc(edge)),u),
                                      $(esc(edge)),
                                      data)
                nothing
            end

            global reactionwrap=function (y, u)
                y.=0
                ## for ii in ..  uu[node.speclist[ii]]=u[ii]
                $(esc(physics)).reaction(rhs($(esc(node)),y),unknowns($(esc(node)),u),$(esc(node)),data)
                ## for ii in .. y[ii]=y[node.speclist[ii]]
                nothing
            end

            global storagewrap= function(y, u)
                y.=0
                $(esc(physics)).storage(rhs($(esc(node)),y),unknowns($(esc(node)),u),$(esc(node)),data)
                nothing
            end

            global sourcewrap=function(y)
                y.=0
                $(esc(physics)).source(rhs($(esc(node)),y),$(esc(node)),data)
                nothing
            end

            global bfluxwrap=function(y, u)
                y.=0
                $(esc(physics)).bflux(rhs($(esc(bedge)),y),unknowns($(esc(bedge)),u),$(esc(bedge)),data)
                nothing
            end

            global breactionwrap=function(y, u)
                y.=0
                $(esc(physics)).breaction(rhs($(esc(bnode)),y),unknowns($(esc(bnode)),u),$(esc(bnode)),data)
                nothing
            end

            global bsourcewrap=function(y)
                y.=0
                $(esc(physics)).bsource(rhs($(esc(bnode)),y),$(esc(bnode)),data)
                nothing
            end

            global bstoragewrap=function(y, u)
                y.=0
                $(esc(physics)).bstorage(rhs($(esc(bnode)),y),unknowns($(esc(bnode)),u),$(esc(bnode)),data)
                nothing
            end

        else
            global issource    = ($(esc(physics)).source    != nofunc2)
            global isreaction  = ($(esc(physics)).reaction  != nofunc2)
            global isbreaction = ($(esc(physics)).breaction != nofunc2)
            global isbsource   = ($(esc(physics)).source    != nofunc2)
            global isbflux     = ($(esc(physics)).bflux     != nofunc2)
            global isbstorage  = ($(esc(physics)).bstorage  != nofunc2)

            global fluxwrap=function(y, u)
                y.=0
                $(esc(physics)).flux( rhs(     $(esc(edge)),y),
                                      unknowns($(esc(edge)),u),
                                      $(esc(edge)))
                nothing
            end

            global reactionwrap=function(y, u)
                y.=0
                ## for ii in ..  uu[node.speclist[ii]]=u[ii]
                $(esc(physics)).reaction(rhs($(esc(node)),y),unknowns($(esc(node)),u),$(esc(node)))
                ## for ii in .. y[ii]=y[node.speclist[ii]]
                nothing
            end

            global storagewrap= function(y, u)
                y.=0
                $(esc(physics)).storage(rhs($(esc(node)),y),unknowns($(esc(node)),u),$(esc(node)))
                nothing
            end

            global sourcewrap=function(y)
                y.=0
                $(esc(physics)).source(rhs($(esc(node)),y),$(esc(node)))
                nothing
            end

            global bfluxwrap=function(y, u)
                y.=0
                $(esc(physics)).bflux(rhs($(esc(bedge)),y),unknowns($(esc(bedge)),u),$(esc(bedge)))
                nothing
            end

            global breactionwrap=function(y, u)
                y.=0
                $(esc(physics)).breaction(rhs($(esc(bnode)),y),unknowns($(esc(bnode)),u),$(esc(bnode)))
                nothing
            end

            global bsourcewrap=function(y)
                y.=0
                $(esc(physics)).bsource(rhs($(esc(bnode)),y),$(esc(bnode)))
                nothing
            end

            global bstoragewrap=function(y, u)
                y.=0
                $(esc(physics)).bstorage(rhs($(esc(bnode)),y),unknowns($(esc(bnode)),u),$(esc(bnode)))
                nothing
            end

        end
    end
end


text:
using Test
using LinearAlgebra
using Nystrom
using StaticArrays

@testset "Basic tests" begin
    pde   = Helmholtz(;dim=3,k=1)
    @gmsh begin
        Geometry.clear_entities!()
        gmsh.initialize()
        GmshSDK.set_meshsize(0.1)
        Geometry.clear_entities!()
        gmsh.model.occ.addSphere(0,0,0,1)
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate(3)
        Ω = GmshSDK.domain(dim=3)
        M = GmshSDK.meshgen(Ω,dim=3)
        return Ω,M
    end
    Γ     = boundary(Ω)
    mesh  = NystromMesh(view(M,Γ),order=1)
    𝐒     = SingleLayerOperator(pde,mesh)
    𝐃     = DoubleLayerOperator(pde,mesh)
    @test Nystrom.kernel_type(𝐒) == Nystrom.SingleLayer()
    @test Nystrom.kernel_type(𝐃) == Nystrom.DoubleLayer()
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
function validate(a::Int, history::Dict{Int,Byte})::Bool
    for g in keys(history)
        if fast_wordle_score(g, a) != history[g]
            return false
        end
    end
    return true

end

"""
    filter_answers(history::Dict{Int,Byte})::Vector{Int}
Return a list of indices into ANS_LIST of words that could be valid given the history
"""
function filter_answers(history::Dict{Int,Byte})::Vector{Int}
    [a for a = 1:NA if validate(a, history)]
end

export filter_answers, validate

"""
    make_counts(g::Int, possibles::Vector{Int})
Given a guess `g` and a list of possible answers `possibles`
create a `Counter` that shows how often each score occurs.
"""
function make_counts(g::Int, possibles::Vector{Int})::Counter{Byte}
    vals = [fast_wordle_score(g, a) for a in possibles]
    return counter(vals)
end

export make_counts

# These are three scoring functions to evaluate a list of possible

"""
    min_max_score(g, possibles)
This score is the largest count returned by `make_counts`.
"""
function min_max_score(g::Int, possibles::Vector{Int})::Int
    c = make_counts(g, possibles)
    return maximum(values(c))
end

"""
    entropy_score(g, possibles)
This is the sum of `x * log(x)` where `x` is a counter value
returned by `make_counts`.
"""
function entropy_score(g::Int, possibles::Vector{Int})
    c = make_counts(g, possibles)
    return sum(x * log(x) for x in values(c))
end

"""
    equi_score(g, possibles)
Always return `1`. This treats all guesses the same.
"""
equi_score(g::Int, possibles::Vector{Int}) = 1

"""
    L2_score(g, possibles)
Sum of `x^2` where `x` is a counter values returned by
`make_counts`.
"""

function L2_score(g::Int, possibles::Vector{Int})
    c = make_counts(g, possibles)
    return sum(x * x for x in values(c))
end

"""
    best_guess(possibles, score_func)
Given a list of possible answers (given by `filter_answers`), return a best guess
based on minimizing `score_func`.
"""
function best_guess(possibles::Vector{Int}, score_func::Function = min_max_score)
    glist = shuffle(collect(1:NG))

    best_g = first(possibles)
    best_s = score_func(best_g, possibles)

    for g in glist
        s = score_func(g, possibles)
        if s < best_s
            best_g = g
            best_s = s
        end
    end
    return best_g
end

export best_guess, min_max_score, entropy_score, equi_score, L2_score

FIRST_GUESS = Dict{Function,Int}()
# @info "Precomputing first guesses"

FIRST_GUESS[min_max_score] = 1535
FIRST_GUESS[equi_score] = 1
FIRST_GUESS[entropy_score] = 10846
FIRST_GUESS[L2_score] = 9958

function make_first_guesses()
    for f in [min_max_score, entropy_score, equi_score, L2_score]
        g = best_guess(collect(1:NA), f)
        FIRST_GUESS[f] = g
        guess = GUESS_LIST[g]
        @info "First guess for $f is $guess"
    end
end

export make_first_guesses


text:
"""
    LaneChangeChoice
A choice of whether to change lanes, and what direction to do it in
"""
const DIR_RIGHT = -1
const DIR_MIDDLE =  0
const DIR_LEFT =  1
struct LaneChangeChoice
    dir::Int # -1, 0, 1
end
Base.show(io::IO, a::LaneChangeChoice) = @printf(io, "LaneChangeChoice(%d)", a.dir)

function get_lane_offset(a::LaneChangeChoice, scene::Scene, roadway::Roadway, vehicle_index::Int)
    if a.dir == DIR_MIDDLE
        return posf(scene[vehicle_index].state).t
    elseif a.dir == DIR_LEFT
        return lane_offset_left(roadway, scene[vehicle_index])
    else
        @assert(a.dir == DIR_RIGHT)
        return lane_offset_right(roadway, scene[vehicle_index])
    end
end

####################

abstract type LaneChangeModel{A} <: DriverModel{A} end


text:
u, v, w = ℝ³(rand(3)), ℝ³(rand(3)), ℝ³(rand(3))
a, b = ℝ³(rand(3)), ℝ³(rand(3))
zero = ℝ³([0.0; 0.0; 0.0])
α, β = rand(2)

## Abstract Vector Space Tests ##

@test isapprox(u + (v + w), (u + v) + w) # 1. associativity of addition
@test isapprox(u + v, v + u) # 2. commutativity of addition
@test isapprox(u + zero, u) # 3. the zero vector
@test isapprox(u - u, zero) # 4. the inverse element
@test isapprox(α * (u + v), α * u + α * v) # 5. distributivity Ι
@test isapprox((α + β) * u, α * u + β * u) # 6. distributivity ΙΙ
@test isapprox(α * (β * u), (α * β) * u) # 7. associativity of scalar multiplication
@test isapprox(1u, u) # 8. the unit scalar 1

## Product Space Tests ##

@test isapprox(dot(u, v), dot(v, u)) # symmetric
@test isapprox(dot(u, α * a + β * b), α * dot(u, a) + β * dot(u, b)) # linear
@test dot(u, u) ≥ 0 && isapprox(dot(u, zero), 0) # positive semidefinite
@test size(outer(u, v)) == (length(u), length(v))

## isapprox

segments = 5
array1 = Array{ℝ³,2}(undef, segments, segments)
for i in 1:segments
    for j in 1:segments
        array1[i, j] = ℝ³(rand(3))
    end
end

array2 = map(x -> 2x, array1)

@test isapprox(array1, array1, atol = TOLERANCE)
@test isapprox(array1, array2, atol = TOLERANCE) == false

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using Octavian, VectorizationBase
const F64 = Octavian.StaticFloat64
function matmul_pack_ab!(C, A, B, ::Val{W₁}, ::Val{W₂}, ::Val{R₁}, ::Val{R₂}) where {W₁, W₂, R₁, R₂}
  M, N = size(C); K = size(B,1)
  zc, za, zb = Octavian.zstridedpointer.((C,A,B))
  nspawn = VectorizationBase.num_cores()
  t = Inf
  GC.@preserve C A B begin
    for _ ∈ 1:2
      t = min(t, @elapsed(Octavian.matmul_pack_A_and_B!(zc, za, zb, Octavian.One(), Octavian.Zero(), M, K, N, Int(nspawn), F64(W₁), F64(W₂), F64(R₁), F64(R₂))))
    end
  end
  return t
end

function bench_size(Cs, As, Bs, ::Val{W₁}, ::Val{W₂}, ::Val{R₁}, ::Val{R₂}) where {W₁, W₂, R₁, R₂}
    if length(first(Cs)) < length(last(Cs))
        matmul_pack_ab!(first(Cs), first(As), first(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    else
        matmul_pack_ab!(last(Cs), last(As), last(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    end
    repeat = 1
    gflop = 0.0
    for _ ∈ 1:repeat
        for (C,A,B) ∈ zip(Cs,As,Bs)
            M, K, N = Octavian.matmul_sizes(C, A, B)
            # sleep(0.5)
            t = matmul_pack_ab!(C, A, B, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
            gf = 2e-9M*K*N / t
            gflop += gf
        end
    end
    gflop / (repeat*length(As))
end
matrix_sizes(s::Int) = (s,s,s)
matrix_sizes(MKN::NTuple{3,Int}) = MKN
size_range(l, u, len) = round.(Int, exp.(range(log(l), stop = log(u), length = len)))
function matrix_range(l, u, len, ::Type{T} = Float64) where {T}
    matrix_range(size_range(l, u, len), T)
end
function matrix_range(S, ::Type{T} = Float64) where {T}
    Alen = 0; Blen = 0; Clen = 0;
    for s ∈ S
        M, K, N = matrix_sizes(s)
        Alen = max(Alen, M*K)
        Blen = max(Blen, K*N)
        Clen = max(Clen, M*N)
    end
    Abuf = rand(T, Alen)
    Bbuf = rand(T, Blen)
    Cbuf = rand(T, Clen)
    As = Vector{Base.ReshapedArray{T, 2, SubArray{T, 1, Vector{T}, Tuple{Base.OneTo{Int}}, true}, Tuple{}}}(undef, length(S))
    Bs = similar(As); Cs = similar(As);
    for (i,s) ∈ enumerate(S)
        M, K, N = matrix_sizes(s)
        As[i] = reshape(view(Abuf, Base.OneTo(M * K)), (M, K))
        Bs[i] = reshape(view(Bbuf, Base.OneTo(K * N)), (K, N))
        Cs[i] = reshape(view(Cbuf, Base.OneTo(M * N)), (M, N))
    end
    Cs, As, Bs
end

T = Float64
min_size = round(Int, sqrt(0.65 * Octavian.VectorizationBase.cache_size(Val(3)) / sizeof(T)))
max_size = round(Int, sqrt( 32  * Octavian.VectorizationBase.cache_size(Val(3)) / sizeof(T)))

SR = size_range(max_size, min_size, 400);
const CsConst, AsConst, BsConst = matrix_range(SR, T);

# using Hyperopt
# ho = @hyperopt for i = 100, sampler=GPSampler(Max),
#     W₁ = LinRange(0.001, 0.3, 1000),
#     W₂ = LinRange(0.01, 2.0, 1000),
#     R₁ = LinRange(0.3, 0.9, 1000),
#     R₂ = LinRange(0.4, 0.99, 1000)
#     print("Params: ", (W₁, W₂, R₁, R₂), "; ")
#     gflop = bench_size(CsConst, AsConst, BsConst, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
#     println(gflop)
#     gflop
# end

# function restart(ho, iterations = 100)
#     olditer = ho.iterations
#     ho2 = Hyperoptimizer(
#         olditer + iterations,
#         ho.params,
#         ho.candidates,
#         ho.history[1:olditer],
#         ho.results[1:olditer],
#         ho.sampler#,
#         # ho.objective
#     )
#     # Hyperopt.optimize(ho2)
#     for nt ∈ ho2
#         (i, W₁, W₂, R₁, R₂) = nt
#         print("Params: ", (W₁, W₂, R₁, R₂), "; ")
#         gflop = bench_size(CsConst, AsConst, BsConst, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
#         println(gflop)
#         push!(ho2.results, gflop)
#     end
#     ho2
# end
# ho2 = restart(ho, 400)

function matmul_objective(params)
    print("Params: ", params, "; ")
    W₁, W₂, R₁, R₂ = params
    gflop = bench_size(CsConst, AsConst, BsConst, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    println(gflop)
    - gflop
end
using Optim
hours = 60.0*60.0; days = 24hours;
init = Float64[Octavian.W₁Default(), Octavian.W₂Default(), Octavian.R₁Default(), Octavian.R₂Default()]
lower = 0.75 .* init;
upper = [1.25init[1], 1.25init[2], 0.75*init[3] + 0.25, 0.75*init[4] + 0.25];
# init = [0.001, 0.9754033943603924, 0.5711159869399494, 0.7547361860432168];

opt = Optim.optimize(
    matmul_objective, init, ParticleSwarm(lower = lower, upper = upper),
    Optim.Options(iterations = 10^6, time_limit = 8hours)
);


text:
@testset "download_filing() (Download Filing)" begin

    file_name = "edgar/data/880794/9999999997-05-050434.txt"
    temp_file = "./temp_filing.txt"
    dest = "./"
    download_filing(file_name, temp_file, dest)
    @test isfile(temp_file)
    rm(temp_file)

end


text:
#=
P02 (*) Find the last but one box of a list.
Example:
* (my-but-last '(a b c d))
(C D)
=#

function my_but_last(coll)
    coll[end - 1:end]
end

@assert my_but_last((:a, :b, :c, :d)) == (:c, :d)
@assert my_but_last((1, 2, 3)) == (2, 3)
@assert my_but_last(((:a, :b), (:c, :d), (:e, :f))) == ((:c, :d), (:e, :f))

println("Tests passed: JL-02.jl")


text:
function estimate(get_sample_from_input, inputs, rel_tolerance, min_samples, max_samples; print_progress = false, return_all_samples = false)
    shuffled_inputs = shuffle(inputs)
    max_samples = min(max_samples, length(inputs))

    i = 0

    estimate(_get_new_estimation, _get_estimation_variance, rel_tolerance, min_samples, max_samples,
        print_progress = print_progress, return_all_samples = return_all_samples) do

        i += 1
        return get_sample_from_input(shuffled_inputs[i])
    end
end

function estimate(get_sample, rel_tolerance, min_samples, max_samples; print_progress = false, return_all_samples = false)
    estimate(get_sample, _get_new_estimation, _get_estimation_variance, rel_tolerance, min_samples, max_samples, print_progress = print_progress, return_all_samples = return_all_samples)
end

function _get_new_estimation(sample::Number, previous_samples, previous_estimation)
    if length(previous_samples) == 0
        sample
    else
        (sum(previous_samples) + sample) / (length(previous_samples) + 1)
    end
end
function _get_new_estimation(sample::Vector, previous_samples, previous_estimation)
    if length(previous_samples) == 0
        mean(sample)
    else
        (sum(previous_samples) + sum(sample)) / (length(previous_samples) + length(sample))
    end
end

function _get_estimation_variance(sample, previous_samples, estimate)
    if length(previous_samples) + length(sample) <= 1 return Inf end

    all_samples = [previous_samples ; sample]

    n = length(all_samples)
    mean = sum(all_samples) / n

    return 1 / (n * (n - 1)) * sum([(s - mean)^2 for s in all_samples]) / max(EPS, mean^2)
end

function estimate(get_sample, get_new_estimation, get_estimation_variance, tolerance, min_samples, max_samples; print_progress = false, return_all_samples = false)
    previous_samples = []
    num_samples = 0

    estimation = 0
    estimation_variance = 2*tolerance

    while (num_samples <= min_samples || tolerance < estimation_variance) && num_samples <= max_samples
        try
            new_sample = get_sample()

            estimation = get_new_estimation(new_sample, previous_samples, estimation)
            estimation_variance = get_estimation_variance(new_sample, previous_samples, estimation)

            _add_sample!(previous_samples, new_sample)
            num_samples = length(previous_samples)

            if print_progress
                println("$estimation_variance \t $new_sample \t $estimation")
            end
        catch e
            if !isa(e, SimulationExpection)
                rethrow(e)
            end
        end
    end

    if return_all_samples
        return estimation, previous_samples
    else
        return estimation
    end
end

_add_sample!(previous_samples::Vector, new_sample) = push!(previous_samples, new_sample)
_add_sample!(previous_samples::Vector, new_sample::Vector) = append!(previous_samples, new_sample)

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark