text:
using MiniFB
# A simple function that renders some noise to a window. Demonstrates how to create a window and render a buffer into it using MiniFB
function noise()
    WIDTH = 800
    HEIGHT = 600
    g_buffer = zeros(UInt32, WIDTH * HEIGHT)
    noise = carry = seed = 0xbeef
    window = mfb_open_ex("Noise Test", WIDTH, HEIGHT, MiniFB.WF_RESIZABLE);
    while mfb_wait_sync(window)
        for i in 1:WIDTH * HEIGHT
            noise = seed;
            noise = noise >> 3;
            noise = noise ^ seed;
            carry = noise & 1;
            noise = noise >> 1;
            seed = seed >> 1;
            seed = seed | (carry << 30);
            noise = noise & 0xFF;
            g_buffer[i] = mfb_rgb(noise, noise, noise);
        end
        state = mfb_update(window, g_buffer);
        if state != MiniFB.STATE_OK
            break;
        end
    end
    mfb_close(window)
end

# Call the function
noise()

# ![](../assets/noise.png)


text:
include("../src/GraphIdx.jl")

using Test
import .GraphIdx.Tree: lowest_common_ancestors, ChildrenIndex, IncidenceIndex

@testset "LCA IncidenceIndex             " begin
    pairs = [(5, 7), (5, 2), (5, 3), (2, 3)]
    idx = IncidenceIndex(7, pairs)
    @test Set(idx[5]) == Set([(7, 1), (2, 2), (3, 3)])
    @test Set(idx[2]) == Set([(5, 2), (3, 4)])
    idx2 = IncidenceIndex(7, [(5, 2)])
    @test Set(idx2[2]) == Set([(5, 1)])
end

@testset "LCA Wikipedia (problematic)    " begin
    pi = [1, 1, 1, 2, 4, 4, 6]
    pairs = [(5, 2)]
    @test lowest_common_ancestors(ChildrenIndex(pi), pi, pairs) == [2]
end

@testset "LCA Wikipedia                  " begin
    #     1  2  3  4  5  6  7
    pi = [1, 1, 1, 2, 4, 4, 6]
    #=
         1
        / \
       /   \
      2     3
      |
      4
     / \
    5   6
        |
        7
    =#
    idx = ChildrenIndex(pi)
    @test Set(idx[1]) == Set([2, 3])
    @test Set(idx[2]) == Set([4])
    @test Set(idx[3]) == Set([])
    @test Set(idx[4]) == Set([5, 6])
    @test Set(idx[5]) == Set([])
    @test Set(idx[6]) == Set([7])
    @test Set(idx[7]) == Set([])
    @test length(idx) == 7
    pairs = [(5, 7), (5, 2), (5, 3), (2, 3)]
    ans = [4, 2, 1, 1]
    @test length(pairs) == length(ans)
    @test lowest_common_ancestors(idx, pi, pairs) == ans
end

ext: jl

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
function expmk(k::Float64, op::AbstractVector{Float64})
    diagm(exp.(k .* op))
end

function expmk(k::Float64, op::AbstractMatrix{Float64})
    exp(k .* op)
end

"""
    symmetric_trotter_path(num_links::Int, tau::Float64, ops::AbstractArray...)

Compute an approximation of `exp(-num_links*tau*sum(ops))` using a symmetric
Trotter decomposition.

The first term in `ops` is placed in the middle of the decomposition.
"""
function symmetric_trotter_path(num_links::Int, tau::Float64, ops::AbstractArray...)
    step = I(size(ops[1])[1])

    for i in length(ops):-1:2
        step *= expmk(-tau/2, ops[i])
    end
    step *= expmk(-tau, ops[1])
    for i in 2:length(ops)
        step *= expmk(-tau/2, ops[i])
    end

    step^num_links
end

"""
    symmetric_trotter_path(num_links::Int, tau::Float64, ops::Dict...)

Compute an approximation of `exp(-num_links*tau*sum(ops))` using a symmetric
Trotter decomposition.

The first term in `ops` is placed in the middle of the decomposition.
"""
function symmetric_trotter_path(num_links::Int, tau::Float64, ops::Dict...)
    path = Dict{BlockLabel,Matrix{Float64}}()

    for label in keys(ops[1])
        step = I(size(ops[1][label])[1])

        for i in length(ops):-1:2
            step *= expmk(-tau/2, ops[i][label])
        end
        step *= expmk(-tau, ops[1][label])
        for i in 2:length(ops)
            step *= expmk(-tau/2, ops[i][label])
        end

        path[label] = step^num_links
    end

    path
end

"""
    wf_pigs(half_path::Matrix{Float64}, trial_idx::Int)

Extract the propagated ground state wavefunction from the middle of a PIGS
path.

To access the middle of the path, we require only one half of the path in
`half_path`. The half path is multiplied by a trial function, which is the
basis vector specified by `trial_idx`.
"""
function wf_pigs(half_path::Matrix{Float64}, trial_idx::Int)
    wf = half_path[:, trial_idx]
    wf ./= sqrt(sum(wf.^2))

    wf
end

mult_by_op(op1::AbstractMatrix{Float64}, op2::AbstractVector{Float64}) = op1 * diagm(op2)
mult_by_op(op1::AbstractMatrix{Float64}, op2::AbstractMatrix{Float64}) = op1 * op2

"""
    energy(path::Dict, ops::Dict...)

Compute the thermal energy of a finite temperature `path` for a Hamiltonian
with terms `ops`.
"""
function energy(path::Dict, ops::Dict...)
    Z = 0.0
    avgs = zeros(Float64, length(ops))

    for label in keys(path)
        Z += tr(path[label])

        for i in 1:length(ops)
            avgs[i] += tr(mult_by_op(path[label], ops[i][label]))
        end
    end

    avgs /= Z
    E = sum(avgs)

    Z, E, avgs
end

"""
    energy_mixed(path::Matrix{Float64}, trial_idx::Int, ops::AbstractArray...)

Compute the ground state energy of a PIGS `path` using the mixed estimator for
a Hamiltonian with terms `ops`.

The path is multiplied on both ends by a trial function, which is the basis
vector specified by `trial_idx`.
"""
function energy_mixed(path::Matrix{Float64}, trial_idx::Int, ops::AbstractArray...)
    avgs = Array{Float64}(undef, length(ops))

    for i in 1:length(ops)
        avgs[i] = mult_by_op(path, ops[i])[trial_idx, trial_idx]
    end

    avgs /= path[trial_idx, trial_idx]
    E0 = sum(avgs)

    E0, avgs
end


text:
using JuMP
using GLPKMathProgInterface
using Ipopt

n = 4;
lowRes = -2;
hiRes  = 2;
m = Model(solver=IpoptSolver())
#@setNLObjective(m, Max,0)

@defVar(m, lowRes <= α[i=1:n]  <=hiRes)
@defVar(m, lowRes <= β1[i=1:n] <=hiRes)
@defVar(m, lowRes <= β2[i=1:n] <=hiRes)
@defVar(m, lowRes <= β3[i=1:n] <=hiRes)
@defVar(m, lowRes <= β4[i=1:n] <=hiRes)

@defVar(m, lowRes <= A011 <= hiRes)
@defVar(m, lowRes <= A012 <= hiRes)
@defVar(m, lowRes <= A013 <= hiRes)
@defVar(m, lowRes <= A014 <= hiRes)
@defVar(m, lowRes <= A021 <= hiRes)
@defVar(m, lowRes <= A022 <= hiRes)
@defVar(m, lowRes <= A023 <= hiRes)
@defVar(m, lowRes <= A024 <= hiRes)
@defVar(m, lowRes <= A031 <= hiRes)
@defVar(m, lowRes <= A032 <= hiRes)
@defVar(m, lowRes <= A033 <= hiRes)
@defVar(m, lowRes <= A034 <= hiRes)
@defVar(m, lowRes <= A041 <= hiRes)
@defVar(m, lowRes <= A042 <= hiRes)
@defVar(m, lowRes <= A043 <= hiRes)
@defVar(m, lowRes <= A044 <= hiRes)

@defVar(m, lowRes <= A111 <= hiRes)
@defVar(m, lowRes <= A112 <= hiRes)
@defVar(m, lowRes <= A113 <= hiRes)
@defVar(m, lowRes <= A114 <= hiRes)
@defVar(m, lowRes <= A121 <= hiRes)
@defVar(m, lowRes <= A122 <= hiRes)
@defVar(m, lowRes <= A123 <= hiRes)
@defVar(m, lowRes <= A124 <= hiRes)
@defVar(m, lowRes <= A131 <= hiRes)
@defVar(m, lowRes <= A132 <= hiRes)
@defVar(m, lowRes <= A133 <= hiRes)
@defVar(m, lowRes <= A134 <= hiRes)
@defVar(m, lowRes <= A141 <= hiRes)
@defVar(m, lowRes <= A142 <= hiRes)
@defVar(m, lowRes <= A143 <= hiRes)
@defVar(m, lowRes <= A144 <= hiRes)

@defVar(m, lowRes <= B011 <= hiRes)
@defVar(m, lowRes <= B012 <= hiRes)
@defVar(m, lowRes <= B013 <= hiRes)
@defVar(m, lowRes <= B014 <= hiRes)
@defVar(m, lowRes <= B021 <= hiRes)
@defVar(m, lowRes <= B022 <= hiRes)
@defVar(m, lowRes <= B023 <= hiRes)
@defVar(m, lowRes <= B024 <= hiRes)
@defVar(m, lowRes <= B031 <= hiRes)
@defVar(m, lowRes <= B032 <= hiRes)
@defVar(m, lowRes <= B033 <= hiRes)
@defVar(m, lowRes <= B034 <= hiRes)
@defVar(m, lowRes <= B041 <= hiRes)
@defVar(m, lowRes <= B042 <= hiRes)
@defVar(m, lowRes <= B043 <= hiRes)
@defVar(m, lowRes <= B044 <= hiRes)

@defVar(m, lowRes <= B111 <= hiRes)
@defVar(m, lowRes <= B112 <= hiRes)
@defVar(m, lowRes <= B113 <= hiRes)
@defVar(m, lowRes <= B114 <= hiRes)
@defVar(m, lowRes <= B121 <= hiRes)
@defVar(m, lowRes <= B122 <= hiRes)
@defVar(m, lowRes <= B123 <= hiRes)
@defVar(m, lowRes <= B124 <= hiRes)
@defVar(m, lowRes <= B131 <= hiRes)
@defVar(m, lowRes <= B132 <= hiRes)
@defVar(m, lowRes <= B133 <= hiRes)
@defVar(m, lowRes <= B134 <= hiRes)
@defVar(m, lowRes <= B141 <= hiRes)
@defVar(m, lowRes <= B142 <= hiRes)
@defVar(m, lowRes <= B143 <= hiRes)
@defVar(m, lowRes <= B144 <= hiRes)

#=
A0 = [A011 A012
     A021 A022]
A1 = [A111 A112
     A121 A122]
B0 = [B011 B012
      B021 B022]
B1 = [B111 B112
      B121 B122]
e = [1 1]'
=#

#@setObjective(m, Max, 0)
@addConstraints(m,begin
                  sum{α[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β2[i], i=1:n}==0
                  sum{β3[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                end)
@addNLConstraints(m,begin
                  B111*β1[1]+B112*β1[1]+B113*β1[1]+B114*β1[1] + B121*β1[2]+B122*β1[2]+B123*β1[2]+B124*β1[2] + B131*β1[3]+B132*β1[3]+B133*β1[3]+B134*β1[3]+ B141*β1[4]+B142*β1[4]+B143*β1[4]+B144*β1[4] ==0
                  B111*β2[1]+B112*β2[1]+B113*β2[1]+B114*β2[1] + B121*β2[2]+B122*β2[2]+B123*β2[2]+B124*β2[2] + B131*β2[3]+B132*β2[3]+B133*β2[3]+B134*β2[3]+ B141*β2[4]+B142*β2[4]+B143*β2[4]+B144*β2[4] ==1
                  B111*β3[1]+B112*β3[1]+B113*β3[1]+B114*β3[1] + B121*β3[2]+B122*β3[2]+B123*β3[2]+B124*β3[2] + B131*β3[3]+B132*β3[3]+B133*β3[3]+B134*β3[3]+ B141*β3[4]+B142*β3[4]+B143*β3[4]+B144*β3[4] ==0
                  B111*β4[1]+B112*β4[1]+B113*β4[1]+B114*β4[1] + B121*β4[2]+B122*β4[2]+B123*β4[2]+B124*β4[2] + B131*β4[3]+B132*β4[3]+B133*β4[3]+B134*β4[3]+ B141*β4[4]+B142*β4[4]+B143*β4[4]+B144*β4[4] ==0
                  A011*α[1]+A012*α[1]+A013*α[1]+A014*α[1]     + A021*α[2]+A022*α[2]+A023*α[2]+A024*α[2]    + A031*α[3]+A032*α[3]+A033*α[3]+A034*α[3] + A041*α[4]+A042*α[4]+A043*α[4]+A044*α[4] ==1//2
                  B011*α[1]+B012*α[1]+B013*α[1]+B014*α[1]     + B021*α[2]+B022*α[2]+B023*α[2]+B024*α[2]    + B031*α[3]+B032*α[3]+B033*α[3]+B034*α[3] + B041*α[4]+B042*α[4]+B043*α[4]+B044*α[4] ==1
                  A111*β1[1]+A112*β1[1]+A113*β1[1]+A114*β1[1] + A121*β1[2]+A122*β1[2]+A123*β1[2]+A124*β1[2] + A131*β1[3]+A132*β1[3]+A133*β1[3]+A134*β1[3] + A141*β1[4]+A142*β1[4]+A143*β1[4]+A144*β1[4]==1
                  A111*β2[1]+A112*β2[1]+A113*β2[1]+A114*β2[1] + A121*β2[2]+A122*β2[2]+A123*β2[2]+A124*β2[2] + A131*β2[3]+A132*β2[3]+A133*β2[3]+A134*β2[3] + A141*β2[4]+A142*β2[4]+A143*β2[4]+A144*β2[4]==0
                  A111*β3[1]+A112*β3[1]+A113*β3[1]+A114*β3[1] + A121*β3[2]+A122*β3[2]+A123*β3[2]+A124*β3[2] + A131*β3[3]+A132*β3[3]+A133*β3[3]+A134*β3[3] + A141*β3[4]+A142*β3[4]+A143*β3[4]+A144*β3[4]==-1
                  A111*β4[1]+A112*β4[1]+A113*β4[1]+A114*β4[1] + A121*β4[2]+A122*β4[2]+A123*β4[2]+A124*β4[2] + A131*β4[3]+A132*β4[3]+A133*β4[3]+A134*β4[3] + A141*β4[4]+A142*β4[4]+A143*β4[4]+A144*β4[4]==0
                end)

@addNLConstraint(m,α[1]  * (B011+B012+B013+B014)*(B011+B012+B013+B014)  + α[2]  * (B021+B022+B123+B124)*(B021+B022+B123+B124) + α[4]  * (B041+B042+B143+B144)*(B041+B042+B143+B144)  == 3//2)
@addNLConstraint(m,β1[1] * (B111+B112+B013+B014)*(B111+B112+B013+B014)  + β1[2] * (B121+B122+B123+B124)*(B121+B122+B123+B124) + β1[4] * (B141+B142+B143+B144)*(B141+B142+B143+B144)  == 1)
@addNLConstraint(m,β2[1] * (B111+B112+B013+B014)*(B111+B112+B013+B014)  + β2[2] * (B121+B122+B123+B124)*(B121+B122+B123+B124) + β2[4] * (B141+B142+B143+B144)*(B141+B142+B143+B144)  == 0)
@addNLConstraint(m,β3[1] * (B111+B112+B013+B014)*(B111+B112+B013+B014)  + β3[2] * (B121+B122+B123+B124)*(B121+B122+B123+B124) + β3[4] * (B141+B142+B143+B144)*(B141+B142+B143+B144)  == -1)
@addNLConstraint(m,β4[1] * (B111+B112+B013+B014)*(B111+B112+B013+B014)  + β4[2] * (B121+B122+B123+B124)*(B121+B122+B123+B124) + β4[4] * (B141+B142+B143+B144)*(B141+B142+B143+B144)  == 2)

@addNLConstraint(m,β1[1] * (B111*(B111+B112+B113+B114) + B112*(B121+B122+B123+B124) + B113*(B131+B132+B133+B134) + B114*(B141+B142+B143+B144))  + β1[2] * (B121*(B111+B112+B113+B114) + B122*(B121+B122+B123+B124) + B123*(B131+B132+B133+B134) + B124*(B141+B142+B143+B144) )  + β1[3] * (B131*(B111+B112+B113+B114) + B132*(B121+B122+B123+B124) + B133*(B131+B132+B133+B134) + B134*(B141+B142+B143+B144)) + β1[4] * (B141*(B111+B112+B113+B114) + B142*(B121+B122+B123+B124) + B143*(B131+B132+B133+B134) + B144*(B141+B142+B143+B144)) == 0)
@addNLConstraint(m,β2[1] * (B111*(B111+B112+B113+B114) + B112*(B121+B122+B123+B124) + B113*(B131+B132+B133+B134) + B114*(B141+B142+B143+B144))  + β2[2] * (B121*(B111+B112+B113+B114) + B122*(B121+B122+B123+B124) + B123*(B131+B132+B133+B134) + B124*(B141+B142+B143+B144) )  + β2[3] * (B131*(B111+B112+B113+B114) + B132*(B121+B122+B123+B124) + B133*(B131+B132+B133+B134) + B134*(B141+B142+B143+B144)) + β2[4] * (B141*(B111+B112+B113+B114) + B142*(B121+B122+B123+B124) + B143*(B131+B132+B133+B134) + B144*(B141+B142+B143+B144)) == 0)
@addNLConstraint(m,β3[1] * (B111*(B111+B112+B113+B114) + B112*(B121+B122+B123+B124) + B113*(B131+B132+B133+B134) + B114*(B141+B142+B143+B144))  + β3[2] * (B121*(B111+B112+B113+B114) + B122*(B121+B122+B123+B124) + B123*(B131+B132+B133+B134) + B124*(B141+B142+B143+B144) )  + β3[3] * (B131*(B111+B112+B113+B114) + B132*(B121+B122+B123+B124) + B133*(B131+B132+B133+B134) + B134*(B141+B142+B143+B144)) + β3[4] * (B141*(B111+B112+B113+B114) + B142*(B121+B122+B123+B124) + B143*(B131+B132+B133+B134) + B144*(B141+B142+B143+B144)) == 0)
@addNLConstraint(m,β4[1] * (B111*(B111+B112+B113+B114) + B112*(B121+B122+B123+B124) + B113*(B131+B132+B133+B134) + B114*(B141+B142+B143+B144))  + β4[2] * (B121*(B111+B112+B113+B114) + B122*(B121+B122+B123+B124) + B123*(B131+B132+B133+B134) + B124*(B141+B142+B143+B144) )  + β4[3] * (B131*(B111+B112+B113+B114) + B132*(B121+B122+B123+B124) + B133*(B131+B132+B133+B134) + B134*(B141+B142+B143+B144)) + β4[4] * (B141*(B111+B112+B113+B114) + B142*(B121+B122+B123+B124) + B143*(B131+B132+B133+B134) + B144*(B141+B142+B143+B144)) == 1)

@addNLConstraint(m,(1/2)*((A111*(B011+B012+B013+B014) + A112*(B021+B022+B023+B024) + A113*(B031+B032+B033+B034) + A114*(B041+B042+B043+B044))*β1[1] +
                          (A121*(B011+B012+B013+B014) + A122*(B021+B022+B023+B024) + A123*(B031+B032+B033+B034) + A124*(B041+B042+B043+B044))*β1[2] +
                          (A131*(B011+B012+B013+B014) + A132*(B021+B022+B023+B024) + A133*(B031+B032+B033+B034) + A134*(B041+B042+B043+B044))*β1[3] +
                          (A141*(B011+B012+B013+B014) + A142*(B021+B022+B023+B024) + A143*(B031+B032+B033+B034) + A144*(B041+B042+B043+B044))*β1[4]) +
                   (1/3)*((A111*(B011+B012+B013+B014) + A112*(B021+B022+B023+B024) + A113*(B031+B032+B033+B034) + A114*(B041+B042+B043+B044))*β3[1] +
                          (A121*(B011+B012+B013+B014) + A122*(B021+B022+B023+B024) + A123*(B031+B032+B033+B034) + A124*(B041+B042+B043+B044))*β3[2] +
                          (A131*(B011+B012+B013+B014) + A132*(B021+B022+B023+B024) + A133*(B031+B032+B033+B034) + A134*(B041+B042+B043+B044))*β3[3] +
                          (A141*(B011+B012+B013+B014) + A142*(B021+B022+B023+B024) + A143*(B031+B032+B033+B034) + A144*(B041+B042+B043+B044))*β3[4]) == 0)

print(m)

status = solve(m)

println("Objective value: ", getObjectiveValue(m))
println("α1 = ", getValue(α))
println("β1 = ", getValue(β1))
println("β2 = ", getValue(β2))
println("β3 = ", getValue(β3))
println("β4 = ", getValue(β4))


text:
function pixelate(s::Simulation, dim::Int64, sparsity_level::T, cushion::T) where T<:Real
	nsteps = s.steptotal[1]
	# either add percent occupation, make new function for it
	im = zeros(dim,dim)
	# ascertain limits of distribution - 4 sigma away!
	localmin = -4*sparsity_level
	localmax = 4*sparsity_level
	# create new bounds for scaling purposes
	globmin = localmin*(1+cushion)
	globmax = localmax*(1+cushion)
	# normalize and scale cell origins to pixel locations
	livingcellIDs = []
	for i in 1:length(s.cells)
		if s.cells[i].d < nsteps
			nothing
		else
			push!(livingcellIDs,i)
		end
	end
	o = zeros(T,length(livingcellIDs),2)
	for i in 1:length(livingcellIDs)
		o[i,1] = s.cells[livingcellIDs[i]].p[1] - globmin
		o[i,2] = s.cells[livingcellIDs[i]].p[2] - globmin
	end
	o .*= (dim / (globmax - globmin))
	# normalize radius to image
	r_im = (1/2.33) * (dim / (localmax-localmin))
	# check in terminal that the radius is reasonable
	#println(r_im)

	for x in 1:length(im[1,:])
		for y in 1:length(im[:,1])
			for i in 1:length(o[:,1])
				h = o[i,1]
				k = o[i,2]
				if ((x-h)^2 + (y-k)^2) < (r_im^2)
					im[y,x] = 1
				else
					nothing
				end
			end
		end
	end
	return im
end

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym


text:
#= semantic checking
e.g. mRNA_F activates g_H → mRNA can not activate the expression of genes
→ be careful! Biology is error-prone, never be too sure.
=#

#=
F: 1) check biological meaning of each token; 2) have all components for each type,
  get ready for IR generation;
I: array of triplets (verb - generalBioSym Arrays - parameters),
   conversion table, sentence specification (encoded in code);
O: can get rid of type and parameters? why shall?
   Array of tuples, tuple = (reaction type, [reactants, products, catalysts]).
   [reactants, products, catalysts] are all in List form
BA
=#
function semanticCheckingForEachTriplet(tripletArr::Array, conversionDict::Dict,
         error::Dict)
  # tripletArr = Array of ((senType, line#), BioSym, paraSet)
  preIRForm = []  # array of tuple ((verb, line#), reactant, product, catalyst)
  sys2userDict = Dict()
  for (key, val) in conversionDict
    sys2userDict[val] = key
  end

  for tp in tripletArr
    err_mess = []  # collection of error messages
    lineNum = tp[1][2]
    tmpIR = []
    if tp[1][1] == "react"
      tmpIR = checkReactType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "uptake"
      tmpIR = checkUptakeType(tp, conversionDict, sys2userDict, err_mess)
    elseif tp[1][1] == "secrete"
      tmpIR = checkSecreteType(tp, conversionDict, sys2userDict, err_mess)
    elseif tp[1][1] == "bind"
      tmpIR = checkBindType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "unbind"
      tmpIR = checkUnbindType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "phosphorylate"
      tmpIR = checkPhosphorylateType(tp, conversionDict, sys2userDict, err_mess)
    elseif tp[1][1] == "dephosphorylate"
      tmpIR = checkDephosphorylateType(tp, conversionDict, sys2userDict, err_mess)
    elseif tp[1][1] == "catalyze"
      tmpIR = checkCatalyzeType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "induce"  # GENE
      tmpIR = checkInduceType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "activate"  # PROTEIN
      tmpIR = checkInduceType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "repress"  # GENE
      tmpIR = checkRepressType(tp, conversionDict, err_mess)
    elseif tp[1][1] == "inhibit"  # PROTEIN
      tmpIR = checkRepressType(tp, conversionDict, err_mess)
    else
      println("**********unknown sentence type: $(tp[1])***************")
      push!(err_mess, "**********unknown sentence type: $(tp[1])***************")
    end
    if length(err_mess) != 0  # error occur
      if !haskey(error, lineNum)
        error[lineNum] = []
      end
      append!(error[lineNum], err_mess)
    else
      if length(tmpIR) != 0  # have something to return
        append!(preIRForm, tmpIR)
      end
    end
  end
  return preIRForm
end

#==========Functions For Each Type========================#

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkDephosphorylateType(senInfo::Tuple, u2sDict::Dict, s2uDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "[]"])
  legalTypeSet2 = Set(["PROTEIN"])
  legalTypeSet3 = Set(["SITE"])
  # check parameter setting
  rever = false
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "reversible"
      rever =  true
    end
  end
  if DS10SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2) # correct reactant
    if isassigned(senInfo[2], 3)  # have site
      if DS10SemanticChecking(senInfo[2][3], u2sDict, legalTypeSet3)  # correct SITE
        product = createProductsList(senInfo[2][2], "-"*senInfo[2][3][1].oriBioName, "", err)
      else  # incorrect site
        println("*semantic error or incorrect symbol format found in $(senInfo[2][3])")
        push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][3])")
        return []
      end
    else  # does not have site --> default product
      product = createProductsList(senInfo[2][2], s2uDict["-pho"], "", err)
    end
    if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1) # correct catalyst of DS1
      if rever  # reversible
        return [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]]), (senInfo[1], [product, senInfo[2][2], senInfo[2][1]])]
      else  # irreversible
        return [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]])]
      end
    elseif DS2SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1) # correct catalyst of DS2
      if rever  # reversible
        returnArr = []
        for cata in senInfo[2][1]
          append!(returnArr, [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]]),
            (senInfo[1], [product, senInfo[2][2], senInfo[2][1]])])
        end
        return returnArr
      else  # irreversible
        returnArr = [(senInfo[1], [senInfo[2][2], product, cata]) for cata in senInfo[2][1]]
        return returnArr
      end
    else  # incorrect catalyst
      println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
    end
  else  # incorrect reactant
    println("***semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkPhosphorylateType(senInfo::Tuple, u2sDict::Dict, s2uDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "[]"])
  legalTypeSet2 = Set(["PROTEIN"])
  legalTypeSet3 = Set(["SITE"])
  # check parameter setting
  rever = false
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "reversible"
      rever =  true
    end
  end
  if DS10SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2) # correct reactant
    if isassigned(senInfo[2], 3)  # have site
      if DS10SemanticChecking(senInfo[2][3], u2sDict, legalTypeSet3)  # correct SITE
        product = createProductsList(senInfo[2][2], "-"*senInfo[2][3][1].oriBioName)
      else  # incorrect site
        println("*semantic error or incorrect symbol format found in $(senInfo[2][3])")
        push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][3])")
        return []
      end
    else  # does not have site --> default product
      product = createProductsList(senInfo[2][2], s2uDict["-pho"])
    end
    if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1) # correct catalyst of DS1
      if rever  # reversible
        return [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]]), (senInfo[1], [product, senInfo[2][2], senInfo[2][1]])]
      else  # irreversible
        return [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]])]
      end
    elseif DS2SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1) # correct catalyst of DS2
      if rever  # reversible
        returnArr = []
        for cata in senInfo[2][1]
          append!(returnArr, [(senInfo[1], [senInfo[2][2], product, senInfo[2][1]]),
            (senInfo[1], [product, senInfo[2][2], senInfo[2][1]])])
        end
        return returnArr
      else  # irreversible
        returnArr = [(senInfo[1], [senInfo[2][2], product, cata]) for cata in senInfo[2][1]]
        return returnArr
      end
    else  # incorrect catalyst
      println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
    end
  else  # incorrect reactant
    println("***semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkCatalyzeType(senInfo::Tuple, u2sDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "[]"])
  legalTypeSet2 = Set(["METABOLITE", "PROTEIN", "GENE"])
  # check parameter setting
  rever = false
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "reversible"
      rever =  true
    end
  end
  if !isassigned(senInfo[2],2) || !isassigned(senInfo[2], 3)
    push!(err, "no enough inputs for \'catalyze\' type")
  else
    if (DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2) &&
      DS1SemanticChecking(senInfo[2][3], u2sDict, legalTypeSet2))  # reactant/product correct
      if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # catalyst correct
        if rever  # reversible
          returnArr = [(senInfo[1], [senInfo[2][2], senInfo[2][3], senInfo[2][1]]),
            (senInfo[1], [senInfo[2][3], senInfo[2][2], senInfo[2][1]])]
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [senInfo[2][2], senInfo[2][3], senInfo[2][1]])]
          return returnArr
        end
      elseif DS2SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # catalyst correct
        if rever  # reversible
          returnArr = []
          for cata in senInfo[2][1]
            append!(returnArr, [(senInfo[1], [senInfo[2][2], senInfo[2][3], cata]),
              (senInfo[1], [senInfo[2][3], senInfo[2][2], cata])])
          end
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [senInfo[2][2], senInfo[2][3], cata]) for cata in senInfo[2][1]]
          return returnArr
        end
      else  # incorrect catalyst
        println("*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
        push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
      end
    else  # incorrect reactants and/or products
      if !DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)
        println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
        push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      end
      if !DS1SemanticChecking(senInfo[2][3], u2sDict, legalTypeSet2)
        println("***semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][3])")
        push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][3])")
      end
    end
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkRepressType(senInfo::Tuple, u2sDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE"])
  # Induce or Activate
  if senInfo[1][1] == "repress"
    legalTypeSet2 = Set(["GENE"])
  else
    legalTypeSet2 = Set(["PROTEIN"])
  end
  if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # regulators correct
    if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # targets correct
      returnArr = [(senInfo[1], [senInfo[2][1], [tar]]) for tar in senInfo[2][2]]
      return returnArr
    else  # semantic error or incorrect targets
      println("*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    end
  elseif DS2SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # regulators correct
    if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # targets correct
      returnArr = [(senInfo[1], [regu, [tar]]) for regu in senInfo[2][1] for tar in senInfo[2][2]]
      return returnArr
    else  # semantic error or incorrect targets
      println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    end
  else  # semantic error or incorrect targets
    println("***semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkInduceType(senInfo::Tuple, u2sDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE"])
  # Induce or Activate
  if senInfo[1][1] == "induce"
    legalTypeSet2 = Set(["GENE"])
  else
    legalTypeSet2 = Set(["PROTEIN"])
  end
  if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # regulators correct
    if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # targets correct
      returnArr = [(senInfo[1], [senInfo[2][1], [tar]]) for tar in senInfo[2][2]]
      return returnArr
    else  # semantic error or incorrect targets
      println("*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    end
  elseif DS2SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # regulators correct
    if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # targets correct
      returnArr = [(senInfo[1], [regu, [tar]]) for regu in senInfo[2][1] for tar in senInfo[2][2]]
      return returnArr
    else  # semantic error or incorrect targets
      println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
    end
  else  # semantic error or incorrect targets
    println("***semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkUnbindType(senInfo::Tuple, u2sDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE", "MRNA"])
  legalTypeSet2 = Set(["PROTEIN", "METABOLITE", "MRNA", "GENE"])
  # check parameter setting
  rever = true
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "irreversible"
      rever =  false
    end
  end
  if DS10SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # reactants correct
    if isassigned(senInfo[2], 2)  # have products
      if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # correct products
        if rever  # reversible reaction
          returnArr = [(senInfo[1], [senInfo[2][1], senInfo[2][2]]), (senInfo[1], [senInfo[2][2], senInfo[2][1]])]
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [senInfo[2][1], senInfo[2][2]])]
          return returnArr
        end
      else  # incorrect product
        println("*semantic error or depreciated expression in $(senInfo[2][2])")
        push!(err, "*semantic error or depreciated expression in $(senInfo[2][2])")
      end
    else  # error: product required
      println("*for \'unbind\', products need to be specified")
      push!(err, "*for \'unbind\', products need to be specified")
    end
  else  # incorrect reactants
    println("**semantic error or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkBindType(senInfo::Tuple, u2sDict::Dict, err::Array)
  # semantic checking, can not all be GENE
  semantic1 = false
  if typeof(senInfo[2][1]) <: Array && typeof(senInfo[2][1][1]) <: generalBioSym
    for gb in senInfo[2][1]
      if u2sDict[gb.bioType] != "GENE"
        semantic1 = true
        break
      end
    end
  end
  legalTypeSet2 = Set(["PROTEIN", "METABOLITE", "MRNA"])
  # check parameter setting
  rever = true
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "irreversible"
      rever =  false
    end
  end
  if semantic1  # reactants correct
    if isassigned(senInfo[2], 2)  # have products
      if DS10SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # correct product
        if rever  # reversible reaction
          returnArr = [(senInfo[1], [senInfo[2][1], senInfo[2][2]]), (senInfo[1], [senInfo[2][2], senInfo[2][1]])]
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [senInfo[2][1], senInfo[2][2]])]
          return returnArr
        end
      else  # incorrect product
        println("*semantic error or depreciated expression in $(senInfo[2][2])")
        push!(err, "*semantic error or depreciated expression in $(senInfo[2][2])")
      end
    else  # generate product
      product = createBindProductsList(senInfo[2][1])
      if rever  # reversible reaction
        returnArr = [(senInfo[1], [senInfo[2][1], product]), (senInfo[1], [product, senInfo[2][1]])]
        return returnArr
      else  # irreversible
        returnArr = [(senInfo[1], [senInfo[2][1], product])]
        return returnArr
      end
    end
  else  # incorrect reactants
    println("**semantic error or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkSecreteType(senInfo::Tuple, u2sDict::Dict, s2uDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE"])
  legalTypeSet2 = Set(["PROTEIN"])
  # check parameter setting
  rever = false
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "reversible"
      rever =  true
    end
  end
  # semantics = DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)
  if DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # reactants correct
    products = createProductsList(senInfo[2][1], s2uDict["_c"], s2uDict["_exc"], err)
    if isassigned(senInfo[2], 2)  # have transporters
      if DS10SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # unambiguous
        # semantics2 = DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)
        if DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)  # transporters correct
          if rever  # reversible reaction
            returnArr = []
            for s in senInfo[2][2]
              append!(returnArr, [(senInfo[1], [senInfo[2][1], products, [s]]),
                (senInfo[1], [products, senInfo[2][1], [s]])])
            end
            return returnArr
          else  # irreversible
            returnArr = []
            for s in senInfo[2][2]
              append!(returnArr, [(senInfo[1], [senInfo[2][1], products, [s]])])
            end
            return returnArr
          end
        else  # transporter(s) not in the required form
          println("*semantic error or incorrect symbol format found in $(senInfo[2][2])")
          push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][2])")
        end
      else  # ambiguous
        println("**depreciated expression due to inherent ambiguity in $(senInfo[2][2])")
        push!(err, "*depreciated expression due to inherent ambiguity in $(senInfo[2][2])")
      end
    else  # does not have transporters
      returnArr = []
      if rever  # reversible
        for (id, sym) in enumerate(senInfo[2][1])
          append!(returnArr, [(senInfo[1], [[sym], [products[id]]]), (senInfo[1], [[products[id]], [sym]])])
        end
      else  # irreversible
        for (id, sym) in enumerate(senInfo[2][1])
          append!(returnArr, [(senInfo[1], [[sym], [products[id]]])])
        end
      end
      return returnArr
    end
  else  # incorrect reactants
    println("***semantic error or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table, reverse conversion table;
O: array of bi-tuple;
BA
=#
function checkUptakeType(senInfo::Tuple, u2sDict::Dict, s2uDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE", "MRNA", "GENE"])
  legalTypeSet2 = Set(["PROTEIN"])
  # check parameter setting
  rever = false
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "reversible"
      rever =  true
    end
  end
  semantics = DS1SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)
  if semantics  # reactants correct
    products = createProductsList(senInfo[2][1], s2uDict["_exc"], s2uDict["_c"], err)
    if isassigned(senInfo[2], 2)  # have transporters
      if DS10SemanticChecking(senInfo[2][1], u2sDict, legalTypeSet1)  # unambiguous
        semantics2 = DS1SemanticChecking(senInfo[2][2], u2sDict, legalTypeSet2)
        if semantics2  # transporters correct
          if rever  # reversible reaction
            returnArr = []
            for s in senInfo[2][2]
              append!(returnArr, [(senInfo[1], [senInfo[2][1], products, [s]]),
                (senInfo[1], [products, senInfo[2][1], [s]])])
            end
            return returnArr
          else  # irreversible
            returnArr = []
            for s in senInfo[2][2]
              append!(returnArr, [(senInfo[1], [senInfo[2][1], products, [s]])])
            end
            return returnArr
          end
        else  # transporter(s) not in the required form
          println("*semantic error or incorrect symbol format found in $(senInfo[2][2])")
          push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][2])")
        end
      else  # ambiguous
        println("**depreciated expression due to inherent ambiguity in $(senInfo[2][2])")
        push!(err, "*depreciated expression due to inherent ambiguity in $(senInfo[2][2])")
      end
    else  # does not have transporters
      returnArr = []
      if rever  # reversible
        for (id, sym) in enumerate(senInfo[2][1])
          append!(returnArr, [(senInfo[1], [[sym], [products[id]]]), (senInfo[1], [[products[id]], [sym]])])
        end
      else  # irreversible
        for (id, sym) in enumerate(senInfo[2][1])
          append!(returnArr, [(senInfo[1], [[sym], [products[id]]])])
        end
      end
      return returnArr
    end
  else  # incorrect reactants
    println("***semantic error or incorrect symbol format found in $(senInfo[2][1])")
    push!(err, "*semantic error or incorrect symbol format found in $(senInfo[2][1])")
  end
  return []
end

#=
F:
I: Tuple of biosymbols and parameters, conversion table;
O: array of bi-tuple;
BA
=#
function checkReactType(senInfo::Tuple, conDict::Dict, err::Array)
  legalTypeSet1 = Set(["PROTEIN", "METABOLITE", "MRNA", "GENE"])
  legalTypeSet2 = Set(["[]"])
  # check parameter setting
  rever = true
  if isdefined(senInfo, 3)  # has parameters, may have model parameters in future
    if senInfo[3][1] == "irreversible"
      rever =  false
    end
  end
  # must have symbols
  if length(senInfo[2]) < 2
    println("no enough input for sentence $(senInfo[1][1])")
    push!(err, "no enough input for sentence $(senInfo[1][1])")
  else
    if length(senInfo[2]) > 2
      println("too many inputs, only the first two groups are handled")
      push!(err, "too many inputs, only the first two groups are handled")
    end
    if DS10SemanticChecking(senInfo[2][1], conDict, legalTypeSet2)  # comes from nothing
      if DS1SemanticChecking(senInfo[2][2], conDict, legalTypeSet1)  # products correct
        if rever  # reversible
          returnArr = []
          for reac in senInfo[2][2]
            append!(returnArr, [(senInfo[1], [senInfo[2][1], [reac]]),
              (senInfo[1], [[reac], senInfo[2][1]])])
          end
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [senInfo[2][1], [reac]]) for reac in senInfo[2][2]]
          return returnArr
        end
      else  # product incorrect
        println("*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
        push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][2])")
      end
    elseif DS10SemanticChecking(senInfo[2][2], conDict, legalTypeSet2)  # goes to nothing
      if DS1SemanticChecking(senInfo[2][1], conDict, legalTypeSet1)  # products correct
        if rever  # reversible
          returnArr = []
          for reac in senInfo[2][1]
            append!(returnArr, [(senInfo[1], [[reac], senInfo[2][2]]),
              (senInfo[1], [senInfo[2][2], [reac]])])
          end
          return returnArr
        else  # irreversible
          returnArr = [(senInfo[1], [[reac], senInfo[2][2]]) for reac in senInfo[2][1]]
          return returnArr
        end
      else  # product incorrect
        println("**semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
        push!(err, "*semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1])")
      end
    else  # incorrect usage of "react" or error
      println("***incorrect use of \'react\', suggest using \'catalyze\';" *
        " OR semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1]) and $(senInfo[2][2])")
      push!(err, "***incorrect use of \'react\', suggest using \'catalyze\';" *
        " OR semantic error or ambiguity or incorrect symbol format found in $(senInfo[2][1]) and $(senInfo[2][2])")
    end
  end
  return []
end

#==========Create Products List=============================#

#=
F: take in an array, create binding compounds, link by hyphen
IOBA
=#
function createBindProductsList(oriList::Array)
  product = generalBioSym()
  product.coeff = 1.0
  product.bioType = oriList[1].bioType  # FIXME compound type -> mRNA-pro, DNA-pro
  product.oriBioName = join([s.oriBioName for s in oriList], "--")
  return [product]
end

#=
F: take in an array, create corresponding products list
IOBA
=#
function createProductsList(oriList::Array, oldTag, newTag, err)
  newList = []
  for arr in oriList
    new_arr = deepcopy(arr)
    if !occursin(oldTag, new_arr.oriBioName)
      println("ERROR: no \'$oldTag\' found in \'$(new_arr.oriBioName)\'")
      push!(err, "ERROR: no \'$oldTag\' found in \'$(new_arr.oriBioName)\'")
    end
    new_arr.oriBioName = replace(new_arr.oriBioName, oldTag => newTag)
    push!(newList, new_arr)
  end
  return newList
end

#=
F: take in an array, create corresponding products list
IOBA
=#
function createProductsList(oriList::Array, addTag::AbstractString)
  newList = []
  for arr in oriList
    new_arr = deepcopy(arr)
    new_arr.oriBioName = new_arr.oriBioName * addTag
    push!(newList, new_arr)
  end
  return newList
end

#==========Data Structure Checking=========================#

#=
F: check an input array is an array with a single/multiple element(s) which
  is/are of generalBioSym type
I:
O: semanticStatus: true if correct, otherwise, false;
BA
=#
function DS1SemanticChecking(inputArr::Array, typeDict::Dict, typeSet)
  if (DS10SemanticChecking(inputArr::Array, typeDict::Dict, typeSet) ||
    DS11SemanticChecking(inputArr::Array, typeDict::Dict, typeSet))
    return true
  else
    return false
  end
end

#=
F: check an input array is an array with a single element which is of generalBioSym type
I:
O: semanticStatus: true for correct, false for incorrect;
BA
=#
function DS10SemanticChecking(inputArr::Array, typeDict::Dict, typeSet)
  semanticStatus = true
  if typeof(inputArr) <: Array && length(inputArr) == 1 && typeof(inputArr[1]) <: generalBioSym
    for gb in inputArr
      if !in(typeDict[gb.bioType], typeSet)
        semanticStatus = false
        # println("semantic error found in $gb")
        break
      end
    end
  else  # incorrect format
    semanticStatus = false
    # println("incorrect symbol format for \'react\' type")
  end
  return semanticStatus
end

#=
F: check an input array is an array with many a single element which is of generalBioSym type
I:
O: semanticStatus: true for correct, false for incorrect;
BA
=#
function DS11SemanticChecking(inputArr::Array, typeDict::Dict, typeSet)
  semanticStatus = true
  if typeof(inputArr) <: Array && length(inputArr) > 1 && typeof(inputArr[1]) <: generalBioSym
    for gb in inputArr
      if !in(typeDict[gb.bioType], typeSet)
        semanticStatus = false
        # println("semantic error found in $gb")
        break
      end
    end
  else  # incorrect format
    semanticStatus = false
    # println("incorrect symbol format for \'react\' type")
  end
  return semanticStatus
end

#=
F: check an input array is an array of arrays with a single/multiple element(s)
   which is/are of generalBioSym type.
I:
O: semanticStatus: true if correct, otherwise, false;
BA
=#
function DS2SemanticChecking(inputArr::Array, typeDict::Dict, typeSet)
  semantic = false
  if (typeof(inputArr) <: Array) && (typeof(inputArr[1]) <: Array)  # data structure 2
    semantic = minimum([DS1SemanticChecking(s, typeDict, typeSet) for s in inputArr])
  end
  return semantic
end

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
module TestExteriorCalculus
using Test

using Catlab
using CombinatorialSpaces.ExteriorCalculus

@present Diffusion2DQuantities(FreeExtCalc2D) begin
  X::Space
  C::Hom(munit(), Form0(X))     # concentration
  ϕ::Hom(munit(), DualForm1(X)) # negative diffusion flux
  k::Hom(Form1(X), Form1(X))    # diffusivity (usually scalar multiplication)
end

@present Diffusion2D <: Diffusion2DQuantities begin
  # Fick's first law
  ϕ == C ⋅ d₀(X) ⋅ k ⋅ ⋆₁(X)
  # Diffusion equation
  C ⋅ ∂ₜ(Form0(X)) == ϕ ⋅ dual_d₁(X) ⋅ ⋆₀⁻¹(X)
end

@present Diffusion2D′ <: Diffusion2DQuantities begin
  # Diffusion equation
  C ⋅ ∂ₜ(Form0(X)) == C ⋅ d₀(X) ⋅ k ⋅ δ₁(X)
end

X, C, ϕ = Diffusion2D[:X], Diffusion2D[:C], Diffusion2D[:ϕ]
@test codom(C) == Form0(X)
@test codom(ϕ) == DualForm1(X)

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.SimulationEngine

Simulation engine for implicit index 1 DAE models with events.

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module SimulationEngine

export simulate!

import Sundials
import ModiaMath
import ModiaMath.DAE
using  StaticArrays

@static if VERSION >= v"0.7.0-DEV.2005"
    @eval using Printf
    using LinearAlgebra
end

# include code
# include("IDA_UserData.jl")  # is included within simulate.jl
include("simulate.jl")

end


text:
#******************************************************************************#
#                                                                              #
#                                                         :::      ::::::::    #
#    tools.jl                                           :+:      :+:    :+:    #
#                                                     +:+ +:+         +:+      #
#    By: mngomane <mngomane@student.42.fr>          +#+  +:+       +#+         #
#                                                 +#+#+#+#+#+   +#+            #
#    Created: 2015/02/12 08:28:05 by mngomane          #+#    #+#              #
#    Updated: 2015/02/12 08:28:05 by mngomane         ###   ########.fr        #
#                                                                              #
#******************************************************************************#

function	error_string(str)
	string(C_RED, str, C_END)
end

function	valid_string(str)
	string(C_GREEN, str, C_END)
end

function	prompt_string()
	string(C_PURPLE, "jshell ~> ", C_END)
end

function	error_prompt()
	string(C_RED, "[ERROR]~> ", C_END, "Available commands: start, exit")
end

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
using Documenter, NonlinearSolve

makedocs(
    sitename="NonlinearSolve.jl",
    authors="Chris Rackauckas",
    modules=[NonlinearSolve],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://nonlinearsolve.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorials/nonlinear.md",
            "tutorials/iterator_interface.md"
        ],
        "Basics" => Any[
            "basics/NonlinearProblem.md",
            "basics/NonlinearFunctions.md",
            "basics/FAQ.md"
        ],
        "Solvers" => Any[
            "solvers/NonlinearSystemSolvers.md",
            "solvers/BracketingSolvers.md"
        ]
    ]
)

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


text:
project_table = DataFrames.readtable("sample_metadata.txt", separator='\t')
pfam_table = apfamet.import_hmm_output(project_table)
pfam_stats_df = apfamet.prepare_pfam_database()
test_rpoB_equiv_table = apfamet.reads_to_rpob_equiv(pfam_table, pfam_stats_df)

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
using NPZ
using Plots
using UnitfulRecipes

include("line.jl")
include("functions.jl")
include("radiation.jl")
include("atmosphere.jl")
include("broadening.jl")
include("populations.jl")
include("characteristics.jl")

pyplot()

"""
    function circle_shape(x, y, r)

Function for a circle with centre coordinates (x, y) and radius r.
"""
function circle_shape(x, y, r)
    θ = LinRange(0, 2π, 500)
    x .+ r*cos.(θ), y .+ r*sin.(θ)
end

function plot_searchlight(k::Vector{Float64},
                          x,
                          y,
                          I,
                          R0,
                          title::String)

    @assert norm(k) ≈ 1 "Wrong direction vector!"

    RES = 300
    x_r = 0.5 - sign(k[1])*k[2]/k[1]
    if x_r < 0
        x_r = 1 - (ceil(x_r) - x_r)
    elseif x_r > 1
        x_r = x_r - floor(x_r)
    end

    y_r = 0.5 - sign(k[1])*k[3]/k[1]
    if y_r < 0
        y_r = 1 - (ceil(y_r) - y_r)
    elseif y_r > 1
        y_r = y_r - floor(y_r)
    end

    heatmap(ustrip.(x), ustrip.(y), ustrip.(transpose(I)),
            dpi=RES, title=title, xaxis="x", yaxis="y",
            aspect_ratio= :equal, clim=(0.,1.))
    plot!(circle_shape(ustrip(x_r), ustrip(y_r), ustrip(R0)),
          aspect_ratio = :equal,
          linecolor=:red,
          lw=2,
          label="")
    savefig("../img/compare_searchlight/$title")

end

function plot_top_intensity(atmos::Atmosphere,
                            line::HydrogenicLine,
                            S_λ::Array{<:UnitsIntensity_λ, 4},
                            α_tot::Array{<:PerLength, 4},
                            θ::Float64,
                            ϕ::Float64,
                            idλ::Int,
                            title::String)

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    I_top = short_characteristics_up(k, S_λ[idλ, :, :, :], S_λ[idλ, 1, :, :],
                                     α_tot[idλ, :, :, :], atmos)

    I_top = ustrip(uconvert.(u"kW*nm^-1*m^-2", I_top[end, 2:end-1, 2:end-1]))

    l1 = line.λidx[1]+1
    l2 = line.λidx[2]
    if l1 <= idλ <= l2
        colors = :gist_yarg
    else
        colors = :magma
    end

    heatmap(ustrip(atmos.x[2:end-1]),
            ustrip(atmos.y[2:end-1]),
            transpose(I_top),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title=title*" at $(round(ustrip(line.λ[idλ]); digits=4)) nm",
            aspect_ratio=:equal,
            c=colors,
            clim=(5,200))

    savefig("../img/compare_line/"*title)
end

function write_top_intensity(atmos::Atmosphere,
                             line::HydrogenicLine,
                             S_λ::Array{<:UnitsIntensity_λ, 4},
                             α_tot::Array{<:PerLength, 4},
                             θ::Float64,
                             ϕ::Float64,
                             fname::String)

    l1 = line.λidx[1]+1
    l2 = line.λidx[2]

    I_top = Array{Float64, 3}(undef, (l2-l1+1,
                                      length(atmos.x[2:end-1]),
                                      length(atmos.y[2:end-1])))

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    println("--Calculating intensity---")
    Threads.@threads for idλ in l1:l2
        intensity = short_characteristics_up(k, S_λ[idλ, :, :, :], S_λ[idλ, 1, :, :],
                                         α_tot[idλ, :, :, :], atmos)

        intensity = transpose(intensity[end, 2:end-1, 2:end-1])
        I_top[idλ,:,:] = ustrip(uconvert.(u"kW*nm^-1*m^-2", intensity))
    end

    println("---Writing to file---")
    npzwrite("../python/linedata/$(fname).npy", I_top)
end

function plot_top_line(atmos::Atmosphere,
                       line::HydrogenicLine,
                       S_λ::Array{<:UnitsIntensity_λ, 4},
                       α_tot::Array{<:PerLength, 4},
                       θ::Float64,
                       ϕ::Float64,
                       title::String)

    nλ, nz, nx, ny = size(S_λ)
    I_λ = zero(S_λ)

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    for iλ in 1:nλ
        I_λ[iλ, :, :, :] = short_characteristics_up(k, S_λ[iλ, :, :, :], S_λ[iλ, 1, :, :],
                                                    α_tot[iλ, :, :, :], atmos)
    end

    start = line.λidx[1]+1
    stop =  line.λidx[2]

    indices = sortperm(line.λ[start:stop])

    for idx in 1:10:size(S_λ)[end]
        for idy in 1:10:size(S_λ)[end-1]
            loc = "_"*string(idx)*"_"*string(idy)
            I_plot = ustrip(uconvert.(u"kW*nm^-1*m^-2", I_λ))[:, end, idx, idy]
            plot(ustrip.(line.λ[start:stop]),
                 ustrip.(I_plot[start:stop]),
                 xaxis="λ [nm]",
                 yaxis="Intensity [kW m^-2 nm^-1]",
                 dpi=300,
                 rightmargin=10Plots.mm,
                 title=title*loc)

            savefig("../img/compare_line/lines/"*title*loc)
        end
    end
end

function plot_top_cont(atmos::Atmosphere,
                       λ::Vector{<:Unitful.Length},
                       S_λ::Array{<:UnitsIntensity_λ, 4},
                       α_tot::Array{<:PerLength, 4},
                       θ::Float64,
                       ϕ::Float64,
                       title::String)

    nλ, nz, nx, ny = size(S_λ)
    I_λ = zero(S_λ)

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    for iλ in 1:nλ
        I_λ[iλ, :, :, :] = short_characteristics_up(k, S_λ[iλ, :, :, :], S_λ[iλ, 1, :, :],
                                                    α_tot[iλ, :, :, :], atmos)
    end

    for idx in 1:10:size(S_λ)[end]
        for idy in 1:10:size(S_λ)[end-1]
            loc = "_"*string(idx)*"_"*string(idy)
            I_plot = ustrip(uconvert.(u"kW*nm^-1*m^-2", I_λ))[:, end, idx, idy]
            plot(ustrip.(λ[:]),
                 ustrip.(I_plot[:]),
                 xaxis="λ [nm]",
                 yaxis="Intensity [kW m^-2 nm^-1]",
                 dpi=300,
                 rightmargin=10Plots.mm,
                 title=title*loc,
                 clim=(5,200))

            savefig("../img/compare_continuum/"*title*loc)
        end
    end
end

"""
    read_quantities(DATA::String; periodic=true)

read quantities from simulation from a hdf5 file
"""
function read_quantities(DATA::String; periodic=true)
    atmos = Atmosphere(get_atmos(DATA; periodic=periodic)...)
    global S_λ, populations
    h5open(DATA, "r") do file
        S_λ = read(file, "source_function")[:, :, :, :]*u"kW*m^-2*nm^-1"
        populations = read(file, "populations")[:, :, :, :]*u"m^-3"
    end

    if periodic
        S_λ = periodic_borders(S_λ)
        populations = periodic_pops(populations)
    end
    return atmos, S_λ, populations
end

"""
    read_irregular(DATA::String)

read quantities from irregular grid simulation from a hdf5 file
"""
function read_irregular(DATA::String)
    local positions, temperature, electron_density, hydrogen_populations, S_λ, populations
    local velocity_z, velocity_x, velocity_y, boundaries
    println("--Extracting simulation results---")
    h5open(DATA, "r") do file
        positions = read(file, "positions")[:, :]*u"m"
        boundaries = read(file, "boundaries")[:]u"m"

        temperature = read(file, "temperature")[:]*u"K"

        electron_density = read(file, "electron_density")[:]*u"m^-3"
        hydrogen_populations = read(file, "hydrogen_populations")[:]*u"m^-3"

        velocity_z = read(file, "velocity_z")[:]u"m*s^-1"
        velocity_x = read(file, "velocity_x")[:]u"m*s^-1"
        velocity_y = read(file, "velocity_y")[:]u"m*s^-1"

        S_λ = read(file, "source_function")[:, :]*u"kW*m^-2*nm^-1"
        populations = read(file, "populations")[:, :]*u"m^-3"
    end

    nn = Matrix{Int}(undef, (1, 1))
    ll = Vector{Int}(undef, 1)

    sites = VoronoiSites(positions, nn, ll, ll, ll, ll, temperature,
                         electron_density, hydrogen_populations, velocity_z,
                         velocity_x, velocity_y, boundaries...,
                         size(positions)[1])

    atmos_size = (430, 256, 256)
    atmos_size = floor.(Int, atmos_size.*1.0)

    println("--Converting grid---")
    atmos, S_λ_grid, populations_grid = Voronoi_to_Raster(sites, atmos_size,
                                                          S_λ, populations;
                                                          periodic=true)

    return atmos, S_λ_grid, populations_grid
end

"""
    plotter(atmos::Atmosphere,
            S_λ::Array{<:UnitsIntensity_λ, 4},
            populations::Array{<:NumberDensity, 4},
            θ::Float64,
            ϕ::Float64)

Compute simulations results for beam angle
"""
function plotter(atmos::Atmosphere,
                 S_λ::Array{<:UnitsIntensity_λ, 4},
                 populations::Array{<:NumberDensity, 4},
                 θ::Float64,
                 ϕ::Float64,
                 title::String)

    line = HydrogenicLine(test_atom(50, 20)..., atmos.temperature)

    LTE_pops = LTE_populations(line, atmos)

    # Find continuum extinction
    α_cont = α_absorption.(line.λ0,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(line.λ0,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])

    γ = γ_constant(line,
                   atmos.temperature,
                   (populations[:, :, :, 1].+populations[:, :, :, 2]),
                   atmos.electron_density)

    damping_λ = Array{Float64, 4}(undef, size(S_λ))
    Threads.@threads for l in eachindex(line.λ)
        damping_λ[l, :, :, :] = damping.(γ, line.λ[l], line.ΔD)
    end

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    profile = compute_voigt_profile(line, atmos, damping_λ, k)

    α_tot = Array{Float64, 4}(undef, size(profile))u"m^-1"
    Threads.@threads for l in eachindex(line.λ)
        α_tot[l, :, :, :] = αline_λ(line,
                                    profile[l, :, :, :],
                                    populations[:, :, :, 2],
                                    populations[:, :, :, 1]) .+ α_cont
    end

    return atmos, line, S_λ, α_tot
end

function plot_convergence(DATA::String, title::String)
    local convergence
    h5open(DATA, "r") do file
        convergence = read(file, "convergence")[:]
    end

    converged = argmin(convergence)
    # println(convergence[1:converged-1])
    plot(convergence[1:converged-1],
         xlabel="iteration",
         ylabel="max rel. diff.",
         title=title,
         dpi=300,
         yscale=:log10)
    savefig("../img/$title")
    # println(title)
end

function write_convergence(DATA::String, title::String)
    local convergence
    h5open(DATA, "r") do file
        convergence = read(file, "convergence")[:]
    end

    fname = split(DATA, "/")[end]
    fname = split(fname, ".")[1]

    converged = argmin(convergence)
    convergence=convergence[1:converged-1]

    npzwrite("../python/convergence/$(title)_convergence.npy", convergence)
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation.ParameterEstimation

using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.ParameterDistributionStorage

using Distributions
using LinearAlgebra
using Random
using ArgParse

#=
s = ArgParseSettings()
@add_arg_table! s begin
    "--relative_weight_option"
        default = "all_but_e"
        arg_type = String
    "--free_parameters"
        default = "CATKEParametersRiDependent"
        arg_type = String
end
args = parse_args(s)
relative_weight_option = args["relative_weight_option"]
free_parameter_type = args["free_parameters"]
=#

relative_weight_option = "all_but_e"
free_parameter_type = "CATKEParametersRiDependent"

relative_weight = Dict(:T => 1.0, :U => 0.5, :V => 0.5, :e => 0.0)

include("calibration_scripts/EKI_hyperparameter_search.jl")

directory = "EKI/$(free_parameter_type)_$(relative_weight_option)/"
isdir(directory) || mkpath(directory)

p = Parameters(RelevantParameters = free_parameter_options[free_parameter_option],
               ParametersToOptimize = free_parameter_options[free_parameter_option])

calibration = InverseProblem(FourDaySuite, p; relative_weights = relative_weight_options["all_but_e"], grid_type=OneDimensionalEnsembleGrid, Nz=64, Δt=10.0);
validation = InverseProblem(merge(TwoDaySuite, SixDaySuite), p; relative_weights = relative_weight_options["all_but_e"], grid_type=OneDimensionalEnsembleGrid, Nz=64, Δt=10.0);
ce = CalibrationExperiment(calibration, validation, p);

loss = ce.calibration.loss
loss_validation = ce.validation.loss
initial_parameters = ce.calibration.default_parameters
parameternames = propertynames(initial_parameters)

plot_stds_within_bounds(loss, loss_validation, initial_parameters, directory, xrange=-3:0.25:5)
v = plot_prior_variance(loss, loss_validation, initial_parameters, directory; xrange=0.1:0.05:1.0)
plot_num_ensemble_members(loss, loss_validation, initial_parameters, directory; xrange=1:5:30)
nl = plot_observation_noise_level(loss, loss_validation, initial_parameters, directory; xrange=-2.0:0.1:3.0)


text:
export COFFOptionalHeaderStandard, COFFImageDataDirectory, COFFDataDirectories,
       COFFOptionalHeaderWindows32, COFFOptionalHeaderWindows64,
       COFFOptionalHeader32, COFFOptionalHeader64

@io struct COFFOptionalHeaderStandard
    Magic::UInt16
    MajorLinkerVersion::UInt8
    MinorLinkerVersion::UInt8
    SizeOfCode::UInt32
    SizeOfInitializedData::UInt32
    SizeOfUninitializedData::UInt32
    AddressOfEntryPoint::UInt32
    BaseOfCode::UInt32
end
const OPTHEADER_STANDARD_MAGIC32 = 0x10b
const OPTHEADER_STANDARD_MAGIC64 = 0x20b

@io struct COFFImageDataDirectory
    VirtualAddress::UInt32
    Size::UInt32
end

@io struct COFFDataDirectories
    ExportTable::COFFImageDataDirectory
    ImportTable::COFFImageDataDirectory
    ResourceTable::COFFImageDataDirectory
    ExceptionTable::COFFImageDataDirectory
    CertificateTable::COFFImageDataDirectory
    BaseRelocatioNTable::COFFImageDataDirectory
    Debug::COFFImageDataDirectory
    Architecture::COFFImageDataDirectory
    GlobalPtr::COFFImageDataDirectory
    TLSTable::COFFImageDataDirectory
    LoadConfigTable::COFFImageDataDirectory
    BoundImport::COFFImageDataDirectory
    IAT::COFFImageDataDirectory
    DelayImportDescriptor::COFFImageDataDirectory
    CLRRuntimeHeader::COFFImageDataDirectory
    Reserverd::COFFImageDataDirectory
end

@io struct COFFOptionalHeaderWindows32
    ImageBase::UInt32
    SectionAlignment::UInt32
    FileAlignment::UInt32
    MajorOperatingSystemVersion::UInt16
    MinorOperatingSystemVersion::UInt16
    MajorImageVersion::UInt16
    MinorImageVersion::UInt16
    MajorSubsystemVersion::UInt16
    MinorSubsystemVersion::UInt16
    Win32VersionValue::UInt32
    SizeOfImage::UInt32
    SizeOfHeaders::UInt32
    CheckSum::UInt32
    Subsystem::UInt16
    DllCharacteristics::UInt16
    SizeOfStackReserve::UInt32
    SizeOfStackCommit::UInt32
    SizeOfHeapReserve::UInt32
    SizeOfHeapCommit::UInt32
    LoaderFlag::UInt32
    NumberOfRvaAndSizes::UInt32
end

@io struct COFFOptionalHeaderWindows64
    ImageBase::UInt64
    SectionAlignment::UInt32
    FileAlignment::UInt32
    MajorOperatingSystemVersion::UInt16
    MinorOperatingSystemVersion::UInt16
    MajorImageVersion::UInt16
    MinorImageVersion::UInt16
    MajorSubsystemVersion::UInt16
    MinorSubsystemVersion::UInt16
    Win32VersionValue::UInt32
    SizeOfImage::UInt32
    SizeOfHeaders::UInt32
    CheckSum::UInt32
    Subsystem::UInt16
    DllCharacteristics::UInt16
    SizeOfStackReserve::UInt64
    SizeOfStackCommit::UInt64
    SizeOfHeapReserve::UInt64
    SizeOfHeapCommit::UInt64
    LoaderFlag::UInt32
    NumberOfRvaAndSizes::UInt32
end

abstract type COFFOptionalHeader end
@io struct COFFOptionalHeader32 <: COFFOptionalHeader
    standard::COFFOptionalHeaderStandard
    BaseOfData::UInt32
    windows::COFFOptionalHeaderWindows32
    directories::COFFDataDirectories
end

@io struct COFFOptionalHeader64 <: COFFOptionalHeader
    standard::COFFOptionalHeaderStandard
    windows::COFFOptionalHeaderWindows64
    directories::COFFDataDirectories
end

function read(io::IO, ::Type{COFFOptionalHeader})
    # First, unpack the standard bits
    standard = unpack(io, COFFOptionalHeaderStandard)

    # If it's a 32-bit header, read that out, otherwise, do the 64-bit dance
    if standard.Magic == OPTHEADER_STANDARD_MAGIC32
        BaseOfData = read(io, UInt32)
        windows = unpack(io, COFFOptionalHeaderWindows32)
        dirs = unpack(io, COFFDataDirectories)
        return COFFOptionalHeader32(standard, BaseOfData, windows, dirs)
    elseif standard.Magic == OPTHEADER_STANDARD_MAGIC64
        windows = unpack(io, COFFOptionalHeaderWindows64)
        dirs = unpack(io, COFFDataDirectories)
        return COFFOptionalHeader64(standard, windows, dirs)
    else
        error("Unknown COFF optional header magic 0x$(hex(standard.Magic))")
    end
end


text:
using MTH229

using Test

## test package
@testset "test packages" begin

    ## Roots
    @test fzero(sin, 3, 4)  ≈ pi
    @test fzero(sin, 3.0)  ≈ pi

    ## Zygote
    f(x) = sin(x)
    @test f'(2)  ≈ cos(2)
    @test f''(2)  ≈ -sin(2)

    ## QuadGK
    @test quadgk(sin, 0, pi)[1]  ≈ 2
    @test abs(quadgk(x -> x^2, 0, 1)[2]) <= 1e-14

    ## SymPy
    @test limit(x -> sin(x)/x, 0) == 1
    @test integrate(sin,  0, pi) == 2

end

@testset "test functions" begin

    f(x) = sin(x)
    c = pi/4
    fn = tangent(f, c)
    @test fn(1)  ≈ f(c) + f'(c)*(1 - c)

    fn = secant(f, pi/6, pi/3)
    @test fn(pi/4) <= f(pi/4)

    out = lim(x -> sin(x)/x, 0)
    @test out[end, 2]  ≈ 1

    @test bisection(sin, 3, 4)  ≈ pi
    @test newton(sin, cos, 3.0)  ≈ pi
    @test newton(sin, 3.0)  ≈ pi

    f₁ = x ->  1 + 100x^2 - x^3
    @test isempty(find_zeros(f₁, -100, 100) )
    f₁s = find_zeros(fisheye(f₁), -pi/2, pi/2)
    @test !isempty(f₁s)
    @test abs(f₁(tan(first(f₁s)))) <= 1e-4

    out = sign_chart(x -> (x-1)*(x-2)/(x-3), 0, 4)
    @test all([o[1] for o ∈ out] .≈[1,2,3])

    @test riemann(sin, 0, pi, 10_000)  ≈ 2
    @test fubini((x,y) -> 1, (x->-sqrt(1-x^2), x->sqrt(1-x^2)), (-1,1)) ≈ pi
end

@testset "2d" begin

    x = [[1,2,3], [4,5,6]]
    @test unzip(x)[1] == [1, 4]
    @test unzip(x)[2] == [2, 5]
    @test unzip(x)[3] == [3, 6]

    @test length(unzip(x -> x, 0, 1)[1])  <= 50 # 21
    @test length(unzip(x-> sin(10pi*x), 0, 1)[1]) >= 50 # 233
    @test length(unzip(x -> x, 0, 1, 10)[1])  == 10

    @test uvec([2,2]) == 1/sqrt(2) * [1,1]

end

@testset "plots" begin

    if isinteractive()

        f(x) = 1/x

        trimplot(f, -1, 1)

        plotif(f, f, -1, 1)
        plotif(f, f', -1, 1)

        signchart(sin, 0, 4pi)

        newton_vis(log, 1/2)

        r(t) = [sin(t), cos(t)]
        ts = range(0, stop=pi/2, length=100)
        plot(unzip(r.(ts))...)
        arrow!(r(pi/4), r'(pi/4))

        V(x,y) = [x, x-y]
        plot()
        vectorfieldplot!(V)

    end

end


text:
#! /usr/bin/julia

# Rosetta Code, Grayscale image

using Color, Images, FixedPointNumbers

const M_RGB_Y = reshape(Color.M_RGB_XYZ[2,:], 3)

function rgb2gray(img::Image)
    g = red(img)*M_RGB_Y[1] + green(img)*M_RGB_Y[2] + blue(img)*M_RGB_Y[3]
    g = clamp(g, 0.0, 1.0)
    return grayim(g)
end

function gray2rgb(img::Image)
    colorspace(img) == "Gray" || return img
    g = map((x)->RGB{Ufixed8}(x, x, x), img.data)
    return Image(g, spatialorder=spatialorder(img))
end

ima = imread("grayscale_image_color.png")
imb = rgb2gray(ima)
imc = gray2rgb(imb)
imwrite(imc, "grayscale_image_rc.png")

ima = imread("grayscale_image_color.png")
imb = convert(Image{Gray{Ufixed8}}, ima)
imwrite(imb, "grayscale_image_julia.png")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
Parts of a single season

# Examples
```jldoctest
julia> PRE
PRE::SeasonPart = 0

julia> REG
REG::SeasonPart = 1

julia> POST
POST::SeasonPart = 2
```
"""
@enum SeasonPart PRE REG POST

import Base.parse

"""
    parse(::Type{SeasonPart}, str::AbstractString)

Parse `str` as a part of a season.

# Examples
```jldoctest
julia> parse(SeasonPart, "PRE")
PRE::SeasonPart = 0

julia> parse(SeasonPart, "REG")
REG::SeasonPart = 1

julia> parse(SeasonPart, "POST")
POST::SeasonPart = 2
```
"""
function Base.parse(::Type{SeasonPart}, str::AbstractString)
    m = Dict([string(e) => e for e in instances(SeasonPart)]...)
    return m[str]
end

"""
An enumeration of the SuperBowls as stylized/marketed (with the season as the value)

# Examples
```jldoctest
julia> SuperBowl(2001)
SB_XXXVI::SuperBowl = 2001

julia> SuperBowl(2015)
SB_50::SuperBowl = 2015

```
"""
@enum SuperBowl begin
    SB_I       = 1966
    SB_II      = 1967
    SB_III     = 1968
    SB_IV      = 1969
    SB_V       = 1970
    SB_VI      = 1971
    SB_VII     = 1972
    SB_VIII    = 1973
    SB_IX      = 1974
    SB_X       = 1975
    SB_XI      = 1976
    SB_XII     = 1977
    SB_XIII    = 1978
    SB_XIV     = 1979
    SB_XV      = 1980
    SB_XVI     = 1981
    SB_XVII    = 1982
    SB_XVIII   = 1983
    SB_XIX     = 1984
    SB_XX      = 1985
    SB_XXI     = 1986
    SB_XXII    = 1987
    SB_XXIII   = 1988
    SB_XXIV    = 1989
    SB_XXV     = 1990
    SB_XXVI    = 1991
    SB_XXVII   = 1992
    SB_XXVIII  = 1993
    SB_XXIX    = 1994
    SB_XXX     = 1995
    SB_XXXI    = 1996
    SB_XXXII   = 1997
    SB_XXXIII  = 1998
    SB_XXXIV   = 1999
    SB_XXXV    = 2000
    SB_XXXVI   = 2001
    SB_XXXVII  = 2002
    SB_XXXVIII = 2003
    SB_XXXIX   = 2004
    SB_XL      = 2005
    SB_XLI     = 2006
    SB_XLII    = 2007
    SB_XLIII   = 2008
    SB_XLIV    = 2009
    SB_XLV     = 2010
    SB_XLVI    = 2011
    SB_XLVII   = 2012
    SB_XLVIII  = 2013
    SB_XLIX    = 2014
    SB_50      = 2015
    SB_LI      = 2016
    SB_LII     = 2017
    SB_LIII    = 2018
    SB_LIV     = 2019
    SB_LV      = 2020
    SB_LVI     = 2021
    # SB_LVII    = 2022
    # SB_LVIII   = 2023
end


text:
function hopping_chain(lx         ::Int,
                       t          ::T=1.0,
                       mu         ::T=0.0;
                       boundary   ::Symbol=:OBC) where {T<:RLorCX}

    hopmatrix = diagm(0 => mu .* ones(T, lx)) +
        diagm(1 => -t .* ones(T, lx-1)) +
        diagm(-1 => conj(-t) .* ones(T, lx-1))

    if boundary == :OBC
        return hopmatrix
    elseif boundary == :PBC
        hopmatrix[1, lx] = -t
        hopmatrix[lx, 1] = conj(-t)
        return hopmatrix
    elseif boundary == :APBC
        hopmatrix[1, lx] = t
        hopmatrix[lx, 1] = conj(t)
        return hopmatrix
    else
        error("unrecognized boundary condition : ", boundary)
    end
    hopmatrix
end

function nnhoppingchain(lx::Int, t1::T, t2::T, mu::T=0.0;
                        boundary::Symbol=:OBC) where {T<:RLorCX}

    hopmatrix = diagm(0 => mu .* ones(T, lx)) +
        diagm(1 => -t1 .* ones(T, lx-1)) +
        diagm(-1 => conj(-t1) .* ones(T, lx-1)) +
        diagm(2 => -t2 .* ones(T, lx-2)) +
        diagm(-2 => conj(-t2) .* ones(T, lx-2))

    if boundary == :OBC
        return hopmatrix
    elseif boundary == :PBC
        hopmatrix[1, lx] = -t1
        hopmatrix[1, lx-1] = -t2
        hopmatrix[2, lx] = -t2
        hopmatrix[lx, 1] = conj(-t1)
        hopmatrix[lx-1, 1] = conj(-t2)
        hopmatrix[lx, 2] = conj(-t2)
        return hopmatrix

    elseif boundary == :APBC
        hopmatrix[1, lx] = t1
        hopmatrix[1, lx-1] = t2
        hopmatrix[2, lx] = t2
        hopmatrix[lx, 1] = conj(t1)
        hopmatrix[lx-1, 1] = conj(t2)
        hopmatrix[lx, 2] = conj(t2)
        return hopmatrix
    else
        error("unrecognized boundary condition : ", boundary)
    end
    hopmatrix
end

function generatebdg(model::UnitCellQModel)
    typeof(model.qtype) <: Fermion ||
        error("BdG Hamiltonian only for FermionType models!")

    ft, f = fermionoperators(model.qtype)
    n_sites = prod(model.lattice.sizes) * model.lattice.unitc.n
    H = zeros(eltype(model.inters[1]), n_sites, n_sites)
    for is in Iterators.product([1:l for l in model.lattice.sizes]...)
        for interaction in model.inters
            if support(interaction) > 2
                println("dropping higher order term")
                continue
            end
            ns = interaction.ucidxs
            offs = interaction.offsets
            indexes = zeros(Int, support(interaction))
            crosses = []
            isinside = true
            for i in 1:support(interaction)
                index, crossings = sitelinearindex(model.lattice, ns[i], offs[i] .+ is)
                if isnothing(index)
                    isinside = false
                    break
                end
                indexes[i] = index
                push!(crosses, crossings)
            end
            # only two-fermionic terms
            if isinside && length(indexes) < 3
                sign = +1
                for n = 1:dimension(model)
                    if model.lattice.bcs[n] == :APBC
                        #println("$indexes, $crosses, $(interaction.amp)")
                        if isodd(crosses[2][n] - crosses[1][n])
                            sign *= -1
                        end
                    end
                end
                for term in interaction.terms
                    if term == (ft, f)
                        i, j = indexes
                        H[i, j] += sign * interaction.amp
                    elseif term == (f, ft)
                        j, i  = indexes
                        H[i, j] += sign * conj(interaction.amp)
                    else
                        error("term not supported yet!")
                    end
                end
            end
        end
    end
    H
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
"""
Contains all functions and structures that require `DifferentialEquations` solvers.
This module provides orbit propagation, and Halo orbit solvers.
"""
module Propagators

# Module Exports
export propagate, halo, monodromy, isperiodic, potential_energy_hessian
export stable_eigenvector, unstable_eigenvector, manifold
export stable_manifold, unstable_manifold
export ODEProblem
export R2BPTic!, CR3BPTic!, CR3BPSTMTic!

# Module Dependencies
using Reexport

using Distributed
using StaticArrays
using LinearAlgebra
using SymbolicUtils
using ComponentArrays
using DifferentialEquations
using Unitful, UnitfulAstro, UnitfulAngles
using AstrodynamicalModels
using Distributed

using ..AstrodynamicsCore

include("R2BP/R2BPPropagators.jl")
include("CR3BP/CR3BPPropagators.jl")
include("CR3BP/Halos.jl")

end


text:
using BinDeps

@BinDeps.setup

include("constants.jl")
include("compile.jl")

# info("libname = $libname")
blas = library_dependency("libblas", alias=["libblas.dll"])
lapack = library_dependency("liblapack", alias=["liblapack.dll"])
depends = JULIA_LAPACK ? [] : [blas, lapack]

# LaPack/BLAS dependencies
if !JULIA_LAPACK
    @static if is_windows()
        # wheel = "numpy/windows-wheel-builder/raw/master/atlas-builds"
        # atlas = "https://github.com/$wheel"
        # atlasdll = "/atlas-3.11.38-sse2-64/lib/numpy-atlas.dll"
        # download("https://raw.githubusercontent.com/$wheel/$atlasdll"),
        #           "$libdir/libatlas.dll")
        ## at the end ...
        # push!(BinDeps.defaults, BuildProcess)
    end
end

csdp = library_dependency("csdp", aliases=[libname], depends=depends)

provides(Sources, URI(download_url), csdp, unpacked_dir=csdpversion)

provides(BuildProcess,
         (@build_steps begin
             GetSources(csdp)
             CreateDirectory(libdir)
             CreateDirectory(builddir)
             @build_steps begin
                  ChangeDirectory(srcdir)
                  patch_int
                  compile_objs
             end
         end),
         [csdp])

# Prebuilt DLLs for Windows
provides(Binaries,
   URI("https://github.com/EQt/winlapack/blob/master/win-csdp-$(Sys.WORD_SIZE).7z?raw=true"),
   [csdp, lapack, blas], unpacked_dir="usr", os = :Windows)

@BinDeps.install Dict(:csdp => :csdp)


text:
export TR_step_computation

function TR_step_computation(h :: Float64,
                             g :: Float64,
                             dN :: Float64,
                             Δ :: Float64)

  if h > 0.0
    if g > 0.0
      d = max(-Δ, dN)
    else
      d = min(dN, Δ)
    end
  else
    if g > 0.0
      d = -Δ
    else
      d = Δ
    end
  end

  return d

end


text:
#
# A value vector indexed by a vector of offsets.
#

# Mutable view over column slices.

type ColumnCursor{T,V<:AbstractVector} <: AbstractVector{T}
    pos::Int
    l::Int
    r::Int
    vals::V

    @inline ColumnCursor(vals) =
        new(0, 1, 1, vals)

    @inline function ColumnCursor(pos, l, r, vals)
        @boundscheck checkbounds(vals, l:r-1)
        new(pos, l, r, vals)
    end
end

ColumnCursor(vals::AbstractVector) =
    ColumnCursor{eltype(vals),typeof(vals)}(vals)

ColumnCursor(pos, l, r, vals::AbstractVector) =
    ColumnCursor{eltype(vals),typeof(vals)}(pos, l, r, vals)

@inline size(cr::ColumnCursor) = (cr.r - cr.l,)

@inline length(cr::ColumnCursor) = cr.r - cr.l

@inline function getindex(cr::ColumnCursor, i::Int)
    @boundscheck checkbounds(cr, i)
    @inbounds val = cr.vals[cr.l + i - 1]
    val
end

@inline function setindex!(cr::ColumnCursor, val, i::Int)
    @boundscheck checkbounds(cr, i)
    @inbounds cr.vals[cr.l + i - 1] = val
    cr
end

linearindexing{CR<:ColumnCursor}(::Type{CR}) = Base.LinearFast()

Base.array_eltype_show_how(::ColumnCursor) = (true, "")

# The column type.

immutable Column{OPT,PLU,O<:AbstractVector{Int},V<:AbstractVector} <: AbstractVector{ColumnCursor}
    len::Int
    offs::O
    vals::V

    function Column(offs::O, vals::V)
        len = length(offs)-1
        return new(len, offs, vals)
    end
end

typealias PlainColumn Column{false,false}
typealias OptionalColumn Column{true,false}
typealias PluralColumn Column{true,true}
typealias NonEmptyPluralColumn Column{false,true}

@inline function Column(offs::OneTo, vals::AbstractVector)
    @boundscheck length(offs) == length(vals) + 1 || error("ill-formed offset vector")
    Column{false,false,typeof(offs),typeof(vals)}(offs, vals)
end

function Column(offs::AbstractVector{Int}, vals::AbstractVector)
    optional = false
    plural = false
    @boundscheck !isempty(offs) && offs[1] == 1 && offs[end] == length(vals)+1 || error("ill-formed offset vector")
    l = 1
    for k = 2:endof(offs)
        @inbounds r = offs[k]
        @boundscheck l <= r || error("ill-formed offset vector")
        optional |= l == r
        plural |= l+1 < r
    end
    Column{optional,plural}(offs, vals)
end

@inline PlainColumn(vals::AbstractVector) =
    let offs = OneTo(length(vals)+1)
        PlainColumn(offs, vals)
    end

@inline PlainColumn(offs::OneTo, vals::AbstractVector) =
    Column{false,false,typeof(offs),typeof(vals)}(offs, vals)

@inline PlainColumn(::AbstractVector{Int}, vals::AbstractVector) =
    PlainColumn(vals)

@inline OptionalColumn(offs::AbstractVector{Int}, vals::AbstractVector) =
    Column{true,false,typeof(offs),typeof(vals)}(offs, vals)

@inline PluralColumn(offs::AbstractVector{Int}, vals::AbstractVector) =
    Column{true,true,typeof(offs),typeof(vals)}(offs, vals)

@inline NonEmptyPluralColumn(offs::AbstractVector{Int}, vals::AbstractVector) =
    Column{false,true,typeof(offs),typeof(vals)}(offs, vals)

convert(::Type{Column}, col::Column) = col

convert(::Type{Column}, data::AbstractVector) =
    let offs = OneTo(length(data)+1)
        Column{false,false}(offs, data)
    end

function convert{T<:Nullable}(::Type{Column}, data::Vector{T})
    len = 0
    for item in data
        if !isnull(item)
            len += 1
        end
    end
    offs = Vector{Int}(length(data)+1)
    vals = Vector{eltype(T)}(len)
    offs[1] = k = 1
    @inbounds for (i, item) in enumerate(data)
        if !isnull(item)
            vals[k] = get(item)
            k += 1
        end
        offs[i+1] = k
    end
    Column{true,false}(offs, vals)
end

function convert{T<:Vector}(::Type{Column}, data::Vector{T})
    len = 0
    for items in data
        len += length(items)
    end
    offs = Vector{Int}(length(data)+1)
    vals = Vector{eltype(T)}(len)
    offs[1] = k = 1
    @inbounds for (i, items) in enumerate(data)
        copy!(vals, k, items)
        k += length(items)
        offs[i+1] = k
    end
    Column{true,true}(offs, vals)
end

# Cursor interface.

@inline cursor(col::Column) =
    ColumnCursor(col.vals)

@inline function cursor(col::Column, pos::Int)
    @boundscheck checkbounds(col.offs, pos:pos+1)
    ColumnCursor(pos, col.offs[pos], col.offs[pos+1], col.vals)
end

@inline function move!(col::Column, cr::ColumnCursor, pos::Int)
    @boundscheck checkbounds(col.offs, pos:pos+1)
    cr.pos = pos
    @inbounds cr.l = col.offs[pos]
    @inbounds cr.r = col.offs[pos+1]
end

@inline done(col::Column, cr::ColumnCursor) = cr.pos >= col.len

@inline function next!(col::Column, cr::ColumnCursor)
    @boundscheck checkbounds(col.offs, cr.pos+1)
    cr.pos += 1
    cr.l = cr.r
    @inbounds cr.r = col.offs[cr.pos+1]
    cr
end

# Array interface.

eltype{OPT,PLU,O,V}(::Type{Column{OPT,PLU,O,V}}) = ColumnCursor{eltype(V),V}

@inline size(col::Column) = (col.len,)

@inline length(col::Column) = col.len

@inline getindex(col::Column, i::Int) =
    cursor(col, i)

function getindex{OPT,PLU,O,V}(col::Column{OPT,PLU,O,V}, idxs::AbstractVector{Int})
    @boundscheck checkbounds(col, idxs)
    offs = Vector{Int}(length(idxs)+1)
    @inbounds offs[1] = top = 1
    i = 1
    @inbounds for idx in idxs
        l = col.offs[idx]
        r = col.offs[idx+1]
        offs[i+1] = top = top + r - l
        i += 1
    end
    perm = Vector{Int}(offs[end]-1)
    i = 1
    @inbounds for idx in idxs
        l = col.offs[idx]
        r = col.offs[idx+1]
        copy!(perm, i, l:r-1)
        i += r - l
    end
    @inbounds vals = col.vals[perm]
    return Column{OPT,PLU}(offs, vals)
end

function getindex{OPT,PLU,O<:OneTo,V}(col::Column{OPT,PLU,O,V}, idxs::AbstractVector{Int})
    offs = OneTo(length(idxs)+1)
    vals = col.vals[idxs]
    return Column{OPT,PLU}(offs, vals)
end

function getindex{OPT,PLU,O,V}(col::Column{OPT,PLU,O,V}, idxs::OneTo)
    len = length(idxs)
    if len == col.len
        return col
    else
        return typeof(col)(
            col.offs[OneTo(len+1)],
            col.vals[OneTo(col.offs[len+1]-1)])
    end
end

function getindex{OPT,PLU,O<:OneTo,V}(col::Column{OPT,PLU,O,V}, idxs::OneTo)
    len = length(idxs)
    if len == col.len
        return col
    else
        return typeof(col)(col.offs[idxs], col.vals[idxs])
    end
end

function vcat{OPT1,OPT2,PLU1,PLU2}(col1::Column{OPT1,PLU1}, col2::Column{OPT2,PLU2})
    offs = Vector{Int}(col1.len+col2.len+1)
    for k = 1:col1.len
        offs[k] = col1.offs[k]
    end
    L = col1.offs[col1.len+1]
    for k = 1:col2.len+1
        offs[col1.len+k] = col2.offs[k] + L
    end
    vals = vcat(col1.vals, col2.vals)
    return Column{OPT1||OPT2,PLU1||PLU2}(offs, vals)
end

vcat{OPT1,OPT2,PLU1,PLU2}(col1::Column{OPT1,PLU1,OneTo{Int}}, col2::Column{OPT2,PLU2,OneTo{Int}}) =
    let offs = OneTo(col1.len+col2.len+1), vals = vcat(col1.vals, col2.vals)
        Column{OPT1||OPT2,PLU1||PLU2}(offs, vals)
    end

linearindexing{C<:Column}(::Type{C}) = Base.LinearFast()

Base.array_eltype_show_how(::Column) = (true, "")

Base.summary(col::Column) =
"$(col.len)-element column of $(typesummary(col))"

typesummary(col::Column) =
    "$(typesummary(col.vals))$(valsigil(col))"

typesummary(vals::AbstractVector) =
    "$(eltype(vals))"

summary{OPT,PLU,O,V}(col::Column{OPT,PLU,O,V}) =
    "$(shapesummary(valtype(col)))$(valsigil(col))"

# Attributes and properties.

valsigil{C<:Column{false,false}}(::Type{C}) = ""
valsigil{C<:Column{true,false}}(::Type{C}) = "?"
valsigil{C<:Column{false,true}}(::Type{C}) = "+"
valsigil{C<:Column{true,true}}(::Type{C}) = "*"
valsigil(col::Column) = valsigil(typeof(col))

Base.keytype{OPT,PLU,O,V}(::Type{Column{OPT,PLU,O,V}}) = eltype(O)
Base.valtype{OPT,PLU,O,V}(::Type{Column{OPT,PLU,O,V}}) = eltype(V)
Base.keytype(col::Column) = keytype(typeof(col))
Base.valtype(col::Column) = valtype(typeof(col))

isplain{C<:Union{Column{false,true},Column{true,false},Column{true,true}}}(::Type{C}) = false
isplain{C<:Column{false,false}}(::Type{C}) = true
isoptional{C<:Union{Column{false,false},Column{false,true}}}(::Type{C}) = false
isoptional{C<:Union{Column{true,false},Column{true,true}}}(::Type{C}) = true
isplural{C<:Union{Column{false,false},Column{true,false}}}(::Type{C}) = false
isplural{C<:Union{Column{false,true},Column{true,true}}}(::Type{C}) = true
isplain(col::Column) = isplain(typeof(col))
isoptional(col::Column) = isoptional(typeof(col))
isplural(col::Column) = isplural(typeof(col))

offsets(col::Column) = col.offs
values(col::Column) = col.vals

# Data view.

typealias ColumnTuple Tuple{Vararg{Column}}
typealias ColumnVector AbstractVector{Column}
typealias Columns Union{ColumnTuple,ColumnVector}

datatype(col::Column{false,false}) =
    eltype(col.vals)

datatype(col::Column{true,false}) =
    Nullable{eltype(col.vals)}

datatype{OPT}(col::Column{OPT,true}) =
    SubArray{eltype(col.vals),1,typeof(col.vals),Tuple{UnitRange{Int}},true}

datatype(cols::Columns) =
    Tuple{map(datatype, cols)...}

function getdata{T}(::Type{T}, col::Column, ::Void)
    @assert col.len == 1
    return getdata(T, col, 1)
end

@inline function getdata{OPT}(::Type{Any}, col::Column{OPT,true}, ::Void)
    @assert col.len == 1
    return col.vals
end

@inline function getdata{T,OPT}(::Type{T}, col::Column{OPT,true}, ::Void)
    return convert(T, col.vals)
end

@inline function getdata(::Type{Any}, col::Column{false,false}, i::Int)
    col.vals[i]
end

@inline function getdata(::Type{Any}, col::Column{true,false}, i::Int)
    l = col.offs[i]
    r = col.offs[i+1]
    T = Nullable{eltype(col.vals)}
    l < r ? convert(T, col.vals[l]) : convert(T, nothing)
end

@inline function getdata{OPT}(::Type{Any}, col::Column{OPT,true}, i::Int)
    l = col.offs[i]
    r = col.offs[i+1]
    view(col.vals, l:r-1)
end

@inline function getdata{T}(::Type{T}, col::Column{false,false}, i::Int)
    convert(T, col.vals[i])
end

@inline function getdata{T}(::Type{T}, col::Column{true,false}, i::Int)
    l = col.offs[i]
    r = col.offs[i+1]
    l < r ? convert(T, col.vals[l]) : convert(T, nothing)
end

@inline function getdata{T,OPT}(::Type{T}, col::Column{OPT,true}, i::Int)
    l = col.offs[i]
    r = col.offs[i+1]
    convert(T, view(col.vals, l:r-1))
end

@inline getdata(::Type{Any}, ::Tuple{}, i) = ()

@inline getdata(::Type{Any}, cols::Tuple{Column, Vararg{Any}}, i) =
    (getdata(Any, cols[1], i), getdata(Any, Base.tail(cols), i)...)

@inline getdata(::Type{Tuple{}}, ::Tuple{}, i) = ()

@inline getdata{T<:Tuple{Any,Vararg{Any}}}(::Type{T}, cols::Tuple{Column, Vararg{Any}}, i) =
    (getdata(Base.tuple_type_head(T), cols[1], i), getdata(Base.tuple_type_tail(T), Base.tail(cols), i)...)

@inline getdata{T<:AbstractRecord}(::Type{T}, cols::Tuple{Vararg{Column}}, i) =
    T(getdata(fieldtypes(T), cols, i)...)

@inline getdata{T}(::Type{T}, cols::ColumnVector, i) =
    getdata(T, (cols...), i)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
function poly_eval(p::MP.AbstractPolynomial{JuMP.AffExpr},
                   a::AbstractVector{Float64})
    vars = MP.variables(p)
    aff = zero(JuMP.AffExpr)
    for term in MP.terms(p)
        mono = MP.monomial(term)
        aff = MA.add_mul!!(aff, mono(vars => a), MP.coefficient(term))
    end
    return aff
end
function sublevel_eval(ell::Sets.Ellipsoid,
                       a::AbstractVector)
    return quad_form(ell.Q, a)
end
function sublevel_eval(set::Union{Sets.PolySet,
                                  Sets.ConvexPolySet},
                       a::AbstractVector)
    return poly_eval(MP.polynomial(set.p), a)
end
function sublevel_eval(set::Sets.Householder, a::AbstractVector, β)
    x = Sets.space_variables(set)
    z = Sets.perspective_variable(set)
    # Avoid large values, with high degree polynomials, it might cause issues
    if false
        scale_a = norm(a)
        scale_β = abs(β)
        scaling = max(scale_a, scale_β) / sqrt(scale_a * scale_β)
    else
        scaling = 1.0
    end
    return Sets.perspective_gauge0(set)(z => β / scaling, x => a / scaling)
end


text:
# Float32 -> Float8 algorithm in analogy to
#
# Float32 -> Float16 algorithm from:
#   "Fast Half Float Conversion" by Jeroen van der Zijp
#   ftp://ftp.fox-toolkit.org/pub/fasthalffloatconversion.pdf
#
# With adjustments for round-to-nearest, ties to even.

function create_base_shifttable(::Type{T}) where {T<:AbstractFloat8}

    basetable = Vector{UInt8}(undef, 512)
    shifttable = Vector{UInt8}(undef, 512)

    if T == Float8
        # elements derive from
        # [1]   2^-6 = Float8(0x01) the smallest representable number (subnormal)
        # [2]   2^-2 = Float8(0x10) the first non-subnormal number
        # [3]   2^4 = 16 > floatmax(Float8) is the smallest power of two that is larger than floatmax(Float8)

        e_limits = [-6,-2,4]

        # shift a 0x1 in the exponent bits created by "significand_mask(Float32) + 0x1"
        # to the first significand bit
        # e_shift_subnorm is 17 for Float8
        e_shift_subnorm = n_significant_bits(Float32)-(n_significant_bits(Float8)-1)+e_limits[2]-1
    elseif T == Float8_4

        # see above
        e_limits = [-9,-6,8]

        # shift a 0x1 in the exponent bits created by "significand_mask(Float32) + 0x1"
        # to the first significand bit
        # e_shift_subnorm is 14 for Float8_4
        e_shift_subnorm = n_significant_bits(Float32)-(n_significant_bits(Float8_4)-1)+e_limits[2]-1
    end

    for i = 0:255                               # all possible exponents for Float32
        e = i - 127                             # subtract Float32 bias
        if e < e_limits[1]                      # Very small numbers map to +- zero
            basetable[i|0x000+1] = zero(T)
            basetable[i|0x100+1] = -zero(T)
            shifttable[i|0x000+1] = n_significant_bits(T)+1
            shifttable[i|0x100+1] = n_significant_bits(T)+1
        elseif e < e_limits[2]                  # Small numbers map to denorms
            basetable[i|0x000+1] = zero(T)
            basetable[i|0x100+1] = -zero(T)
            shifttable[i|0x000+1] = -e+e_shift_subnorm
            shifttable[i|0x100+1] = -e+e_shift_subnorm
        elseif e < e_limits[3]                  # Normal numbers just lose precision
            basetable[i|0x000+1] = ((e+bias(T)) << n_significant_bits(T))
            basetable[i|0x100+1] = ((e+bias(T)) << n_significant_bits(T)) | sign_mask(T)
            shifttable[i|0x000+1] = n_significant_bits(Float32)-n_significant_bits(T)
            shifttable[i|0x100+1] = n_significant_bits(Float32)-n_significant_bits(T)
        elseif e < 128                          # Large numbers map to Infinity
            basetable[i|0x000+1] = inf8(T)
            basetable[i|0x100+1] = -inf8(T)
            shifttable[i|0x000+1] = n_significant_bits(T)+1
            shifttable[i|0x100+1] = n_significant_bits(T)+1
        else                                    # Infinity and NaN's stay Infinity and NaN's
            basetable[i|0x000+1] = inf8(T)
            basetable[i|0x100+1] = -inf8(T)
            shifttable[i|0x000+1] = n_significant_bits(Float32)-n_significant_bits(T)
            shifttable[i|0x100+1] = n_significant_bits(Float32)-n_significant_bits(T)
        end
    end

    return basetable, shifttable
end

const basetable8, shifttable8 = create_base_shifttable(Float8)
const basetable8_4, shifttable8_4 = create_base_shifttable(Float8_4)

# function Float8(val::Float32)
#
#     f = reinterpret(UInt32, val)
#
#     if isnan(val)       #TODO retain the significant bits for NaN?
#         return nan8(Float8)
#     end
#
#     # exponent as Int64
#     i = f >> n_significant_bits(Float32) + 1
#     @inbounds sh = shifttable8[i]
#     f &= significand_mask(Float32)
#
#     # If `val` is subnormal, the tables are set up to force the
#     # result to 0, so the significand has an implicit `1` in the
#     # cases we care about.
#
#     f |= significand_mask(Float32) + 0x1
#     @inbounds h = (basetable8[i] + (f >> sh) & significand_mask(Float8)) % UInt8
#
#     # rounding
#     nextbit = (f >> (sh-1)) & 1
#     if nextbit != 0 && (h & exponent_mask(Float8)) != exponent_mask(Float8)
#         # Round halfway to even or check lower bits
#         if h&1 == 1 || (f & ((1<<(sh-1))-1)) != 0
#             h += one(UInt8)
#         end
#     end
#     return reinterpret(Float8, h)
# end

function Float8_4(val::Float32)

    f = reinterpret(UInt32, val)

    if isnan(val)       #TODO retain the significant bits for NaN?
        return nan8(Float8_4)
    end

    # exponent as Int64
    i = f >> n_significant_bits(Float32) + 1
    @inbounds sh = shifttable8_4[i]
    f &= significand_mask(Float32)

    # If `val` is subnormal, the tables are set up to force the
    # result to 0, so the significand has an implicit `1` in the
    # cases we care about.

    f |= significand_mask(Float32) + 0x1
    @inbounds h = (basetable8_4[i] + (f >> sh) & significand_mask(Float8_4)) % UInt8

    # rounding
    nextbit = (f >> (sh-1)) & 1
    if nextbit != 0 && (h & exponent_mask(Float8_4)) != exponent_mask(Float8_4)
        # Round halfway to even or check lower bits
        if h&1 == 1 || (f & ((1<<(sh-1))-1)) != 0
            h += one(UInt8)
        end
    end
    return reinterpret(Float8_4, h)
end


text:
using ExaTron
using MPI
MPI.Init()

comm = MPI.COMM_WORLD

datafile = "./data/"*ARGS[1]
pq_val = parse(Float64, ARGS[2])
va_val = parse(Float64, ARGS[3])
max_iter = parse(Int, ARGS[4])
gpu = parse(Bool, ARGS[5])

# Warm-up
@time ExaTron.admm_rect_gpu_mpi(datafile; iterlim=1, rho_pq=pq_val, rho_va=va_val, use_gpu=gpu, use_polar=true, scale=1e-4)
# Computation
@time ExaTron.admm_rect_gpu_mpi(datafile; iterlim=max_iter, rho_pq=pq_val, rho_va=va_val, use_gpu=gpu, use_polar=true, scale=1e-4)

MPI.Finalize()


text:
#Write a function that computes the running total of a list.

list = [1,2,3,4,5,6,7]

function run_sum!(x)
    sum(x)
end

for i in 1:20
    l = 1:i
    println(l)
    println(sum(l))
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using Statistics

@testsuite "statistics" (AT, eltypes)->begin
    for ET in eltypes
        if !(ET in [Float16, Float32, Float64])
            continue
        end
        @testset "std" begin
            @test compare(std, AT, rand(ET, 10))
            @test compare(std, AT, rand(ET, 10,1,2))
            @test compare(std, AT, rand(ET, 10,1,2); corrected=true)
            @test compare(std, AT, rand(ET, 10,1,2); dims=1)
        end

        @testset "var" begin
            @test compare(var, AT, rand(ET, 10))
            @test compare(var, AT, rand(ET, 10,1,2))
            @test compare(var, AT, rand(ET, 10,1,2); corrected=true)
            @test compare(var, AT, rand(ET, 10,1,2); dims=1)
            @test compare(var, AT, rand(ET, 10,1,2); dims=[1])
            @test compare(var, AT, rand(ET, 10,1,2); dims=(1,))
            @test compare(var, AT, rand(ET, 10,1,2); dims=[2,3])
            @test compare(var, AT, rand(ET, 10,1,2); dims=(2,3))
        end

        @testset "mean" begin
            @test compare(mean, AT, rand(ET, 2, 2))
            @test compare(mean, AT, rand(ET, 2, 2); dims=2)
            @test compare(mean, AT, rand(ET, 2, 2, 2); dims=[1,3])
            @test compare(x->mean(sin, x), AT, rand(ET, 2,2))
            @test compare(x->mean(sin, x; dims=2), AT, rand(ET, 2,2))
            @test compare(x->mean(sin, x; dims=[1,3]), AT, rand(ET, 2,2,2))
        end
    end

    for ET in eltypes
        # Doesn't work with ComplexF32 in oneAPI for some reason.
        if !(ET in [Float32, Float64, Float16, ComplexF64])
            continue
        end
        @testset "cov" begin
            s = 100
            @test compare(cov, AT, rand(ET, s))
            @test compare(cov, AT, rand(ET, s, 2))
            @test compare(cov, AT, rand(ET, s, 2); dims=2)
            if ET <: Real
                @test compare(cov, AT, rand(ET(1):ET(100), s))
            end
        end

        @testset "cor" begin
            s = 100
            @test compare(cor, AT, rand(ET, s))
            @test compare(cor, AT, rand(ET, s, 2))
            @test compare(cor, AT, rand(ET, s, 2); dims=2)
            if ET <: Real
                @test compare(cor, AT, rand(ET(1):ET(100), s))
            end
        end
    end
end