text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestSVD

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted

@testset "Simple svdvals / svd tests" begin
    ≊(x,y) = isapprox(x,y,rtol=1e-15)

    m = [2, 0]
    @test @inferred(svdvals(m)) ≊ [2]
    @test @inferred(svdvals!(float(m))) ≊ [2]
    for sf in (@inferred(svd(m)), @inferred(svd!(float(m))))
        @test sf.S ≊ [2]
        @test sf.U'sf.U ≊ [1]
        @test sf.Vt'sf.Vt ≊ [1]
        @test sf.U*Diagonal(sf.S)*sf.Vt' ≊ m
    end
    F = @inferred svd(m, full=true)
    @test size(F.U) == (2, 2)
    @test F.S ≊ [2]
    @test F.U'F.U ≊ Matrix(I, 2, 2)
    @test F.Vt'*F.Vt ≊ [1]
    @test @inferred(svdvals(3:4)) ≊ [5]
    A = Matrix(1.0I, 2, 2)
    Z = svd(Hermitian(A); full=true)
    @test Z.S ≈ ones(2)
    @test Z.U'Z.U ≈ I(2)

    m1 = [2 0; 0 0]
    m2 = [2 -2; 1 1]/sqrt(2)
    m2c = Complex.([2 -2; 1 1]/sqrt(2))
    @test @inferred(svdvals(m1))  ≊ [2, 0]
    @test @inferred(svdvals(m2))  ≊ [2, 1]
    @test @inferred(svdvals(m2c)) ≊ [2, 1]

    sf1 = @inferred svd(m1)
    sf2 = @inferred svd(m2)
    @test sf1.S ≊ [2, 0]
    @test sf2.S ≊ [2, 1]
    # U & Vt are unitary
    I22 = Matrix(I, 2, 2)
    @test sf1.U*sf1.U'   ≊ I22
    @test sf1.Vt*sf1.Vt' ≊ I22
    @test sf2.U*sf2.U'   ≊ I22
    @test sf2.Vt*sf2.Vt' ≊ I22
    # SVD not uniquely determined, so just test we can reconstruct the
    # matrices from the factorization as expected.
    @test sf1.U*Diagonal(sf1.S)*sf1.Vt' ≊ m1
    @test sf2.U*Diagonal(sf2.S)*sf2.Vt' ≊ m2

    @test ldiv!([0., 0.], svd(Matrix(I, 2, 2)), [1., 1.]) ≊ [1., 1.]
    @test inv(svd(Matrix(I, 2, 2))) ≈ I
    @test inv(svd([1 2; 3 4])) ≈ [-2.0 1.0; 1.5 -0.5]
    @test inv(svd([1 0 1; 0 1 0])) ≈ [0.5 0.0; 0.0 1.0; 0.5 0.0]
    @test_throws SingularException inv(svd([0 0; 0 0]))
    @test inv(svd([1+2im 3+4im; 5+6im 7+8im])) ≈ [-0.5 + 0.4375im 0.25 - 0.1875im; 0.375 - 0.3125im -0.125 + 0.0625im]
end

n = 10

Random.seed!(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int)
    aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    asym = aa' + aa                 # symmetric indefinite
    for a in (aa, view(aa, 1:n, 1:n))
        usv = svd(a)
        @testset "singular value decomposition" begin
            @test usv.S === svdvals(usv)
            @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ a
            @test convert(Array, usv) ≈ a
            @test usv.Vt' ≈ usv.V
            @test_throws ErrorException usv.Z
            b = rand(eltya,n)
            @test usv\b ≈ a\b
            @test Base.propertynames(usv) == (:U, :S, :V, :Vt)
            @test size(usv) == size(a)
            if eltya <: BlasFloat
                svdz = svd!(Matrix{eltya}(undef,0,0))
                @test svdz.U ≈ Matrix{eltya}(I, 0, 0)
                @test svdz.S ≈ real(zeros(eltya,0))
                @test svdz.Vt ≈ Matrix{eltya}(I, 0, 0)
            end
        end
        @testset "singular value decomposition of adjoint/transpose" begin
            for transform in (adjoint, transpose)
                usv = svd(transform(a))
                @test usv.S === svdvals(usv)
                @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ transform(a)
                @test convert(Array, usv) ≈ transform(a)
                @test usv.Vt' ≈ usv.V
                @test_throws ErrorException usv.Z
                b = rand(eltya,n)
                @test usv\b ≈ transform(a)\b
            end
        end
        @testset "Generalized svd" begin
            a_svd = a[1:div(n, 2), :]
            gsvd = svd(a,a_svd)
            @test Base.propertynames(gsvd) == (:alpha, :beta, :vals, :S, :D1, :D2, :R0, :U, :V, :Q, :a, :b, :k, :l, :R)
            @test gsvd.U*gsvd.D1*gsvd.R*gsvd.Q' ≈ a
            @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ a_svd
            @test usv.Vt' ≈ usv.V
            @test_throws ErrorException usv.Z
            @test_throws ErrorException gsvd.Z
            @test gsvd.vals ≈ svdvals(a,a_svd)
            α = eltya == Int ? -1 : rand(eltya)
            β = svd(α)
            @test β.S == [abs(α)]
            @test svdvals(α) == abs(α)
            u,v,q,d1,d2,r0 = svd(a,a_svd)
            @test u ≈ gsvd.U
            @test v ≈ gsvd.V
            @test d1 ≈ gsvd.D1
            @test d2 ≈ gsvd.D2
            @test q ≈ gsvd.Q
            @test gsvd.a.^2 + gsvd.b.^2 ≈ fill(1, length(gsvd.a))
            @test gsvd.alpha.^2 + gsvd.beta.^2 ≈ ones(eltya, length(gsvd.a))
            #testing the other layout for D1 & D2
            b = rand(eltya,n,2*n)
            c = rand(eltya,n,2*n)
            gsvd = svd(b,c)
            @test gsvd.U*gsvd.D1*gsvd.R*gsvd.Q' ≈ b
            @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ c
        end
    end
    @testset "singular value decomposition of Hermitian/real-Symmetric" begin
        for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,))
            usv = svd(T(asym))
            @test usv.S === svdvals(usv)
            @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ T(asym)
            @test convert(Array, usv) ≈ T(asym)
            @test usv.Vt' ≈ usv.V
            @test_throws ErrorException usv.Z
            b = rand(eltya,n)
            @test usv\b ≈ T(asym)\b
        end
    end
    if eltya <: LinearAlgebra.BlasReal
        @testset "Number input" begin
            x, y = randn(eltya, 2)
            @test svd(x)    == svd(fill(x, 1, 1))
            @test svdvals(x)    == first(svdvals(fill(x, 1, 1)))
            @test svd(x, y) == svd(fill(x, 1, 1), fill(y, 1, 1))
            @test svdvals(x, y) ≈  first(svdvals(fill(x, 1, 1), fill(y, 1, 1)))
        end
    end
    if eltya != Int
        @testset "isequal, ==, and hash" begin
            x, y   = rand(eltya), convert(eltya, NaN)
            Fx, Fy = svd(x), svd(y)
            @test   Fx == Fx
            @test !(Fy == Fy)
            @test isequal(Fy, Fy)
            @test hash(Fx)          == hash(Fx)
            @test hash(Fx, UInt(1)) == hash(Fx, UInt(1))
            @test hash(Fy)          == hash(Fy)
            @test hash(Fy, UInt(1)) == hash(Fy, UInt(1))
        end
    end
end

@testset "SVD Algorithms" begin
    ≊(x,y) = isapprox(x,y,rtol=1e-15)

    x = [0.1 0.2; 0.3 0.4]

    for alg in [LinearAlgebra.QRIteration(), LinearAlgebra.DivideAndConquer()]
        sx1 = svd(x, alg = alg)
        @test sx1.U * Diagonal(sx1.S) * sx1.Vt ≊ x
        @test sx1.V * sx1.Vt ≊ I
        @test sx1.U * sx1.U' ≊ I
        @test all(sx1.S .≥ 0)

        sx2 = svd!(copy(x), alg = alg)
        @test sx2.U * Diagonal(sx2.S) * sx2.Vt ≊ x
        @test sx2.V * sx2.Vt ≊ I
        @test sx2.U * sx2.U' ≊ I
        @test all(sx2.S .≥ 0)
    end
end

@testset "REPL printing of SVD" begin
    svdd = svd(randn(3, 3))
    svdstring = sprint((t, s) -> show(t, "text/plain", s), svdd)
    ustring = sprint((t, s) -> show(t, "text/plain", s), svdd.U)
    sstring = sprint((t, s) -> show(t, "text/plain", s), svdd.S)
    vtstring = sprint((t, s) -> show(t, "text/plain", s), svdd.Vt)
    @test svdstring == "$(summary(svdd))\nU factor:\n$ustring\nsingular values:\n$sstring\nVt factor:\n$vtstring"
end

@testset "REPL printing of Generalized SVD" begin
    a = randn(3, 3)
    b = randn(3, 3)
    svdd = svd(a, b)
    svdstring = sprint((t, s) -> show(t, "text/plain", s), svdd)
    ustring = sprint((t, s) -> show(t, "text/plain", s), svdd.U)
    qstring = sprint((t, s) -> show(t, "text/plain", s), svdd.Q)
    vstring = sprint((t, s) -> show(t, "text/plain", s), svdd.V)
    d1string = sprint((t, s) -> show(t, "text/plain", s), svdd.D1)
    d2string = sprint((t, s) -> show(t, "text/plain", s), svdd.D2)
    r0string = sprint((t, s) -> show(t, "text/plain", s), svdd.R0)
    @test svdstring == "$(summary(svdd))\nU factor:\n$ustring\nV factor:\n$vstring\nQ factor:\n$qstring\nD1 factor:\n$d1string\nD2 factor:\n$d2string\nR0 factor:\n$r0string"
end

@testset "c-tor with varying input eltypes" begin
    A = randn(Float64, 10, 10)
    U, S, V = svd(A)
    Ut = convert.(Float16, U)
    Vt = convert.(Float32, V)
    svdc = SVD{ComplexF32}(Ut, S, Vt)
    @test svdc isa SVD{ComplexF32}
    Uc, Sc, Vc = svdc
    @test Uc * diagm(0=>Sc) * transpose(V) ≈ complex.(A) rtol=1e-3
end

@testset "adjoint of SVD" begin
    n = 5
    B = randn(5, 2)

    @testset "size(b)=$(size(b))" for b in (B[:, 1], B)
        @testset "size(A)=$(size(A))" for A in (
            randn(n, n),
            # Wide problems become minimum norm (in x) problems similarly to LQ
            randn(n + 2, n),
            randn(n - 2, n),
            complex.(randn(n, n), randn(n, n)))

            F = svd(A)
            x = F'\b
            @test x ≈ A'\b
            @test length(size(x)) == length(size(b))
        end
    end
end

end # module TestSVD

max_stars_repo_name: GiggleLiu/julia

max_issues_repo_name: GiggleLiu/julia

max_forks_repo_name: GiggleLiu/julia


text:
### Initialize population
function initialize_population!(f::Function, lb, ub, method::Clearing, options::Options)

	# Parameters
	n = length(lb)

	# Pre-allocate pop
	method.population = zeros(method.nind * 2, n + 3)

	# Init pop
	if options.multithreads
		Threads.@threads for i in 1:method.nind
		    method.population[i,1:n] = lb .+ rand(n) .* (ub .- lb)
		    method.population[i,n+1] = unidrnd(3)
		    method.population[i,n+2] = 0
		    method.population[i,n+3] = f(method.population[i,1:n])
		end
	else
		for i in 1:method.nind
		    method.population[i,1:n] = lb .+ rand(n) .* (ub .- lb)
		    method.population[i,n+1] = unidrnd(3)
		    method.population[i,n+2] = 0
		    method.population[i,n+3] = f(method.population[i,1:n])
		end
	end
end

### Update population
function update_population!(f::Function, lb, ub, method::Clearing, options::Options)

	# Selection
	selection!(method)

	# Crossover
	crossover!(lb, ub, method)

	# Mutation
	mutation!(lb, ub, method)

	# Population assessment
	assessment!(f, method, options)

	# Eclaircissement
	clearing!(lb, ub, method)

	# Density indexes to the niches
	set_density!(method)
end

### Selection
function selection!(method::Clearing)
	# Parameter
	n = size(method.population, 2) - 3

    for i in 1:method.nind
		# Select randomly two parents from the parent pool
		i1 = unidrnd(method.nind)
        i2 = unidrnd(method.nind)
		# Individuals must be different...
        while i1 == i2
          i2 = unidrnd(method.nind)
        end
		# If i2 is stronger than i1, select i2
        method.population[i2,n+3] >= method.population[i1,n+3] && method.population[i2,n+2] > method.population[i1,n+2] ? i1 = i2 : nothing
		# Copy to the new generation
        method.population[method.nind+i,:] = method.population[i1,:]
    end
end

### Crossover
function crossover!(lb, ub, method::Clearing)
	# Parameter
	n = length(lb)

    for i in (method.nind + 1) : (2 * method.nind - 1)
		# Autoadaptative crossover
		if method.crossover == "auto"
			# Select randomly the crossover type
			bitrand(1)[1] ? select = method.population[i,n+1] : select = method.population[i+1,n+1]
			if select == 1.
				BLX!(i, lb, ub, method)
			elseif select == 2.
				vSBX!(i, lb, ub, method)
		    elseif select == 3.
				PNX!(i, lb, ub, method)
		    end
		# BLX crossover
		elseif method.crossover == "BLX"
			BLX!(i, lb, ub, method)
		# vSBX crossover
		elseif method.crossover == "vSBX"
			vSBX!(i, lb, ub, method)
		# PNX crossover
		elseif method.crossover == "PNX"
			PNX!(i, lb, ub, method)
		else
			println("Crossover method unknown...");
		end
	 end
end
# BLX-dec crossover
function BLX!(i, lb, ub, method::Clearing)
	# Parameter
	n = length(lb)
	# Chose randomly between [dec, 2 * dec + 1]
	δ = .- method.dec .+ rand(n) .* (2 .* method.dec .+ 1)
	#
	E1 = method.population[i,1:n] .+ δ .* (method.population[i+1,1:n] .- method.population[i,1:n])
	E2 = method.population[i+1,1:n] .+ δ .* (method.population[i,1:n] .- method.population[i+1,1:n])
	# Copy the children in the population with overshoot tests
	method.population[i,1:n] = max.(min.(ub, E1), lb)
	method.population[i+1,1:n] = max.(min.(ub, E2), lb)
	# Update crossover type
	method.population[i,n+1], method.population[i+1,n+1] = 1., 1.
end
# vSBX crossover
function vSBX!(i, lb, ub, method::Clearing)
	# Parameter
	n = length(lb)
	E1, E2 = zeros(n), zeros(n)

	for j in 1:n
		# Select randomly u between ]0,1[
		u = rand(1)[1]
		while u == 0. || u == 1.
			u = rand(1)[1]
		end
		#
		if  u <= 0.5
			β = 1. / (2. * u) ^ (1. / (method.η + 1.))
			E1[j] = 0.5 * ((1. + β) * method.population[i,j] + (1. - β) * method.population[i+1,j])
			E2[j] = 0.5 * ((1. - β) * method.population[i,j] + (1. + β) * method.population[i+1,j])
		else
			β = 1. / (2. * (1. - u)) ^( 1. / (method.η + 1.))
			E1[j] = 0.5 * (3. - β) * method.population[i,j] - (1. - β) * method.population[i+1,j]
			E2[j] = 0.5 * ((β - 1.) * method.population[i,j] + (3. - β) * method.population[i+1,j])
		end
	end
	# Copy the children in the population with overshoot tests
	method.population[i,1:n] = max.(min.(ub, E1), lb)
	method.population[i+1,1:n] = max.(min.(ub, E2), lb)
	# Update crossover type
	method.population[i,n+1], method.population[i+1,n+1] = 2., 2.
end
# PNX crossover
function PNX!(i, lb, ub, method::Clearing)
	# Parameter
	n = length(lb)
	#
	E1 = method.population[i,1:n] .+ abs.(method.population[i,1:n] .- method.population[i+1,1:n]) .* randn(n)
	E2 = method.population[i+1,1:n] .+ abs.(method.population[i,1:n] .- method.population[i+1,1:n]) .* randn(n)
	# Copy the children in the population with overshoot tests
	method.population[i,1:n] = max.(min.(ub, E1), lb)
	method.population[i+1,1:n] = max.(min.(ub, E2), lb)
	# Update crossover type
	method.population[i,n+1], method.population[i+1,n+1] = 3., 3.
end

### Mutations
function mutation!(lb, ub, method::Clearing)
	# Parameters
	n = length(lb)

    for i in (method.nind + 1) : (2 * method.nind)
		# Mutation poly
	    for j in 1:n
	        if rand(1)[1] < method.pm
	            u = rand(1)[1]
	            if u < 0.5
	                δ = 2. * u ^ (1. / (method.nm + 1.)) - 1.
				else
					δ = 1. - (2. * (1. - u)) ^ (1. / (method.nm + 1.))
	            end
				# Copy + overshoot
	            method.population[i,j] = min(max(method.population[i,j] + (ub[j]-lb[j]) * δ, lb[j]), ub[j])
	        end
	    end
		# Crossover gene mutation
	    if rand(1)[1] < method.pXgene
	       method.population[i,n+1] = unidrnd(3)
	    end
	end
end

### Assessment
function assessment!(f::Function, method::Clearing, options::Options)
    # Parameters
	n = size(method.population, 2) - 3
	# Asses the objective function
	if options.multithreads
	    Threads.@threads for i in (method.nind + 1) : (2 * method.nind)
	        method.population[i,n+3] = f(method.population[i,1:n])
	    end
	else
		for i in (method.nind + 1) : (2 * method.nind)
	        method.population[i,n+3] = f(method.population[i,1:n])
	    end
	end
end

### Clearing
function clearing!(lb, ub, method::Clearing)
   # Parameters
   n = length(lb)
   # Sort rows in ascending order from the fitness values
   method.population = sortrows!(method.population, n+3, order=true)
   #
   for i in 1 : (2 * method.nind)
      	if method.population[i,n+3] > -Inf
            ND = 1
            for j in (i + 1) : (2 * method.nind)
               	if method.population[j,n+3] > -Inf && norm((method.population[i,1:n] .- method.population[j,1:n]) ./ (ub .- lb), method.norm) < method.sigma
                 	ND < method.kapa ? ND = ND + 1 : method.population[j,n+3] = -Inf
           		end
        	end
      	end
   	end
end

### Density indexes to the niches
function set_density!(method::Clearing)
	# Parameters
    n = size(method.population, 2) - 3

	for i in 1 : (2 * method.nind)

		method.population[i,n+2] = 0.

		for j in 1:n

			method.population = sortrows!(method.population, j)

			δ = method.population[2 * method.nind,j] - method.population[1,j]

			if δ > 0.
				method.population[1, n+2] = Inf
				method.population[2 * method.nind, n+2] = Inf
				for k = 2 : (2 * method.nind - 1)
				 	method.population[k, n+2] < Inf ? method.population[k,n+2] = method.population[k,n+2] + (method.population[k+1,j] - method.population[k-1,j]) / δ : nothing
				end
			end
		end
	end
	# Sort in ascending order from the fitness and niches indexes
	# Pop = sortrows(Pop,[-(n+3) -(n+2)])
	method.population = sortrows!(method.population, n+3, order=true)
end

ext: jl


text:
module RPA
using Distributed: @distributed
using LinearAlgebra:diagm,eigen,Hermitian,cholesky,I
using SharedArrays

export seq2ind, hubbardvertexpp, vertexpp2ph,fermifunc,chiq,chiq0,chiqsimp, chikqm,projchi, kmesh2D,kqloc, kmesh

"""
seq2ind(x::Int,ind...)

find index (i,j,l) of `x`, ind=nspin,norb,nlat correspond to (spin,orbit,sublattice),i.e. x=i+(j-1)*nspin+(l-1)*nspin*norb
"""
function seq2ind(seq::Int,ind...)
    n=length(ind)
    ind1=cumprod(collect(ind))
    if ind1[end]<seq
        error("error:seq>product(ind)")
        return nothing
    end
    res = Int[]
    resm=Int[seq]
    for i=n:-1:1
        v1,v2=fldmod1(resm[end],ind1[i])
        push!(res,v1)
        push!(resm,v2)
    end
    push!(res,resm[end])
    return reverse(res[2:end])
end

"""
hubbardvertexpp(U::Float64,JH::Float64,U1::Float64,JP::Float64,nlat::Int=1,norb::Int=1,nspn::Int=2) -> Array{Float64,4}

return the particle-particle channel vertex Vpp_{ij,mn} i.e.
```math
(<ij|V|mn>-<ij|V|nm>)c⁺_{i} c⁺_{j} c_{n}  c_{m},
```
where i denotes the sublattice,orbital and spin degrees of freedom.
"""
function hubbardvertexpp(U::Float64,JH::Float64,U1::Float64,JP::Float64,nlat::Int=1,norb::Int=1,nspn::Int=2;exchg=true)
    n=nlat*norb*nspn
    v=zeros(n,n,n,n)
    for i1=1:n
        s1,o1,l1=seq2ind(i1,nspn,norb,nlat)
        for i2=1:n
            s2,o2,l2=seq2ind(i2,nspn,norb,nlat)
            for i3=1:n
                s3,o3,l3=seq2ind(i3,nspn,norb,nlat)
                for i4=1:n
                    s4,o4,l4=seq2ind(i4,nspn,norb,nlat)
                    if (l1==l2==l3==l4) && (s2==s4)&&(s1==s3)
                        if s2 != s1
                            if o1==o2==o3==o4
                                v[i4,i3,i2,i1]=U
                            elseif (o1==o2 !=o3==o4)
                                v[i4,i3,i2,i1]=JP
                            end
                        end
                        if (o2==o4 !=o1==o3)
                            v[i4,i3,i2,i1]=U1
                        end
                        if (o2==o3 !=o1==o4)
                            v[i4,i3,i2,i1]=JH
                        end
                    end
                end
            end
        end
    end
    if exchg
        return v-permutedims(v,[1,2,4,3])
    else
        return v
    end
end

"""
vertepp2ph(v::Array{Float64,4}) -> Array{Float64,4}

return the particle-hole vertex Vph_{ij,mn}=Vpp_{in,jm},Vph_{+-,-+}=Vpp_{++,--}
"""
vertexpp2ph(v::Array{Float64,4})= permutedims(v,[1,3,4,2])

"""
``1/N\\sum_{kk'q\\alpha\beta m n}[V^{ph}_{\\alpha\\beta,mn}(q)-V^{ph}_{\\alpha m,\\beta n}(q)]c^\\dagger_{k+q\\alpha}c_{k\\beta}c^\\dagger_{k'n}c_{k'+q m}``
"""
function vphq(q::Array{Float64,1})
    return nothing
end

"""
    fermifunc(e::T,temperature::T=1e-12,mu::T=0.0) where {T<:Real} -> Float64

Fermi distribution function.
"""
function fermifunc(e::T,temperature::T=1e-12,mu::T=0.0) where {T<:Real}
    if temperature > 1e-10
        f=(e-mu)/temperature
        if f>20
            f=0.0
        elseif f<-20
            f=1.0
        else
            f=1/(1+exp(f))
        end
    else
        f = e-mu
        if f>0.0
            f=0.0
        else
            f=1.0
        end
    end
    return f
end

"""
chikq0(model::Function,k::Array{Float64,1},q::Array{Float64,1},omega::Float64;eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...)

return chi0(k,q)_{ij,mn}= chi0(k,q)_{-+,+-}, chi0(k,q)_{12,34}== -<c^\\dagger_{k+q,2}c_{k ,1}c^\\dagger_{k,3}c_{k+q,4}>
`model` to get the matrix of tight-bind Hamiltonian,i.e.model(k::Array{Float64,1};kwargs...)
"""
function chikq0(model::Function,k::Array{Float64,1},q::Array{Float64,1},omega::Float64;eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...)
    Fk=eigen(model(k;kwargs...))
    Fkq=eigen(model(k+q;kwargs...))
    n=length(Fk.values)

    Ek=Fk.values
    Ekq=Fkq.values

    temp=diagm([-(fermifunc(Ek[i],tem,mu)-fermifunc(Ekq[j],tem,mu))/(omega+eta*im+Ek[i]-Ekq[j]) for i=1:n for j=1:n])
    u=kron(conj(Fk.vectors),Fkq.vectors)
    return u*temp*adjoint(u)
end

function chikq01(Ek::Vector{Float64},Ekq::Vector{Float64},uk::Matrix{ComplexF64},ukq::Matrix{ComplexF64},omega::Float64;eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0)
    n=length(Ek)
    temp=diagm([-(fermifunc(Ek[i],tem,mu)-fermifunc(Ekq[j],tem,mu))/(omega+eta*im+Ek[i]-Ekq[j]) for i=1:n for j=1:n])
    u=kron(conj(uk),ukq)
    return u*temp*adjoint(u)
end

"""
    chikq03(model::Function,k::Array{Float64,1},q::Array{Float64,1},omega::Float64=0.0;tem::Float64=1e-12,mu::Float64=0.0,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing,kwargs...) -> Tuple{Vector,Matrix,Matrix}
"""
function chikq03(model::Function,k::Array{Float64,1},q::Array{Float64,1},omega::Float64=0.0;tem::Float64=1e-12,mu::Float64=0.0,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing,kwargs...)
    @assert tem < 1e-10 "chikq03 error: now only support the zero temperature (<1e-10) $(tem). "
                Fk=eigen(model(k;kwargs...))
                Fkq=eigen(model(k+q;kwargs...))
                n=length(Fk.values)
                if bnd == nothing
                    bnd0 = 1:n
                else
                    bnd0 = bnd
                end
                Ek=Fk.values
                Ekq=Fkq.values
                ee1 = Float64[]
                ee2 = Float64[]
                lit1 = Int[]
                lit2 = Int[]
                one1=Float64[]
                one2=Float64[]

                for i = bnd0
                    for j = bnd0
                        ft = (fermifunc(Ek[i],tem,mu)-fermifunc(Ekq[j],tem,mu))
                        if abs(ft) >eps()
                            if Ek[i]-Ekq[j]<0.0
                                push!(lit1,(i-1)*n +j )
                                push!(ee1,(omega+Ek[i]-Ekq[j])/ft)
                                push!(one1,1.0)
                            else Ek[i]-Ekq[j]>0.0
                                    push!(lit2,(i-1)*n+j)
                                    push!(ee2,(omega+Ek[i]-Ekq[j])/ft)
                                    push!(one2,-1.0)
                            end
                        end
                    end
                end
                u=kron(conj(Fk.vectors),Fkq.vectors)
                ee=vcat(ee1,ee2)
        # @assert  length(ee1) == length(ee2) "chikq03 error: the number $(length(ee1)) of positive energy is not equal to that $(length(ee2)) of negative energy"

                uu = zeros(ComplexF64,n*n,length(ee))
                uu[:,1:length(ee1)]=u[:,lit1]
                uu[:,length(ee1)+1:end]=u[:,lit2]
                return ee, uu, vcat(one1,one2)
end
"""
chiq(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,3},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...) where T<:Number -> Tuple{Array{ComplexF64,4},Array{ComplexF64,4}}
chiq(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,2},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...) where T<:Number -> Tuple{Array{ComplexF64,4},Array{ComplexF64,4}}

Get the particle-hole susceptibilities χ⁰(ω,q) and χ(ω,q).
`model`::Function to get the matrix of tight-bind model, model(k::Array{Float64,1};kwargs...);
`km` store the k points of BZ;
`viq` store the q ponts;
`vph` is the bare particle-hole vertex; vph[ndim,ndim] just for hubbard interactions
`omegam` store the energy points;
`eta` is the magnitude of broaden;
`tem` is the temperature;
`mu` is the chemical potential.
`kwargs` is transfered to `model` function.
"""
function chiq(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,3},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...) where T<:Number
    nk=size(km,2)
    nq=size(viq,2)
    nw=length(omegam)
    ndim=size(vph,1)
    chi0 = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    chi = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    idmat=Matrix{Float64}(I,ndim,ndim)

    @sync @distributed for j=1:nq*nw #iq=1:nq
        iq,iw = fldmod1(j,nw)
        q = viq[:,iq]
        ω=omegam[iw]
        chi0[:,:,iw,iq] = @distributed (+) for i =1:nk
            1/nk*chikq0(model,km[:,i],q,ω;eta=eta,tem=tem,mu=mu,kwargs...)
        end
        chi[:,:,iw,iq] = chi0[:,:,iw,iq]*inv(idmat+vph[:,:,iq]*chi0[:,:,iw,iq])
    end
    return chi0,chi
end

function chiq(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,2},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...) where T<:Number
    nk=size(km,2)
    nq=size(viq,2)
    nw=length(omegam)
    ndim=size(vph,1)
    chi0 = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    chi = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    idmat=Matrix{Float64}(I,ndim,ndim)

    @sync @distributed for j=1:nq*nw
        iq,iw = fldmod1(j,nw)
        q = viq[:,iq]
        ω=omegam[iw]
        chi0[:,:,iw,iq] = @distributed (+) for i =1:nk
             1/nk*chikq0(model,km[:,i],q,ω;eta=eta,tem=tem,mu=mu,kwargs...)
        end
            chi[:,:,iw,iq] = chi0[:,:,iw,iq]*inv(idmat+vph*chi0[:,:,iw,iq])
    end

    return chi0,chi
end

"""
    kmesh2D(ak::Vector{Float64},bk::Vector{Float64},nab::Tuple{Int,Int},orign::Vector{Float64}=Float64[0.0,0.0]) ->Array{Float64,2}
Construct the k mesh in 2-D BZ. size(result) =(2,na,nb)
`ak`: one of two reciprocal lattice vectors;
`bk`: another reciprocal lattice vector;
`nab`: (na,nb), na (nb) is the number of k points along the `ak` (`bk`) direction;
`orign`:the location of orignal point.
"""
function kmesh2D(ak::Vector{Float64},bk::Vector{Float64},nab::Vector{Int},orign::Vector{Float64}=Float64[0.0,0.0])
    na=nab[1]
    nb=nab[2]
    res = zeros(Float64,2,na,nb)
    for i=1:nb
        for j=1:na
            res[:,j,i] = ak/na*(j-1)+bk/nb*(i-1)+orign
        end
    end
    return res
end
"""
    kmesh(abck::Matrix{Float64},,nab::Vector{Int},orign::Vector{Float64}) ->Array{Float64,Ndim}
Construct the k mesh in BZ. size(result) =(dim,nab...)
`abck`: store reciprocal lattice vectors;
`nab`: store the number of k ponts along the corresponding direction. e.g. (na,nb,nc), na (nb,nc) is the number of k points along the `ak` (`bk`,`ck`) direction;
`orign`:the location of orignal point.
"""
function kmesh(abck::Matrix{Float64},nab::Vector{Int},orign::Vector{Float64})
    @assert size(abck,2) == size(nab,1) "kmesh error: size(abck,2) != size(nab,1)"
    @assert size(abck,1) == size(orign,1) "kmesh error: size(abck,1) != size(orign,1)"
    n = size(nab,1)
    if n == 2
        res = zeros(Float64,2,nab[1],nab[2])
        for i=1:nab[2]
            for j=1:nab[1]
                res[:,j,i] = abck[:,1]/nab[1]*(j-1)+abck[:,2]/nab[2]*(i-1)+orign
            end
        end
    elseif n==3
        res = zeros(Float64,3,nab[1],nab[2],nab[3])
        for k = 1:nab[3]
            for i=1:nab[2]
                for j=1:nab[1]
                    res[:,j,i,k] = abck[:,1]/nab[1]*(j-1)+abck[:,2]/nab[2]*(i-1)+ abck[:,3]/nab[3]*(k-1)+orign
                end
            end
        end
    elseif n==1
        res = zeros(Float64,1,nab[1])
        for i = 0:nab[1]
            res[:,i] = abck[:,1]/nab[1]*i + orign
        end
    end
    return res
end
"""
    kqloc(kq::Vector{Float64},dab::Matrix{Float64},nab::Vector{Int},orign::Vector{Float64})->Tuple{Int,Vector{Int}}

To find the location k+q point in the array (result of `kmesh`or `kmesh2D` function ). The first element Tuple{Int,Vector{Int}} is the linear index,the second element of Tuple is the cartesian index.
"""
function kqloc(kq::Vector{Float64},dab::Matrix{Float64},nab::Vector{Int},orign::Vector{Float64})
    @assert length(kq) == size(dab,1) "kqloc error: length(kq)$(length(kq)) == size(dab,1)$(size(dab,1))"
    @assert size(dab,2) == size(nab,1) "kqloc error: size(dab,2) [$(size(dab,2))] == size(nab,1) [$(size(nab,1))]"
    x = dab\(kq-orign)
    intx = round.(Int,x)
    res = map(mod,intx,nab)
    res = res.+1
    nabc = cumprod(nab)
    n = length(nabc)
    if n==2
        (res0=(res[2]-1)*nabc[1] + res[1])
    elseif n==3
        (res0=(res[3]-1)*nabc[2]+(res[2]-1)*nabc[1] + res[1])
    elseif n==1
         (res0 == res[1])
    else
        @assert n>=4 "kqloc error: dimension of bz should be smaller than 4."
    end
    return  res0,res
end
"""
    chiq(eigenvc::Array{ComplexF64,3},eigenval::Array{Float64,2},km::Array{Float64,2},nab::Vector{Int},viq::Array{Float64,2},vph::Array{T,2},omegam::Vector{Float64};eta::Float64=0.01,orign::Vector{Float64}=[0.0,0.0],tem::Float64=1e-12,mu::Float64=0.0) where T<:Number -> Tuple{Array{ComplexF64,4},Array{ComplexF64,4}}

"""
function chiq(eigenvc::Array{ComplexF64,3},eigenval::Array{Float64,2},km::Array{Float64,2},nab::Vector{Int},viq::Array{Float64,2},vph::Array{T,2},omegam::Vector{Float64};eta::Float64=0.01,orign::Vector{Float64}=[0.0,0.0],tem::Float64=1e-12,mu::Float64=0.0) where T<:Number
    nk=size(km,2)
    nq=size(viq,2)
    nw=length(omegam)
    ndim=size(vph,1)
    chi0 = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    chi = SharedArray{ComplexF64}(ndim,ndim,nw,nq)
    idmat=Matrix{Float64}(I,ndim,ndim)
    ndim0 = size(nab,1)
    dab = zeros(Float64,size(km,1),ndim0)
    @assert size(orign,1) == ndim0 "chiq error: size(orign,1)$(size(orign,1)) == size(nab,1)$(ndim0)"

    dab[:,1] = km[:,2] - km[:,1]
    ndim0 <=3 && ndim0>=2 &&  (dab[:,2] = km[:,1+nab[1]] - km[:,1])
    ndim0 == 3 && (dab[:,3] = km[:,1+nab[1]*nab[2]] - km[:,1])

    @sync @distributed for j=1:nq*nw
        iq,iw = fldmod1(j,nw)
        q = viq[:,iq]
        ω=omegam[iw]

        chi0[:,:,iw,iq] = @distributed (+) for i =1:nk
            kq=km[:,i]+q
            ikq,=kqloc(kq,dab,nab,orign)
            1/nk*chikq01(eigenval[:,i],eigenval[:,ikq],eigenvc[:,:,i],eigenvc[:,:,ikq],ω;eta=eta,tem=tem,mu=mu)
        end
            chi[:,:,iw,iq] = chi0[:,:,iw,iq]*inv(idmat+vph*chi0[:,:,iw,iq])
    end

    return chi0,chi
end

"""
chiq0(model::Function,km::Array{Float64,2},viq::Array{Float64,2},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...) -> Array{Float64,4}

Get the particle-hole susceptibilities χ⁰(ω,q).
`model`::Function to get the matrix of tight-bind model;
`km` store the k points of BZ;
`viq` store the q ponts;
`omegam` store the energy points;
`eta` is the magnitude of broaden;
`tem` is the temperature;
`mu` is the chemical potential.
'kwargs' is transfered to `model` function
"""
function chiq0(model::Function,km::Array{Float64,2},viq::Array{Float64,2},omegam::Vector{Float64};eta::Float64=0.01,tem::Float64=1e-12,mu::Float64=0.0,kwargs...)
                    nk=size(km,2)
                    nq=size(viq,2)
                    nw=length(omegam)
                    ndim=size(vph,1)
                    chi0=zeros(ComplexF64,ndim,ndim,nw,nq)
                    #idmat=Matrix{Float64}(I,ndim,ndim)
                    for iq=1:nq
                        for iw=1:nw
                            for i =1:nk
                                chi0[:,:,iw,iq] += 1/nk*chikq0(model,km[:,i],viq[:,iq],omegam[iw];eta=eta,tem=tem,mu=mu,kwargs...)
                            end
                        end
                    end
                    return chi0
                end
"""
chiqsimp(vph::Array{T,3},chi0::Array{ComplexF64,4}) where T<:Number -> Array{ComplexF64,4}
chiqsimp(vph::Array{T,2},chi0::Array{ComplexF64,4}) where T<:Number -> Array{ComplexF64,4}

Get the susceptibility χ_{αβ}(ω,q).vph::Array{T,2} for hubbard interactions,
"""
function chiqsimp(vph::Array{T,3},chi0::Array{ComplexF64,4}) where T<:Number
                    nq=size(chi0,4)
                    nw=size(chi0,3)
                    ndim=size(vph,1)
                    chi=zeros(ComplexF64,ndim,ndim,nw,nq)
                    idmat=Matrix{Float64}(I,ndim,ndim)
                    for iq=1:nq
                        for iw=1:nw
                            chi[:,:,iw,iq] = chi0[:,:,iw,iq]*inv(idmat+vph[:,:,iq]*chi0[:,:,iw,iq])
                        end
                    end
                    return chi
                end
 function chiqsimp(vph::Array{T,2},chi0::Array{ComplexF64,4}) where T<:Number
                    nq=size(chi0,4)
                    nw=size(chi0,3)
                    ndim=size(vph,1)
                    chi=zeros(ComplexF64,ndim,ndim,nw,nq)
                    idmat=Matrix{Float64}(I,ndim,ndim)
                    for iq=1:nq
                        for iw=1:nw
                            chi[:,:,iw,iq] = chi0[:,:,iw,iq]*inv(idmat+vph[:,:]*chi0[:,:,iw,iq])
                        end
                    end
                    return chi
                end

"""
    chikqm(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,5};tem::Float64=1e-12,mu::Float64=0.0,isval::Bool=true,η::Float64=1e-6,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing, kwargs...) where T<:Number -> Tuple{Array{Array{ComplexF64,1},1},Array{Matrix{ComplexF64},1},Array{Matrix{ComplexF64},1}}
    chikqm(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,2};tem::Float64=1e-12,mu::Float64=0.0,isval::Bool=true,η::Float64=1e-6,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing, kwargs...) where T<:Number -> Tuple{Array{Array{ComplexF64,1},1},Array{Matrix{ComplexF64},1},Array{Matrix{ComplexF64},1}}

Get the eigenvalues, eigenvectors, and unitary transformation (from orbital to band) of particle-hole susceptibilities χ_{αβ}(ω,k1,k2,q). Now only the zero-temperature case is supported.
`model`::Function to get the matrix of tight-bind model,model(k::Array{Float64,1};kwargs...);
`km` store the k points of BZ;
`viq` store the q ponts;
`vph` store the particle-hole vertex,e.g. vph[ndim,ndim,nk,nk,nq] where sqrt(nidm) is the number of degrees of freedom in the unit cell,nk=size(km,2),nq=size(viq,2). or vph[ndim,ndim] for hubbard interactions;
`eta` is the magnitude of broaden;
`tem` is the temperature;
`mu` is the chemical potential.
`isval` : only need the eigenvalues ( isval=true,zero temperature ), isval=false denotes that the cholesky method is used.
`η`:small number to advoid the semipositive Hamiltonian,i.e. Hamiltonian+diagm([η,η,...])
`bnd`:select the bands to calculate the χ₀
`kwargs` store the keys which are transfered to `model` function.
"""
function chikqm(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,5};tem::Float64=1e-12,mu::Float64=0.0,isval::Bool=true,η::Float64=1e-6,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing, kwargs...) where T<:Number
    nk=size(km,2)
    nq=size(viq,2)
    #ndim = size(vph,1)
    val = Array{ComplexF64,1}[]
    vec = Array{ComplexF64,2}[]
    f2c = Array{ComplexF64,2}[]

    for iq=1:nq
       # res2 = Int[]
        chi0inv = Float64[]
        um = Array{ComplexF64,2}[]
        list=Int[1]
        g1=Float64[]
        for i =1:nk
            temp0,tempu,g0=chikq03(model,km[:,i],viq[:,iq],0.0;tem=tem,mu=mu,bnd=bnd,kwargs...)
            push!(um,tempu)
            push!(chi0inv,temp0...)
            push!(g1,g0...)
            n0 =size(tempu,2)
            push!(list,n0)
        end
        list0 = cumsum(list)
        wk = zeros(ComplexF64,list0[end]-1,list0[end]-1)
        for i =1:nk
            for j=1:nk
                wk[list0[j]:list0[j+1]-1,list0[i]:list0[i+1]-1]=transpose(conj(um[j]))*vph[:,:,j,i,iq]*um[i]/nk
            end
        end
         temp1=wk

        temp3=Hermitian(-diagm(chi0inv)*1.0+1.0*temp1,:U)

        g=diagm(g1)
        if isval
            temp2=eigen(temp3*g)
        else
            nd = size(temp3,1)
            ck=cholesky(temp3+diagm([η for i=1:nd]), Val(false); check = true)
            temp2=eigen(ck.U*g*ck.L)
            temp4=inv(ck.U)*temp2.vectors*diagm(sqrt.(abs.(temp2.values)))
            push!(vec,temp4)
        end
        push!(val,temp2.values)
        push!(f2c,hcat(um...))

    end
    return val,vec,f2c
end

function chikqm(model::Function,km::Array{Float64,2},viq::Array{Float64,2},vph::Array{T,2};tem::Float64=1e-12,mu::Float64=0.0,isval::Bool=true,η::Float64=1e-6,bnd::Union{UnitRange{Int},StepRange{Int,Int},Vector{Int},Nothing}=nothing, kwargs...) where T<:Number
    nk=size(km,2)
    nq=size(viq,2)
    #ndim = size(vph,1)
    val = Array{ComplexF64,1}[]
    vec = Array{ComplexF64,2}[]
    f2c = Array{ComplexF64,2}[]

    for iq=1:nq
       # res2 = Int[]
        chi0inv = Float64[]
        um = Array{ComplexF64,2}[]
        list=Int[1]
        g1=Float64[]
        for i =1:nk
            temp0,tempu,g0=chikq03(model,km[:,i],viq[:,iq],0.0;tem=tem,mu=mu,bnd=bnd,kwargs...)
            push!(um,tempu)
            push!(chi0inv,temp0...)
            push!(g1,g0...)
            n0 =size(tempu,2)
            push!(list,n0)
        end
        list0 = cumsum(list)
        wk = zeros(ComplexF64,list0[end]-1,list0[end]-1)
        for i =1:nk
            for j=1:nk
                wk[list0[j]:list0[j+1]-1,list0[i]:list0[i+1]-1]=transpose(conj(um[j]))*vph*um[i]/nk
            end
        end
         temp1=wk

        temp3=Hermitian(-diagm(chi0inv)*1.0+1.0*temp1,:U)

# println(iq," g1test:",g1)
        g=diagm(g1)
        if isval
            temp2=eigen(temp3*g)
        else
            nd = size(temp3,1)
            ck=cholesky(temp3+diagm([η for i=1:nd]), Val(false); check = true)
            temp2=eigen(ck.U*g*ck.L)
            temp4=inv(ck.U)*temp2.vectors*diagm(sqrt.(abs.(temp2.values)))
            push!(vec,temp4)
        end
        push!(val,temp2.values)
        push!(f2c,hcat(um...))

    end
    return val,vec,f2c
end

"""
    projchi(val::Array{Array{T,1},1},vec::Array{Array{T1,2},1},f2c::Array{ Array{ComplexF64,2},1},omegam::Vector{Float64},eta::Float64=1e-2) where {T<:Number,T1<:Number}->Array{ComplexF64,4}

Get the ``\\chi_{ij,nm}(\\omega,q)```. The eigenvalues,eigenvectors, and unitary matrix are obtained by method `chikqm`.
"""
function projchi(val::Array{Array{T,1},1},vec::Array{Array{T1,2},1},f2c::Array{ Array{ComplexF64,2},1},omegam::Vector{Float64};eta::Float64=1e-2) where {T<:Number,T1<:Number}
                    nq = length(val)
                    nw = length(omegam)
                    ndim,n=size(f2c[1])

                    chi=zeros(ComplexF64,ndim,ndim,nw,nq)
                    for iq =1:nq
                        for iw = 1:nw
                            tm = vec[iq]
                            tmu = f2c[iq]*tm
                            chi[:,:,iw,iq]=tmu*diagm([1/(-omegam[iw]-im*eta+val[iq][i]) for i=1:n])*inv(diagm(sign.(val[iq])))*adjoint(tmu)
                        end
                    end
                    return chi
                end

 """
    projchiim(val::Array{Array{T,1},1},vec::Array{Array{T1,2},1},f2c::Array{ Array{ComplexF64,2},1},omegam::Vector{Float64},eta::Float64=1e-2) where {T<:Number,T1<:Number}->Array{ComplexF64,4}

Get the image of ``\\chi_{ij,nm}(\\omega,q)```. The eigenvalues,eigenvectors, and unitary matrix are obtained by method `chikqm`.
"""
function projchiim(val::Array{Array{T,1},1},vec::Array{Array{T1,2},1},f2c::Array{ Array{ComplexF64,2},1},omegam::Vector{Float64};eta::Float64=1e-2) where {T<:Number,T1<:Number}
                    nq = length(val)
                    nw = length(omegam)
                    ndim=size(f2c[1],1)
                    chi=zeros(ComplexF64,ndim,ndim,nw,nq)

                    for iq =1:nq
                        for iw = 1:nw
                            for i =1:length(val[iq])
                                temp=f2c[iq]*vec[iq][:,i]
                                chi[:,:,iw,iq] +=temp*eta/((omegam[iw]-val[iq][i])^2+eta^2)*adjoint(temp)
                            end
                        end
                    end
                    return chi
                end

end # module


text:
include(joinpath(BASE_DIR, "test", "data_5bus_pu.jl"))
include(joinpath(BASE_DIR, "test", "data_14bus_pu.jl"))

checksys = false

@testset "Test System constructors from .jl files" begin
    tPowerSystem = System(nothing)
    nodes_5 = nodes5()
    nodes_14 = nodes14()

    for i in nodes_5
          nodes_5[i].angle = deg2rad(nodes_5[i].angle)
    end

    # Components with forecasts cannot be added to multiple systems, so clear them on each
    # test.

    sys5 = System(nodes_5, thermal_generators5(nodes_5), loads5(nodes_5), nothing, nothing,  100.0, nothing,
                  nothing; runchecks = checksys)
    clear_components!(sys5)

    sys5b = System(nodes_5, thermal_generators5(nodes_5), loads5(nodes_5), nothing, battery5(nodes_5),  100.0, nothing,
                   nothing; runchecks = checksys)
    clear_components!(sys5b)

    # GitHub issue #234 - fix forecasts5 in data file, use new format
    #_sys5b = PowerSystems._System(nodes_5, thermal_generators5(nodes_5), loads5(nodes_5), nothing, battery5(nodes_5),
    #                              100.0, forecasts5, nothing, nothing)
    #sys5b = System(_sys5b)

    sys5bh = System(nodes_5, vcat(thermal_generators5(nodes_5), hydro_generators5(nodes_5)), loads5(nodes_5), branches5(nodes_5),
                    battery5(nodes_5),  100.0, nothing, nothing; runchecks = checksys)
    clear_components!(sys5bh)

    sys5bh = System(; buses=nodes_5,
                    generators=vcat(thermal_generators5(nodes_5), hydro_generators5(nodes_5)),
                    loads=loads5(nodes_5),
                    branches=branches5(nodes_5),
                    storage=battery5(nodes_5),
                    basepower=100.0,
                    services=nothing,
                    annex=nothing,
                    runchecks = checksys)
    clear_components!(sys5bh)

    # Test Data for 14 Bus

    # GitHub issue #234 - fix forecasts5 in data file, use new format
    #_sys14 = PowerSystems._System(nodes_14, thermal_generators14, loads14, nothing, nothing,
    #                            100.0, Dict{Symbol,Vector{<:Forecast}}(),nothing,nothing)
    #sys14 = System(_sys14)

    for i in nodes_14
          nodes_14[i].angle = deg2rad(nodes_14[i].angle)
    end

    sys14b = PowerSystems.System(nodes_14, thermal_generators14(nodes_14), loads14(nodes_14), nothing,
                                 battery14(nodes_14), 100.0, nothing, nothing; runchecks = checksys)
    clear_components!(sys14b)
    sys14b = PowerSystems.System(nodes_14, thermal_generators14(nodes_14), loads14(nodes_14), branches14(nodes_14),
                                 battery14(nodes_14), 100.0, nothing, nothing; runchecks = checksys)
    clear_components!(sys14b)
end

@testset "Test System constructor from Matpower" begin
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))
end


text:
##  Fast Fourier transform of a real-valued signal
##  + simple power spectrum
##  + Hilbert transform
##  Author Bernd Blasius (bernd.blasius@gmail.com)

# this files will need the Julia Package FFTW.jl

function power(x,fs)
    # simple power spectrum
    # x: time series
    # f: sampling frequency
    n = length(x)
    xbar = fftshift(fft(x))
    f = (-n/2:n/2-1)*(fs/n)    # zero-centered frequency range
    power= abs.(xbar).^2/n     # zero-centered power

    f, power
end

function power_real(x,fs)
    # simple power spectrum using rfft
    # x: time series
    # f: sampling frequency
    n = length(x)
    xbar = rfft(x)
    n1 = length(xbar)
    f = (0:n1-1) * (fs/n)     #frequency range
    power = abs.(xbar).^2/n    #power

    f, power
end

function hilbert(x,fs)
    # calculate the Hilbert transform of a real valued signal
    n = length(x)
    xbar = rfft(x)
    irfft(im*xbar,n) # Hilbert transform
end

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
const qnm = PyNULL()

function __init__()
    copy!(qnm, pyimport("qnm"))
    qnm.download_data()
end

function qnmfunctionnew(s,l,m,n,a; qnm=qnm)
    grav_freq = qnm.modes_cache(s=s,l=l,m=m,n=n)
    ω, Alm, Cllʼ = grav_freq(a=a)
    qnmfunction(Custom; s=s,l=l,m=m,n=n,a=a,ω=ω,Alm=Alm,Cllʼ=Cllʼ)
end

max_stars_repo_name: Potatoasad/KerrQNMShifts

max_issues_repo_name: Potatoasad/KerrQNMShifts

max_forks_repo_name: Potatoasad/KerrQNMShifts


text:
module L1L2

    using DataFrames
    using LinearAlgebra
    using ..DataMod, ..ManualModelMod

    export l1l2

    # Soft thresholding
    function S(x::Float64, λ::Float64)::Float64
        if x >= λ
            return x - λ
        elseif x <= -λ
            return x + λ
        else
            return 0
        end
    end
    S(xs::Vector{Float64}, λ::Float64)::Vector{Float64} = [S(x, λ) for x in xs]
    S(xs::Matrix{Float64}, λ::Float64)::Vector{Float64} = S(vec(xs), λ)

    # Calculate next weight vector
    function nextβ(β::Vector{Float64}, τplus2ϵλ::Float64, λγ::Float64,
                   τIminusXᵀX::Matrix{Float64}, XᵀY::Vector{Float64})::Vector{Float64}
        return S(τIminusXᵀX * β + XᵀY, λγ) / τplus2ϵλ
    end

    # Calculate max iteration
    function lmax(β₁::Vector{Float64}, β₀::Vector{Float64}, κ::Float64,
                  κ₀::Float64, λ::Float64, ϵ::Float64, ξ::Float64)::Int
        # numnum = norm(β₁ - β₀) * (κ + κ₀ + 4 * ϵ * λ)
        # numdenom = (2 * κ₀ + 4 * ϵ * λ) * ξ
        # denomnum = κ + κ₀ + 4 * ϵ * λ
        # denomdenom = κ - κ₀ # THIS GOES WRONG IF κ == κ₀
        # num = log(numnum / numdenom)
        # denom = log(denomnum / denomdenom)
        # return round(Int, num / denom + 1, RoundDown)
        return 10000
    end

    function weight_loop(data::Data, λ::Float64, γ::Float64, ϵ::Float64, ξ::Float64)
        X = hcat(ones(Float64, size(data.xmat, 1)), data.xmat)
        Ys = Vector{Float64}[data.df[!, y] for y in data.ys]
        XᵀX = X'X
        τ = norm(XᵀX)
        τplus2ϵλ = τ + 2 * ϵ * λ
        λγ = λ * γ
        τIminusXᵀX = τ * I - XᵀX
        XᵀYs = [X'Y for Y in Ys]
        β₀ = zeros(length(data.xs) + 1)
        βs = [nextβ(β₀, τplus2ϵλ, λγ, τIminusXᵀX, XᵀY) for XᵀY in XᵀYs]
        κ, κ₀ = τ, τ # should change
        max_iter = maximum([lmax(β, β₀, κ, κ₀, λ, ϵ, ξ) for β in βs])::Int
        @inbounds for i in 1:max_iter
            βs = Vector{Float64}[nextβ(βs[y], τplus2ϵλ, λγ, τIminusXᵀX, XᵀYs[y]) for y in 1:length(data.ys)]
        end
        return βs
    end

    function l1l2(data::Data; λ::Float64 = 0.5, γ::Float64 = 1.0, ϵ::Float64 = 1.0,
                  ξ::Float64 = 0.001)::Vector{ManualModel}
        # λ, γ, ϵ, ξ = (0.5, 1.0, 1.0, 0.001)
        weights = weight_loop(data, λ, γ, ϵ, ξ)
        return ManualModel.(weights, true, data.ys, data)
    end

end

using .L1L2
export l1l2

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
# This software is licensed under the MIT License (MIT).

import Pkg

srcdir = @__DIR__

Pkg.activate(srcdir)

# Need to have PackageCompiler installed
Pkg.add("PackageCompiler", preserve = Pkg.PRESERVE_ALL)

prj_file = Pkg.project().path
prj_dir = dirname(prj_file)

@info("Building custom system image in project \"$prj_dir\".")

@info("Instantiating project.")
Pkg.instantiate()

@info("Precompiling packages.")
Pkg.precompile()

excluded = [
  "PyCall",
  "SymPy",
  "StatsPlots",
  "ImageMagick",
  "StatsBase",
  "PackageCompiler",
  "Ipopt",
  "SpecialFunctions",
]

prj_sysimage = get(Pkg.TOML.parsefile(prj_file), "deps", Dict{String,Any}())
sysimg_pkgs = filter(x -> x ∉ excluded, sort(collect(keys(prj_sysimage))))
@info("Package to include in system image: $(join(sysimg_pkgs, " "))")

@info("Building system image.")

import PackageCompiler, Libdl

sysimage_path = joinpath(srcdir, "JuliaSysimage." * Libdl.dlext)

PackageCompiler.create_sysimage(
    Symbol.(sysimg_pkgs),
    sysimage_path = sysimage_path,
    precompile_execution_file = joinpath(srcdir, "precompile_exec.jl"),
    cpu_target = PackageCompiler.default_app_cpu_target(),
    replace_default = false
)

@info "Created custom Julia system image"

max_stars_repo_name: prismia-chat/flux-binder

max_issues_repo_name: prismia-chat/flux-binder

max_forks_repo_name: prismia-chat/flux-binder


text:
abstract type StorageType{T<:Real} end

Base.eltype(::StorageType{T}) where {T<:Real} = T

is_of_storage_type(::Any, ::StorageType) = false

is_of_storage_type(::T, ::StorageType{T}) where {T<:Real} = true

is_of_storage_type(x::AbstractArray{<:Real}, s::StorageType) = false

is_of_storage_type(x::AbstractArray, s::StorageType) = all(is_of_storage_type.(x, Ref(s)))

# A Tuple of objects are of a particular storage type if each element of the tuple is of
# that storage type.
function is_of_storage_type(xs::Union{Tuple, NamedTuple}, s::StorageType)
    return all(map(x -> is_of_storage_type(x, s), xs))
end

is_of_storage_type(x::Gaussian, s::StorageType) = is_of_storage_type((x.m, x.P), s)

is_of_storage_type(::Nothing, ::StorageType) = true

#
# Indicate to use `StaticArray`s to store model parameters.
#

struct SArrayStorage{T<:Real} <: StorageType{T} end

SArrayStorage(T) = SArrayStorage{T}()

storage_type(::SArray{<:Any, T}) where {T} = SArrayStorage(T)

is_of_storage_type(::SArray{<:Any, T}, ::SArrayStorage{T}) where {T<:Real} = true

#
# Indicate to use `Array`s to store model parameters.
#

struct ArrayStorage{T<:Real} <: StorageType{T} end

ArrayStorage(T) = ArrayStorage{T}()

storage_type(::Array{T}) where {T} = ArrayStorage(T)

is_of_storage_type(::Array{T}, ::ArrayStorage{T}) where {T<:Real} = true

#
# Working with scalars.
#

struct ScalarStorage{T<:Real} <: StorageType{T} end

ScalarStorage(T) = ScalarStorage{T}()

storage_type(::T) where {T<:Real} = T

is_of_storage_type(::T, ::ScalarStorage{T}) where {T<:Real} = true


text:
include("Scanner.jl")

struct SyntaxException <: Exception
  msg::String
end

@enum MPType begin
  MInt
  MReal
  MBool
  MString
  MIntArray
  MRealArray
  MBoolArray
  MStringArray
  MIntRef
  MRealRef
  MBoolRef
  MStringRef
  MNothing
  MError
  MPassed
  MUndefined
end

ref_types = Vector{MPType}([
  MIntRef,
  MRealRef,
  MBoolRef,
  MStringRef
])

array_types = Vector{MPType}([
  MIntArray,
  MRealArray,
  MBoolArray,
  MStringArray
])

scalar_types = Vector{MPType}([
  MInt,
  MReal,
  MBool,
  MString
])

token_class_to_scalar_type = Dict{TokenClass,MPType}(
  kw_int => MInt,
  kw_real => MReal,
  kw_boolean => MBool,
  kw_string => MString
)

token_class_to_ref_type = Dict{TokenClass,MPType}(
  kw_int => MIntRef,
  kw_real => MRealRef,
  kw_boolean => MBoolRef,
  kw_string => MStringRef
)

scalar_to_ref_type = Dict{MPType,MPType}(
  MInt => MIntRef,
  MReal => MRealRef,
  MBool => MBoolRef,
  MString => MStringRef,
  MIntRef => MIntRef,
  MRealRef => MRealRef,
  MBoolRef => MBoolRef,
  MStringRef => MStringRef,
  MIntArray => MIntArray,
  MRealArray => MRealArray,
  MBoolArray => MBoolArray,
  MStringArray => MStringArray
)

scalar_to_array_type = Dict{MPType,MPType}(
  MInt => MIntArray,
  MReal => MRealArray,
  MBool => MBoolArray,
  MString => MStringArray
)

mutable struct ImplicitParameter
  name::String
  type::MPType
  unique_id::String
end

mutable struct SymTableEntry
  var_name::String
  var_type::MPType
  val_identifier::String
  scope_level::Int
  is_implicit::Bool
end
SymTableEntry(var_name::String, var_type::MPType, val_identifier::String, scope_level::Int) =
  SymTableEntry(var_name::String, var_type::MPType, val_identifier::String, scope_level::Int, false)

abstract type Node end
# Pseudonodes are nodelike structs that assist in parsing
# but are not included in the abstract syntax tree (AST).
abstract type Pseudonode <: Node end
abstract type Statement <: Node end
abstract type Value <: Node end
abstract type Factor <: Value end
abstract type If <: Statement end

mutable struct BinaryOperation <: Value
  op::TokenClass
  left::Value
  right::Value
  line::Int
  type::MPType
end
BinaryOperation(op, left, right, line) = BinaryOperation(op, left, right, line, MUndefined)

mutable struct UnaryOperation <: Value
  op::TokenClass
  operand::Value
  line::Int
  type::MPType
end
UnaryOperation(op, operand, line) = UnaryOperation(op, operand, line, MUndefined)

mutable struct ImmediateInt <: Value
  value::Int
  type::MPType
end
ImmediateInt(value::Int) = ImmediateInt(value, MUndefined)

mutable struct ImmediateReal <: Value
  value::Float64
  type::MPType
end
ImmediateReal(value::Float64) = ImmediateReal(value, MUndefined)

mutable struct ImmediateString <: Value
  value::String
  type::MPType
end
ImmediateString(value::String) = ImmediateString(value, MUndefined)

mutable struct ImmediateBool <: Value
  value::Bool
  type::MPType
end
ImmediateBool(value::Bool) = ImmediateBool(value, MUndefined)

mutable struct Block <: Statement
  statements::Vector{Statement}
  is_subroutine_block::Bool
  subroutine::Union{Node,Nothing}
  line::Int
  type::MPType
end
Block(statements::Vector{Statement}, line::Int) = Block(statements, false, nothing, line, MUndefined)

mutable struct TypeOfVarOrValue <: Node
  scalar_type::MPType
  true_type::MPType
  size::Value
  line::Int
  type::MPType
end
TypeOfVarOrValue(scalar_type::MPType, true_type::MPType, size::Value, line::Int) =
  TypeOfVarOrValue(scalar_type, true_type, size, line, MUndefined)

function is_array_type(t::TypeOfVarOrValue)
  return t.size isa ImmediateInt && t.size.value == 1
end

mutable struct Parameter <: Node
  name::String
  scalar_type::MPType
  is_var_par::Bool
  is_array::Bool
  size::Value       # Array size can only be determined at run time
  line::Int
  type::MPType
end
Parameter(name::String,
  scalar_type::MPType,
  is_var_par::Bool,
  is_array::Bool,
  size::Value,
  line::Int) =
  Parameter(name,
    scalar_type,
    is_var_par,
    is_array,
    size,
    line,
    MUndefined)

# Function or procedure.
mutable struct Subroutine <: Node
  name::String
  params::Vector{Parameter}
  ret_type::TypeOfVarOrValue
  body::Block
  calls_subroutines::Vector{String}
  llvm_function_name::String
  is_called::Bool   # No codegen for subroutines that are never called
  line::Int
  type::MPType
end
Subroutine(name, params, ret_type, body, line) =
  Subroutine(
    name,
    params,
    ret_type,
    body,
    Vector{String}(),
    "",
    false,
    line,
    MUndefined)

mutable struct Definitions <: Node
  defs::Vector{Subroutine}
  line::Int
  type::MPType
  end
Definitions(defs, line) = Definitions(defs, line, MUndefined)

mutable struct Program <: Node
  definitions::Definitions
  main::Block
  line::Int
  type::MPType
  end
Program(definitions, main, line) = Program(definitions, main, line, MUndefined)

mutable struct CallStatement <: Statement
  identifier::String
  arguments::Vector{Value}
  implicit_params::Vector{SymTableEntry}
  subroutine::Union{Subroutine,Nothing}
  call_id::Int
  fingerprint::String
  scope_level::Int
  line::Int
  type::MPType
end
CallStatement(identifier, arguments, call_id, line) =
  CallStatement(identifier, arguments, Vector{SymTableEntry}(), nothing, call_id, "", 0, line, MUndefined)

mutable struct Declaration <: Statement
  var_type::TypeOfVarOrValue
  names::Vector{Token}
  unique_ids::Vector{String}
  line::Int
  type::MPType
end
Declaration(var_type, names, line) = Declaration(var_type, names, Vector{String}(), line, MUndefined)

mutable struct Variable <: Node
  identifier::String
  is_array_access::Bool
  array_index::Value
end

mutable struct Assignment <: Statement
  variable_name::String
  var_unique_id::String
  var_type::MPType
  is_array_access::Bool
  array_index::Value  # Array index. Can use ImmediateInt(-1) for non-array access assignments.
  value::Value
  line::Int
  type::MPType
end
Assignment(variable_name::String, is_array_access::Bool, index::Value, value::Value, line::Int) =
  Assignment(variable_name, "", MUndefined, is_array_access, index, value, line, MUndefined)

mutable struct IfThen <: If
  condition::Value
  then_stmt::Statement
  line::Int
  type::MPType
end
IfThen(condition, then_stmt, line) = IfThen(condition, then_stmt, line, MUndefined)

mutable struct IfThenElse <: If
  condition::Value
  then_stmt::Statement
  else_stmt::Statement
  line::Int
  type::MPType
end
IfThenElse(condition, then_stmt, else_stmt, line) = IfThenElse(condition, then_stmt, else_stmt, line, MUndefined)

mutable struct While <: Statement
  condition::Value
  do_stmt::Statement
  line::Int
  type::MPType
end
While(condition, do_stmt, line) = While(condition, do_stmt, line, MUndefined)

mutable struct Read <: Statement
  variables::Vector{Variable}
  var_types::Vector{MPType}
  line::Int
  type::MPType
end
Read(variable, line) = Read(variable, line, MUndefined)

mutable struct Write <: Statement
  arguments::Vector{Value}
  line::Int
  type::MPType
end
Write(arguments, line) = Write(arguments, line, MUndefined)

mutable struct Assert <: Statement
  argument::Value
  line::Int
  type::MPType
end
Assert(argument, line) = Assert(argument, line, MUndefined)

mutable struct Return <: Statement
  value::Value
  line::Int
  type::MPType
end
Return(value, line) = Return(value, line, MUndefined)

mutable struct SizeFactor <: Factor
  array::Factor
  line::Int
  type::MPType
  end
SizeFactor(array, line) = SizeFactor(array, line, MUndefined)

mutable struct NotFactor <: Factor
  argument::Factor
  line::Int
  type::MPType
end
NotFactor(argument, line) = NotFactor(argument, line, MUndefined)

mutable struct ParenFactor <: Factor
  expression::Value
  line::Int
  type::MPType
end
ParenFactor(expression, line) = ParenFactor(expression, line, MUndefined)

mutable struct LiteralFactor <: Factor
  token::Token
  unique_id::String  # Used to declare string literals
  line::Int
  type::MPType
end
LiteralFactor(token, index, line) = LiteralFactor(token, index, line, MUndefined)

mutable struct VariableFactor <: Factor
  identifier::String
  variable_entry::Union{SymTableEntry,Nothing}
  line::Int
  type::MPType
end
VariableFactor(identifier::String, line::Int) = VariableFactor(identifier, nothing, line, MUndefined)

mutable struct GetRef <: Value
  operand::VariableFactor
  line::Int
  type::MPType
end
GetRef(operand::VariableFactor, line::Int) = GetRef(operand, line, MUndefined)

mutable struct ArrayAccessFactor <: Factor
  identifier::String
  index::Value
  variable_entry::Union{SymTableEntry,Nothing}
  line::Int
  type::MPType
end
ArrayAccessFactor(identifier, index, line) = ArrayAccessFactor(identifier, index, nothing, line, MUndefined)

mutable struct CallFactor <: Factor
  identifier::String
  arguments::Vector{Value}
  implicit_params::Vector{SymTableEntry}
  subroutine::Union{Subroutine,Nothing}
  call_id::Int
  fingerprint::String
  scope_level::Int
  line::Int
  type::MPType
end
CallFactor(identifier, arguments, call_id, line) =
  CallFactor(identifier, arguments, Vector{String}(), nothing, call_id, "", 0, line, MUndefined)

Call = Union{CallFactor,CallStatement}

relational_operators = [
  equals,
  not_equal,
  less_than,
  less_than_or_equal,
  greater_than,
  greater_than_or_equal
]

adding_operators = [
  plus,
  minus,
  kw_or
]

mult_operators = [
  times,
  divide,
  modulo,
  kw_and
]

literals = [
  int_literal,
  string_literal,
  real_literal,
  kw_true,
  kw_false
]

mutable struct ParsingContext
  tokens::Vector{Token}
  index::Int
  signatures::Vector{Subroutine}
  string_literals::Vector{LiteralFactor}
  all_var_names_and_types::Vector{Tuple{String,MPType}}
  id_counter::Int
end
ParsingContext(tokens::Vector{Token}) =
  ParsingContext(tokens, 1, Vector{Subroutine}(), Vector{LiteralFactor}(), Vector{Tuple{String,MPType}}(), 0)

function create_node_id(pc::ParsingContext)
  id = pc.id_counter
  pc.id_counter += 1
  return id
end

function next_token(pc::ParsingContext)
  if pc.index > length(pc.tokens)
    throw(SyntaxException(
      "Failed to parse program. Do all your statements have a terminating semicolon?"))
  end
  return pc.tokens[pc.index]
end

nextclass(pc::ParsingContext) = next_token(pc::ParsingContext).class

# Convenience function. Returns the token, its class and its line number.
function token_class_line(pc::ParsingContext)
  token = next_token(pc::ParsingContext)
  return token, token.class, token.line
end

"""
Consumes the next token, checking that the token class matches.
"""
function match_term(terminal::TokenClass, pc::ParsingContext, current_unit::String)
  DEBUG && println("match_term called in $(current_unit) with terminal $(terminal), next is $(nextclass(pc))")
  token, class, line = token_class_line(pc)
  does_match = terminal == class
  if !does_match
    if terminal == semicolon && class == eoi
      throw(SyntaxException(
        "Unexpected end of input. Did you forget a semicolon?"))
    elseif terminal == identifier && class ∈ values(keywords)
      throw(SyntaxException(
        "Expected a variable identifier but got keyword \"$(token.lexeme)\" on line $(line)."))
    else
      throw(SyntaxException(
        "Failed to parse $(current_unit) on or around line $(line)."))
    end
  end
  pc.index += 1
  return does_match
end

"""
The main function of the parser. Maintains an index ("next") to point
at the token to be processed next. Builds the AST by calling mutually
recursive parsing functions, starting with statements().
"""
parse_input(input::Vector{Token}) = parse_input(input, ParsingContext(input))

function parse_input(input::Vector{Token}, pc::ParsingContext)
  if length(input) < 2
    throw(SyntaxException(
      "The empty string is not a valid program."
    ))
  end
  program(pc)
end

function program(pc::ParsingContext)
  current_unit = "program"
  DEBUG && println("this is $(current_unit), next is ", next_token(pc))
  match_term(kw_program, pc, current_unit)
  match_term(identifier, pc, current_unit)
  match_term(semicolon, pc, current_unit)
  token, class, line = token_class_line(pc)
  defs = Definitions(Vector{Subroutine}(), line)
  if class ∈ [kw_function, kw_procedure]
    defs::Definitions = definitions(pc)
    push!(pc.signatures, defs.defs...)
  end
  main = block(pc)
  match_term(dot, pc, current_unit)
  return Program(defs, main, line)
end

function procedure(pc::ParsingContext)
  current_unit = "procedure"
  token, class, line = token_class_line(pc)
  match_term(kw_procedure, pc, current_unit)
  name = next_token(pc)
  match_term(identifier, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  params = parameters(pc)
  match_term(close_paren, pc, current_unit)
  match_term(semicolon, pc, current_unit)
  body::Block = block(pc)
  body.is_subroutine_block = true
  match_term(semicolon, pc, current_unit)
  return_type = TypeOfVarOrValue(MNothing, MNothing, ImmediateInt(0), line)
  return Subroutine(name.lexeme, params, return_type, body, line)
end

function func(pc::ParsingContext)
  current_unit = "function"
  token, class, line = token_class_line(pc)
  match_term(kw_function, pc, current_unit)
  name = next_token(pc)
  match_term(identifier, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  params = parameters(pc)
  match_term(close_paren, pc, current_unit)
  match_term(colon, pc, current_unit)
  ret_type = var_type(pc)
  match_term(semicolon, pc, current_unit)
  body::Block = block(pc)
  body.is_subroutine_block = true
  match_term(semicolon, pc, current_unit)
  return Subroutine(name.lexeme, params, ret_type, body, line)
end

function parameters(pc::ParsingContext)
  current_unit = "parameters"
  token, class, line = token_class_line(pc)
  params = Vector{Parameter}()
  while nextclass(pc) != close_paren
    push!(params, parameter(pc))
    if nextclass(pc) == comma
      match_term(comma, pc, current_unit)
      if nextclass(pc) == close_paren
        throw(SyntaxException(
          "A parameter list on line $(line) ends in a comma."
        ))
      end
    end
  end
  return params
end

function parameter(pc::ParsingContext)
  current_unit = "parameter"
  token, class, line = token_class_line(pc)
  is_var_par = (class == kw_var)
  is_var_par && match_term(kw_var, pc, current_unit)
  name = next_token(pc)
  match_term(identifier, pc, current_unit)
  match_term(colon, pc, current_unit)
  is_array = false
  if nextclass(pc) == kw_array
    is_array = true
  end
  v_type::TypeOfVarOrValue = var_type(pc)
  return Parameter(name.lexeme, v_type.scalar_type, is_var_par, is_array, v_type.size, line)
end

function definitions(pc::ParsingContext)
  current_unit = "definitions"
  token, class, line = token_class_line(pc)
  defs = Vector{Subroutine}()
  while nextclass(pc) ∈ [kw_procedure, kw_function]
    def = (nextclass(pc) == kw_procedure ? procedure(pc) : func(pc))
    push!(defs, def)
  end
  return Definitions(defs, line)
end

function block(pc::ParsingContext)
  current_unit = "block"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(kw_begin, pc, current_unit)
  stmts = Vector{Statement}()
  while nextclass(pc) != kw_end
    push!(stmts, statement(pc))
    if nextclass(pc) == semicolon
      match_term(semicolon, pc, current_unit)
    end
  end
  match_term(kw_end, pc, current_unit)
  return Block(stmts, line)
end

function statement(pc::ParsingContext)
  current_unit = "statement"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  if class == kw_var
    return var_declaration(pc)
  end
  if class == kw_begin
    return block(pc)
  end
  if class == kw_if
    return if_statement(pc)
  end
  if class == kw_while
    return while_statement(pc)
  end
  return simple_statement(pc)
end

function var_declaration(pc::ParsingContext)
  current_unit = "var_declaration"
  token, class, line = token_class_line(pc)
  match_term(kw_var, pc, current_unit)
  names = Vector{Token}()
  while nextclass(pc) != colon
    push!(names, next_token(pc))
    match_term(identifier, pc, current_unit)
    if nextclass(pc) == comma
      match_term(comma, pc, current_unit)
      if nextclass(pc) == colon
        throw(SyntaxException(
          "Did not expect comma before colon in declaration on line $(line)."
        ))
      end
    end
  end
  match_term(colon, pc, current_unit)
  v_type = var_type(pc)

  # Collect name-type pairs into parsing context
  true_type = v_type.true_type
  for name in names
    entry = (name.lexeme, true_type)
    if entry ∉ pc.all_var_names_and_types
      push!(pc.all_var_names_and_types, entry)
    end
  end

  return Declaration(v_type, names, line)
end

function simple_statement(pc::ParsingContext)
  current_unit = "simple statement"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  if class == kw_return
    return return_statement(pc)
  end
  if class == kw_read
    return read_statement(pc)
  end
  if class == kw_writeln
    return write_statement(pc)
  end
  if class == kw_assert
    return assert_statement(pc)
  end
  if class == identifier
    # The statement is a call or an assignment.
    # Disambiguate by using lookahead.
    if pc.tokens[pc.index+1].class == open_paren
      return call_statement(pc)
    end
    if pc.tokens[pc.index+1].class ∈ [assign, open_sqr_bracket]
      return assignment(pc)
    end
  end
  throw(SyntaxException(
    "Failed to parse statement on line $(line).")
  )
end

function read_statement(pc::ParsingContext)
  current_unit = "read statement"
  token, class, line = token_class_line(pc)
  match_term(kw_read, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  vars = Vector{Value}()
  while nextclass(pc) != close_paren
    push!(vars, variable(pc))
  end
  match_term(close_paren, pc, current_unit)
  return Read(vars, line)
end

function return_statement(pc::ParsingContext)
  current_unit = "return statement"
  token, class, line = token_class_line(pc)
  match_term(kw_return, pc, current_unit)
  if nextclass(pc) ∈ [semicolon, kw_end]
    return Return(Nothing, line)
  end
  return Return(expr(pc), line)
end

function assignment(pc::ParsingContext)
  current_unit = "assignment"
  token, class, line = token_class_line(pc)
  var::Variable = variable(pc)
  match_term(assign, pc, current_unit)
  value = expr(pc)
  return Assignment(var.identifier, var.is_array_access, var.array_index, value, line)
end

# A function or procedure call as a statement.
# The return value is discarded.
function call_statement(pc::ParsingContext)
  current_unit = "call statement"
  token, class, line = token_class_line(pc)
  match_term(identifier, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  args = arguments(pc, token.lexeme)
  match_term(close_paren, pc, current_unit)
  return CallStatement(token.lexeme, args, create_node_id(pc), line)
end

function write_statement(pc::ParsingContext)
  current_unit = "write statement"
  token, class, line = token_class_line(pc)
  match_term(kw_writeln, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  args = arguments(pc, "writeln")
  match_term(close_paren, pc, current_unit)
  return Write(args, line)
end

function if_statement(pc::ParsingContext)
  current_unit = "if_statement"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(kw_if, pc, current_unit)
  condition = expr(pc)
  match_term(kw_then, pc, current_unit)
  then_stmt = statement(pc)
  if nextclass(pc) != kw_else
    return IfThen(condition, then_stmt, line)
  end
  match_term(kw_else, pc, current_unit)
  else_stmt = statement(pc)
  return IfThenElse(condition, then_stmt, else_stmt, line)
end

function while_statement(pc::ParsingContext)
  current_unit = "while_statement"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(kw_while, pc, current_unit)
  condition = expr(pc)
  match_term(kw_do, pc, current_unit)
  do_stmt = statement(pc)
  return While(condition, do_stmt, line)
end

function arguments(pc::ParsingContext, subroutine_name::String)
  current_unit = "arguments"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  subroutine = getfirst(sr -> sr.name == subroutine_name, pc.signatures)
  if subroutine === nothing  # No signature found for subroutine
    if subroutine_name ∈ ["read", "writeln"]
      return read_writeln_args(pc, subroutine_name)
    end
  end
  args = Vector{Value}()
  while nextclass(pc) != close_paren
    push!(args, expr(pc))
    if nextclass(pc) == comma match_term(comma, pc, current_unit) end
  end
  return args
end

function read_writeln_args(pc, subroutine_name::String)
  current_unit = "read_writeln_args"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  args = Vector{Value}()
  while nextclass(pc) != close_paren
    push!(args, expr(pc))
    if nextclass(pc) == comma match_term(comma, pc, current_unit) end
  end
  return args
end

function expr(pc::ParsingContext)
  current_unit = "expression"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  left_simple = simple_expr(pc)
  token, class, line = token_class_line(pc)
  if class ∈ relational_operators
    match_term(class, pc, current_unit)
    right_simple = simple_expr(pc)
    # return RelationalExpression(left_simple, right_simple, token, line)
    return BinaryOperation(class, left_simple, right_simple, line)
  end
  return left_simple
end

function simple_expr(pc::ParsingContext)
  current_unit = "simple expression"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  first_term_sign = Token(plus, "+", line)
  if class ∈ [plus, minus]
    first_term_sign = token
    match_term(class, pc, current_unit)
  end
  first_term = term(pc)
  if first_term_sign.class == minus
    first_term = UnaryOperation(minus, first_term, first_term_sign.line)
  end
  ret_expr = first_term
  while nextclass(pc) ∈ adding_operators
    operation = next_token(pc)
    match_term(operation.class, pc, current_unit)
    new_term = term(pc)
    ret_expr = BinaryOperation(operation.class, ret_expr, new_term, operation.line)
  end
  return ret_expr
end

function term(pc::ParsingContext)
  current_unit = "term"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  factors = Array{Tuple{Factor,Token},1}()
  first_factor = factor(pc)
  ret_term = first_factor
  # push!(factors, (first_factor, Token(identity_op, "identity", line)))
  while nextclass(pc) ∈ mult_operators
    operator = next_token(pc)
    match_term(operator.class, pc, current_unit)
    ret_term = BinaryOperation(operator.class, ret_term, factor(pc), operator.line)
    # push!(factors, (factor(pc), operator))
  end
  # return Term(factors, line)
  return ret_term
end

function factor(pc::ParsingContext)
  current_unit = "factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  fact = Nothing
  if class == open_paren
    fact = paren_factor(pc)
  elseif class == kw_not
    fact = not_factor(pc)
  elseif class ∈ literals
    fact = literal_factor(pc)
  elseif class == identifier
    # The factor is a call or a variable.
    # Disambiguate by using lookahead.
    if pc.tokens[pc.index+1].class == open_paren
      fact = call_factor(pc)
    else
      fact = variable_factor(pc)
    end
  end
  token, class, line = token_class_line(pc)
  if class == dot
    match_term(dot, pc, current_unit)
    match_term(kw_size, pc, current_unit)
    return SizeFactor(fact, line)
  end
  return fact
end

function call_factor(pc::ParsingContext)
  current_unit = "call_factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  id = token
  match_term(identifier, pc, current_unit)
  match_term(open_paren, pc, current_unit)
  args = arguments(pc, id.lexeme)
  match_term(close_paren, pc, current_unit)
  return CallFactor(id.lexeme, args, create_node_id(pc), line)
end

function variable(pc::ParsingContext)
  current_unit = "variable"
  token, class, line = token_class_line(pc)
  id = token
  match_term(identifier, pc, current_unit)
  token, class, _ = token_class_line(pc)
  idx = ImmediateInt(-1)
  is_array_access = false
  if class == open_sqr_bracket
    is_array_access = true
    match_term(open_sqr_bracket, pc, current_unit)
    idx = expr(pc)
    match_term(close_sqr_bracket, pc, current_unit)
  end
  return Variable(id.lexeme, is_array_access, idx)
end

function variable_factor(pc::ParsingContext)
  current_unit = "variable_factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  var::Variable = variable(pc)
  if var.is_array_access
    return ArrayAccessFactor(var.identifier, var.array_index, line)
  end
  return VariableFactor(var.identifier, line)
end

# A boolean negation ("not") factor
function not_factor(pc::ParsingContext)
  current_unit = "not_factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(kw_not, pc, current_unit)
  return NotFactor(factor(pc), line)
end

function literal_factor(pc::ParsingContext)
  current_unit = "literal_factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(class, pc, current_unit)
  unique_id = "@_$(pc.index)_"
  factor = LiteralFactor(token, unique_id, line)
  if class == string_literal push!(pc.string_literals, factor) end
  return factor
end

# Parenthesized expression as factor
function paren_factor(pc::ParsingContext)
  current_unit = "paren_factor"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  match_term(open_paren, pc, current_unit)
  expression = expr(pc)
  match_term(close_paren, pc, current_unit)
  return ParenFactor(expression, line)
end

function var_type(pc::ParsingContext)
  current_unit = "var_type"
  token, class, line = token_class_line(pc)
  DEBUG && println("this is $(current_unit), next is ", token)
  if class == kw_array
    match_term(kw_array, pc, current_unit)
    match_term(open_sqr_bracket, pc, current_unit)
    len = expr(pc)
    match_term(close_sqr_bracket, pc, current_unit)
    match_term(kw_of, pc, current_unit)
    s_type = next_token(pc)
    match_term(s_type.class, pc, current_unit)
    scalar_type = token_class_to_scalar_type[s_type.class]
    return TypeOfVarOrValue(scalar_type, scalar_to_array_type[scalar_type], len, line)
  end
  match_term(class, pc, current_unit)
  return TypeOfVarOrValue(token_class_to_scalar_type[class], token_class_to_scalar_type[class], ImmediateInt(0), line)
end

# function var_as_arg(pc::ParsingContext)
#   current_unit = "var_as_par"
#   token, class, line = token_class_line(pc)
#   DEBUG && println("this is $(current_unit), next is ", token)
#   match_term(identifier, pc, current_unit)
#   return VarAsArgument(token.lexeme, line)
# end

# parse_input(source::String) = parse_input(scan_input(source))

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
@doc doc"""
mlr_beta{Float64}(y::Vector{Float64}, n::Int64=10)

Moving linear regression intercept (column 1) and slope (column 2)
""" ->
function mlr_beta{Float64}(y::Vector{Float64}, n::Int64=10)
    @assert n<length(y) && n>0 "Argument n out of bounds."
    out = zeros(Float64, (length(y),2))
    out[1:n-1,:] = NaN
    xi = collect(1.0:n)
    xbar = mean(xi)
    ybar = runmean(y, n, false)
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        out[i,2] = cov(xi,yi) / var(xi)
        out[i,1] = ybar[i] - out[i,2]*xbar
    end
    return out
end

@doc doc"""
mlr_slope{Float64}(y::Vector{Float64}, n::Int64=10)

Moving linear regression slope
""" ->
function mlr_slope{Float64}(y::Vector{Float64}, n::Int64=10)
    @assert n<length(y) && n>0 "Argument n out of bounds."
    out = zeros(y)
    out[1:n-1] = NaN
    xi = collect(1.0:n)
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        out[i] = cov(xi,yi) / var(xi)
    end
    return out
end

@doc doc"""
mlr_intercept{Float64}(y::Vector{Float64}, n::Int64=10)

Moving linear regression y-intercept
""" ->
function mlr_intercept{Float64}(y::Vector{Float64}, n::Int64=10)
    @assert n<length(y) && n>0 "Argument n out of bounds."
    out = zeros(y)
    out[1:n-1] = NaN
    xi = collect(1.0:n)
    xbar = mean(xi)
    ybar = runmean(y, n, false)
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        out[i] = ybar[i] - xbar*(cov(xi,yi)/var(xi))
    end
    return out
end

@doc doc"""
mlr{Float64}(y::Vector{Float64}, n::Int64=10)

Moving linear regression estimates
""" ->
function mlr{Float64}(y::Vector{Float64}, n::Int64=10)
    b = mlr_beta(y, n)
    return b[:,1] + b[:,2]*float(n)
end

@doc doc"""
mlr_se{Float64}(y::Vector{Float64}, n::Int64=10)

Moving linear regression standard error of estimate
""" ->
function mlr_se{Float64}(y::Vector{Float64}, n::Int64=10)
    yhat = mlr(y, n)
    r = zeros(Float64, n)
    out = zeros(y)
    out[1:n-1] = NaN
    nf = float(n)
    @inbounds for i = n:length(y)
        r = y[i-n+1:i] - yhat[i]
        out[i] = sqrt(sum(r.^2)/nf)
    end
    return out
end

@doc doc"""
mlr_ub{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)

Moving linear regression upper bound
""" ->
function mlr_ub{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)
    return y + se*mlr_se(y,n)
end

@doc doc"""
mlr_lb{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)

Moving linear regression lower bound
""" ->
function mlr_lb{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)
    return y - se*mlr_se(y,n)
end

@doc doc"""
mlr_bands{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)

Moving linear regression bands

`Output:`

Column 1: Lower bound

Column 2: Regression estimate

Column 3: Upper bound
""" ->
function mlr_bands{Float64}(y::Vector{Float64}, n::Int64=10, se::Float64=2.0)
    out = zeros(Float64, (length(y),3))
    out[1:n-1,:] = NaN
    out[:,2] = mlr(y, n)
    out[:,1] = mlr_lb(y, n, se)
    out[:,1] = mlr_ub(y, n, se)
    return out
end

@doc doc"""
mlr_rsq{Float64}(y::Vector{Float64}, n::Int64=10; adjusted::Bool=false)

Moving linear regression R-squared (and adjusted R-squared)
""" ->
function mlr_rsq{Float64}(y::Vector{Float64}, n::Int64=10; adjusted::Bool=false)
    yhat = mlr(y, n)
    rsq = runcor(y, yhat, n, false) .^ 2
    if adjusted
        return rsq - (1.0-rsq)*(1.0/(float(n)-2.0))
    else
        return rsq
    end
end


text:
function set_index!(p::Profile; obs="", var="")
    if obs != ""
        # make obs unique
    end

    if var != ""
        # make var unique
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module LiquidSDR

include("Modems.jl")
include("Filters.jl")

using Reexport
@reexport using  .Modems, .Filters

end


text:
#===============================================================================

    Various random walkers, including:
        ✔ uni- and multivariate Uniform random walk
        ✔ multivariate Gaussian random walk
        ✔ mixture of two Gaussian random walks
    All of the above allow for incorporating a restriction on the positivity
    of any chosen subset of parameter's coordinates.

===============================================================================#

"""
    RandomWalk <: TransitionKernel

Supertype for all random walkers.
"""
abstract type RandomWalk <: TransitionKernel end

#===============================================================================
                            Uniform random walk
===============================================================================#

@doc raw"""
    mutable struct UniformRandomWalk{T,S} <: RandomWalk
        ϵ::T
        pos::S
    end

A uniform random walker that moves all **unrestricted** coordinates according to
the additive rule
```math
x_i+U,
    \quad\mbox{with}
        \quad U\sim\texttt{Unif}([-ϵ_i,ϵ_i]),
            \qquad i\in\texttt{indices},
```
and moves all coordinates **restricted to be positive** according to the rule
```math
x_i\exp(U),
    \quad\mbox{with}
        \quad U\sim\texttt{Unif}([-ϵ_i,ϵ_i]),
            \qquad i\in\texttt{indices}.
```
"""
mutable struct UniformRandomWalk{T,S} <: RandomWalk
    ϵ::T
    pos::S

    function UniformRandomWalk(ϵ, pos=custom_zero(ϵ, Bool))
        @assert all(ϵ .> 0.0)
        new{typeof(ϵ),typeof(pos)}(ϵ, pos)
    end
end

Base.eltype(::UniformRandomWalk{T}) where T = eltype(T)
Base.length(rw::UniformRandomWalk) = length(rw.ϵ)

"""
    Base.rand(rw::UniformRandomWalk, θ)

Sample a new state for a random walker `rw` that is currently in a state `θ`.
"""
Base.rand(rw::UniformRandomWalk, θ) = rand(Random.GLOBAL_RNG, rw, θ)

function Base.rand(
        rng::Random.AbstractRNG,
        rw::UniformRandomWalk,
        θ
    )
    _rand(u) = rand(rng, u)
    U = _rand.(Uniform.(-rw.ϵ, rw.ϵ))
    θ.* ( exp.(U).*rw.pos .+ 1.0.*.!rw.pos ) .+ U.*.!rw.pos
end

function Random.rand!(rng::Random.AbstractRNG, rw::UniformRandomWalk, θ, θ°)
    θ° .= rand(rng, rw, θ)
end

Random.rand!(rw::UniformRandomWalk, θ, θ°) = rand!(Random.GLOBAL_RNG, rw, θ, θ°)

"""
    Distributions.logpdf(rw::UniformRandomWalk, θ, θ°)

Evaluate the log-probability density function for a move of a random walker
from θ to θ°.
"""
function Distributions.logpdf(rw::UniformRandomWalk, θ, θ°)
    mapreduce(
        i->( rw.pos[i] ? -log(2.0*rw.ϵ[i])-log(θ°[i]) : 0.0 ),
        +,
        1:length(rw.ϵ)
    )
end

#===============================================================================
                            Gaussian random walk
===============================================================================#
@doc raw"""
    mutable struct GaussianRandomWalk{T} <: RandomWalk
        Σ::Symmetric{T,Array{T,2}}
        pos::Vector{Bool}
    end

A Gaussian random walker that first transforms all coordinates restricted to be
positive to a log-scale via
```math
θ_i\leftarrow\log(θ_i),\quad i\in\{i: θ_i\texttt{_restricted_to_be_positive}\},
```
then moves according to
```math
θ°∼N(θ, Σ),\quad \forall i,
```
and finally, transforms the restricted coordinates back to the original scale
via
```math
θ°_i\leftarrow\exp(θ°_i),\quad i\in\{i: θ_i\texttt{_restricted_to_be_positive}\}.
```
Restrictions are specified in the vector `pos` (and no restrictions correspond
to all entries in `pos` being `false`).
"""
mutable struct GaussianRandomWalk{T} <: RandomWalk
    Σ::Symmetric{T,Array{T,2}}
    pos::Vector{Bool}

    function GaussianRandomWalk(Σ, pos=nothing)
        @assert size(Σ, 1) == size(Σ, 2)
        Σ = reshape(Σ, (size(Σ, 1), size(Σ, 2)))
        T = eltype(Σ)
        pos = (pos === nothing) ? fill(false, size(Σ, 1)) : pos
        new{T}(Symmetric(Σ), pos)
    end
end

Base.eltype(::GaussianRandomWalk{T}) where T = T
Base.length(rw::GaussianRandomWalk) = size(rw.Σ, 1)
remove_constraints!(rw::GaussianRandomWalk, θ) = (θ[rw.pos] .= log.(θ[rw.pos]))
reimpose_constraints!(rw::GaussianRandomWalk, θ) = (
    θ[rw.pos] .= exp.(θ[rw.pos])
)

Base.rand(rw::GaussianRandomWalk, θ) = rand(Random.GLOBAL_RNG, rw, θ)

function Base.rand(rng::Random.AbstractRNG, rw::GaussianRandomWalk, θ)
    remove_constraints!(rw, θ)
    θ° = rand(MvNormal(θ, rw.Σ))
    reimpose_constraints!(rw, θ°)
    reimpose_constraints!(rw, θ)
    θ°
end

function Random.rand!(rw::GaussianRandomWalk, θ, θ°)
    rand!(Random.GLOBAL_RNG, rw, θ, θ°)
end

function Random.rand!(rng::Random.AbstractRNG, rw::GaussianRandomWalk, θ, θ°)
    θ° .= rand(rng, rw, θ)
end

_logjacobian(rw::GaussianRandomWalk, θ) = -sum(log.(θ[rw.pos]))

function Distributions.logpdf(rw::GaussianRandomWalk, θ, θ°)
    logJ = _logjacobian(rw, θ°)
    remove_constraints!(rw, θ)
    remove_constraints!(rw, θ°)
    lpdf = logpdf(MvNormal(θ, rw.Σ), θ°) + logJ
    reimpose_constraints!(rw, θ)
    reimpose_constraints!(rw, θ°)
    lpdf
end

#===============================================================================
            A mixture of two multidimensional Gaussian random walks
===============================================================================#
@doc raw"""
    mutable struct GaussianRandomWalkMix{T} <: RandomWalk
        gsn_A::GaussianRandomWalk{T}
        gsn_B::GaussianRandomWalk{T}
        λ::Float64
    end

A mixture of two Gaussian random walkers. It performs updates according to
```math
X_{i+1} = X_i + λZ_A + (1-λ)Z_B,
```
where
```math
Z_A ∼ N(0,Σ_A),\quad\mbox{and}\quad Z_B ∼ N(0, Σ_B),
```
(with $X_i$ being already appropriately log-transformed if needed).
"""
mutable struct GaussianRandomWalkMix{T} <: RandomWalk
    gsn_A::GaussianRandomWalk{T}
    gsn_B::GaussianRandomWalk{T}
    λ::Float64

    function GaussianRandomWalkMix(Σ_A, Σ_B, λ=0.5, pos=nothing)
        @assert 0.0 <= λ <= 1.0
        @assert size(Σ_A) == size(Σ_B)
        gsn_A = GaussianRandomWalk(Σ_A, pos)
        gsn_B = GaussianRandomWalk(Σ_B, pos)
        T = eltype(gsn_A)
        new{T}(gsn_A, gsn_B, λ)
    end
end

Base.eltype(::GaussianRandomWalkMix{T}) where T = T
Base.length(rw::GaussianRandomWalkMix) = length(rw.gsn_A)

Base.rand(rw::GaussianRandomWalkMix, θ) = rand(Random.GLOBAL_RNG, rw, θ)

function Base.rand(rng::Random.AbstractRNG, rw::GaussianRandomWalkMix, θ)
    rw_i = pick_kernel(rng, rw)
    rand(rng, rw_i, θ)
end

Random.rand!(rw::GaussianRandomWalkMix, θ, θ°) = rand!(Random.GLOBAL_RNG, rw, θ, θ°)

function Random.rand!(rng::Random.AbstractRNG, rw::GaussianRandomWalkMix, θ, θ°)
    rw_i = pick_kernel(rng, rw)
    rand!(rng, rw_i, θ, θ°)
end

function pick_kernel(rng::Random.AbstractRNG, rw::GaussianRandomWalkMix)
    rand(rng, Bernoulli(rw.λ)) ? rw.gsn_B : rw.gsn_A
end

function Distributions.logpdf(rw::GaussianRandomWalkMix, θ, θᵒ)
    log( (1-rw.λ)*exp(logpdf(rw.gsn_A, θ, θᵒ))
          + rw.λ *exp(logpdf(rw.gsn_B, θ, θᵒ)) )
end


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @save

has_TM_flag = 0

# load synthetic data -
path_to_input = joinpath(_PATH_TO_DATA, "synthetic_input_array.dat")
synthetic_input_data = CSV.read(path_to_input, DataFrame)

path_to_output = joinpath(_PATH_TO_DATA, "synthetic_output_array.dat")
synthetic_output_data = CSV.read(path_to_output, DataFrame)

# initialize storage for the training data -
training_data = Vector{Tuple{Vector{Float32},Vector{Float32}}}()

# build a model architecture -
deep_coag_model = Chain(Dense(12, 16, σ), Dense(16, 5));

# setup a loss function -
loss(x, y) = Flux.Losses.mae(deep_coag_model(x), y; agg = mean)

# pointer to params -
ps = Flux.params(deep_coag_model)

# # use old school gradient descent -
opt = Momentum(0.1, 0.95)

# main training loop -
(P, D) = size(synthetic_input_data)
for i ∈ 1:P

    for j = 1:P

        # leave one out -
        if (i != j)

            # get input and output data -
            input_data = convert.(Float32, Vector(synthetic_input_data[j, :]))
            output_data = convert.(Float32, Vector(synthetic_output_data[j, :]))
            data_example = (input_data, output_data)

            # capture -
            push!(training_data, data_example)
        end
    end

    # ok, so have the training data for this case -
    # train -
    @epochs 12000 Flux.train!(loss, ps, training_data, opt)

    # save -
    model_name = "deep_coag_model-L$(i)O-TF-TM-$(has_TM_flag)-synthetic-V2-V3.bson"
    model_file_path = joinpath(_PATH_TO_MODELS, model_name)
    @save model_file_path deep_coag_model

    # need to empty the training data -
    empty!(training_data)
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module Fill_To_Max_Time

export fill_to_max_time

using AdaptiveStressTesting: ASTAction
using RLESUtils, Obj2Dict, RNGWrapper, Loggers, Obj2DataFrames

using ..DefineSave
using ..SaveHelpers

function fill_to_max_time(filename::AbstractString)
    d = trajLoad(filename)

    # disable ending on nmac
    sim_params = get_sim_params(d)
    sim_params.end_on_nmac = false
    set_sim_params!(d, sim_params)

    action_seq = get_action_seq(d)
    ast_params = get_ast_params(d)
    rsg_length = ast_params.rsg_length

    max_steps = sim_params.max_steps
    steps_to_append = max_steps - length(action_seq)  #determine number of missing steps

    # steps_to_append > 0 check is automatically handled by comprehension
    seed = steps_to_append #arbitrarily use this as a seed
    actions_to_append = ASTAction[ ASTAction(rsg_length, seed) for t = 1:steps_to_append ] #append hash of t
    action_seq = vcat(action_seq, actions_to_append)

    set_action_seq!(d, action_seq)
    outfilename = trajSave(string(getSaveFileRoot(filename), "_filled"), d)
    println("File: ", filename, "; Steps appended: ", steps_to_append)

    outfilename
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
# Contained functions in file: bldup!, difun!

function bldup!(r,dy)
#=**************************************************************************
     This is a function in an ozone NOx formaldehyde mechanism, adapted from Fortran
     by Obin Sturm (UC Davis) in October 2019.
     Written in Julia Version 1.2.0, tested up to 1.4.2

     This module was recieved from Mike Kleeman (UC Davis) in October 2018.
     The original comments and file description are below.

     MODEL DEPENDENT SUBROUTINE TO CALCULATE FORMATION
     RATES OF NON-REACTING SPECIES

     Inputs:
     r   - reaction rates for the mechanism (get from difun)
         1) NO2 + HV = NO + O
         2) O + O2 = O3
         3) O3 + NO = NO2 + O2
         4) HCHO + HV = 2 HO2. + CO
         5) HCHO + HV = H2 + CO
         6) HCHO + HO. = HO2. + CO + H2O
         7) HO2. + NO = HO. + NO2
         8) HO. + NO2 = HNO3
         9) HO2H + HV = 2 HO.
        10) HO2H + HO. = H2O + HO2.

      Outputs:
         dy  - formation rates of the buildup species [ppm min-1]
**************************************************************************=#
    # HNO3
    dy[1] = r[8]

    # CO
    dy[2] = r[4]+r[5]+r[6]

    # HO
    dy[3] = r[5]

    return dy
end

function  difun!(c,a,s,rk,r,fr,lr)
    #=**************************************************************************
    This is a function in an ozone NOx formaldehyde mechanism, adapted from Fortran
    by Obin Sturm (UC Davis) in October 2019.

    This module was recieved from Mike Kleeman (UC Davis) in October 2018.
    The original comments and file description are below.
    MODEL DEPENDENT SUBROUTINE TO CALCULATE RATES OF
    FORMATION AND LOSS FOR ACTIVE SPECIES

        dA/dt = fr - lr*a

    Inputs
     c   - constant species
            c(1) = O2 (set to atmospheric value = 2.090E+05ppm)
            c(2) = HV (leave set to 1)
     a   - active species concentrations
            a(1) = O3
            a(2) = NO
            a(3) = NO2
            a(4) = HCHO
            a(5) = HO2.
            a(6) = HO2H
     rk  - reaction constants for the mechanism
            1) NO2 + HV = NO + O
            2) O + O2 = O3
            3) O3 + NO = NO2 + O2
            4) HCHO + HV = 2 HO2. + CO
            5) HCHO + HV = H2 + CO
            6) HCHO + HO. = HO2. + CO + H2O
            7) HO2. + NO = HO. + NO2
            8) HO. + NO2 = HNO3
            9) HO2H + HV = 2 HO.
           10) HO2H + HO. = H2O + HO2.

    Outputs:
     s   - steady state species concentrations
            s(1) = HO.
            s(2) = O

     r   - reaction rates for mechanism listed above
     fr  - formation rates for each of the active species (ppm min-1)
     lr  - loss rates for each of the active species (ppm min-1)
    **************************************************************************=#

    # DEFINE RATES OF REACTIONS, PARTIAL RATES
    # FOR REACTIONS WITH STEAdy STATE REACTANTS

     r[1] = rk[1]*a[3]*c[2]
     r[2] = rk[2]*c[1]
     r[3] = rk[3]*a[1]*a[2]
     r[4] = rk[4]*a[4]*c[2]
     r[5] = rk[5]*a[4]*c[2]
     r[6] = rk[6]*a[4]
     r[7] = rk[7]*a[5]*a[2]
     r[8] = rk[8]*a[3]
     r[9] = rk[9]*a[6]*c[2]
     r[10] = rk[10]*a[6]

     # CALCULATE STEADY STATE CONCENTRATIONS AND
     # RATES OF REACTION WITH SS SPECIES

     # HO.
     s[1] = (r[7]+r[9]*2.0)/(1.0e-20+r[6]+r[8]+r[10])
     r[6] = s[1]*r[6]
     r[8] = s[1]*r[8]
     r[10] = s[1]*r[10]

     # O
     s[2] = r[1]/(1e-20+r[2])
     r[2] = s[2]*r[2]

     # DEFINE FORMATION RATES FOR ACTIVE SPECIES

     # O3
     fr[1] = r[2]

     # NO
     fr[2] = r[1]

     # NO2
     fr[3] = r[3] + r[7]

     # HCHO
     fr[4] = 0.0

     # HO2.
     fr[5] = 2*r[4] + r[6] + r[10]

     # HO2H
     fr[6] = 0.0

     # DEFINE LOSS RATES OF ACTIVE species

     # O3
     lr[1] = rk[3]*a[2]

     # NO
     lr[2] = rk[3]*a[1] + rk[7]*a[5]

     # NO2
     lr[3] = rk[1]*c[2] + rk[8]*s[1]

     # HCHO
     lr[4] =  rk[4]*c[2]+rk[5]*c[2]+rk[6]*s[1]

     # HO2.
     lr[5] = rk[7]*a[2]

     # HO2H
     lr[6] =  rk[9]*c[2]+rk[10]*s[1]
     return s,r,fr,lr
 end