text:
using TopologyPreprocessing
using Test
using Eirene

#%%
@testset "BettiCurves.jl" begin
    sample_distance_matrix1 = [0  1  25 4  5  9  13 17;
                                  1  0  2  26 6  10 14 18;
                                  25 2  0  3  7  11 15 19;
                                  4  26 3  0  8  12 16 20;
                                  5  6  7  8  0  21 27 24;
                                  9  10 11 12 21 0  22 28;
                                  13 14 15 16 27 22 0  23;
                                  17 18 19 20 24 28 23 0 ]
    sample_distance_matrix2 = [1  1  41 4  5  9  13 17 25 33;
                               1  1  2  42 6  10 14 18 26 34;
                               41 2  1  3  7  11 15 19 27 35;
                               4  42 3  1  8  12 16 20 28 36;
                               5  6  7  8  1  21 43 24 29 37;
                               9  10 11 12 21 1  22 44 30 38;
                               13 14 15 16 43 22 1  23 31 39;
                               17 18 19 20 24 44 23 1  32 40;
                               25 26 27 28 29 30 31 32 1  45;
                               33 34 35 36 37 38 39 40 45 1;]
    # get_bettis
    for matrix = [sample_distance_matrix1, sample_distance_matrix2]
        for max_B_dim = 1:4
            eirene_results = eirene(matrix, model="vr", maxdim = max_B_dim)
            all_bettis = get_bettis(eirene_results, max_B_dim)

            @test length(all_bettis) == max_B_dim
            @test all_bettis isa Vector{Array{Float64,2}}
        end

        for max_B_dim = 1:4, min_B_dim = 1:3
            if min_B_dim > max_B_dim
                @debug "Continue at " min_B_dim, max_B_dim
                continue
            end
           eirene_results = eirene(matrix, model="vr", maxdim = max_B_dim)
           all_bettis = get_bettis(eirene_results, max_B_dim, min_dim=min_B_dim)

           @test length(all_bettis) == max_B_dim - (min_B_dim-1)
           @test all_bettis isa Vector{Array{Float64,2}}
        end
    end

    # normalise_bettis
    #  as betticurve results
    for matrix = [sample_distance_matrix1, sample_distance_matrix2]
        for max_B_dim = 1:4, min_B_dim = 1:3
            if min_B_dim > max_B_dim
                @debug "Continue at " min_B_dim, max_B_dim
                continue
            end
            eirene_results = eirene(matrix, model="vr", maxdim = max_B_dim)
            betti_result = betticurve(eirene_results, dim=max_B_dim)

            normed_all_bettis = normalise_bettis(betti_result)
            @test typeof(normed_all_bettis) == typeof(betti_result)
            @test length(normed_all_bettis) != max_B_dim
            @test size(normed_all_bettis) == size(betti_result)
            @test normed_all_bettis isa Array{Float64,2}

            # Betti values are unchanged:
            @test normed_all_bettis[:,2] == betti_result[:,2]
            # Max val is 1
            @test findmax(normed_all_bettis[:,1])[1] == 1.
        end
    end

    #  as get_bettis results
    for matrix = [sample_distance_matrix1, sample_distance_matrix2]
        for max_B_dim = 1:4, min_B_dim = 1:3
            if min_B_dim > max_B_dim
                @debug "Continue at " min_B_dim, max_B_dim
                continue
            end
            eirene_results = eirene(matrix, model="vr", maxdim = max_B_dim)
            all_bettis = get_bettis(eirene_results, max_B_dim)

            normed_all_bettis = normalise_bettis(all_bettis)
            @test typeof(normed_all_bettis) == typeof(all_bettis)
            @test length(normed_all_bettis) == max_B_dim
            @test normed_all_bettis isa Vector{Array{Float64,2}}

            # Betti values are unchanged:
            @test normed_all_bettis[max_B_dim][:,2] == all_bettis[max_B_dim][:,2]
            # Max val is 1
            @test findmax(normed_all_bettis[max_B_dim][:,1])[1] == 1.
        end
    end

    # get_vectorized_bettis
    let max_B_dim = 5,
        min_B_dim = 1,
        eirene_results = eirene(sample_distance_matrix1, model="vr", maxdim = max_B_dim)

        let eirene_bettis = get_bettis(eirene_results, max_B_dim, min_dim=min_B_dim),
            vectorized_bettis = get_vectorized_bettis(eirene_results, max_B_dim, min_dim=min_B_dim)

            @test size(vectorized_bettis)[2] == max_B_dim - (min_B_dim-1)
            for d in min_B_dim:max_B_dim
                @test vectorized_bettis[:,d] == eirene_bettis[d][:,2]
            end
        end
    end
end

@testset "BettiCurves.jl -> plot bettis" begin
    # TODO remove tests which test Plots.plot function and not plot_bettis functionality
    sample_distance_matrix1 = [0  1  25 4  5  9  13 17;
                                  1  0  2  26 6  10 14 18;
                                  25 2  0  3  7  11 15 19;
                                  4  26 3  0  8  12 16 20;
                                  5  6  7  8  0  21 27 24;
                                  9  10 11 12 21 0  22 28;
                                  13 14 15 16 27 22 0  23;
                                  17 18 19 20 24 28 23 0 ]
    sample_distance_matrix2 = [1  1  41 4  5  9  13 17 25 33;
                               1  1  2  42 6  10 14 18 26 34;
                               41 2  1  3  7  11 15 19 27 35;
                               4  42 3  1  8  12 16 20 28 36;
                               5  6  7  8  1  21 43 24 29 37;
                               9  10 11 12 21 1  22 44 30 38;
                               13 14 15 16 43 22 1  23 31 39;
                               17 18 19 20 24 44 23 1  32 40;
                               25 26 27 28 29 30 31 32 1  45;
                               33 34 35 36 37 38 39 40 45 1;]

    # plot_bettis tests for get_bettis:
    let max_B_dim = 5,
        min_B_dim = 1,
        eirene_results = eirene(sample_distance_matrix1, model="vr", maxdim = max_B_dim)
        all_bettis = get_bettis(eirene_results, max_B_dim)

        p = plot_bettis(all_bettis);
        @test length(p.series_list) == max_B_dim-(min_B_dim-1)
        @test p.attr[:plot_title] == ""

        @test_throws DomainError plot_bettis(all_bettis, min_dim=max_B_dim+1)

        for (dim_index, dim)= enumerate(min_B_dim:max_B_dim)
            @test p.series_list[dim_index][:label] == "β$(dim)"
            if !isnan(all_bettis[dim_index][:,1][1])
                @test p.series_list[dim_index][:x] == all_bettis[dim_index][:,1]
                @test p.series_list[dim_index][:y] == all_bettis[dim_index][:,2]
            end
        end

        # for dim = min_B_dim:max_B_dim
        #     p = plot_bettis(all_bettis, min_dim = dim);
        #     @test length(p.series_list) == max_B_dim-(dim-1)
        # end

        let p1 = plot_bettis(all_bettis, betti_labels=false)
            for (dim_index, dim)= enumerate(min_B_dim:max_B_dim)
                @test p1.series_list[dim][:label] == "y$(dim)"
                if !isnan(all_bettis[dim_index][:,1][1])
                    @test p1.series_list[dim][:x] == all_bettis[dim][:,1]
                    @test p1.series_list[dim][:y] == all_bettis[dim][:,2]
                end
            end
        end

        let lw=4,
            p1 = plot_bettis(all_bettis, betti_labels=true, lw=lw)
            for (dim_index, dim)= enumerate(min_B_dim:max_B_dim)
                @test p1.series_list[dim_index][:label] == "β$(dim)"
                @test p1.series_list[dim_index][:linewidth] == lw
            end
        end

        let plt_title = "test_title",
            p1 = plot_bettis(all_bettis, title=plt_title, lw=9, xlabel="2")
            @test_skip p1.attr[:plot_title] == plt_title # why plot-title is not returning the title?
            for dim = min_B_dim:max_B_dim
                @test p1.series_list[dim][:label] == "β$(dim)"
            end
        end

        let plt_title = "test_title",
            lw = 9,
            p1 =  plot_bettis(all_bettis, title=plt_title, lw=lw, xlabel="2", default_labels=false)
            @test_skip p1.attr[:plot_title] == plt_title # why plot-title is not returning the title?
            for dim = min_B_dim:max_B_dim
                @test p1.series_list[dim][:label] == "β$(dim)"
                @test p1.series_list[dim][:linewidth] == lw
                # @test for xlabel
                # @test for no  label
            end
        end

    end
    # plot_bettis tests for get_vectorized_bettis:
    let max_B_dim = 5,
        min_B_dim = 1,
        eirene_results = eirene(sample_distance_matrix1, model="vr", maxdim = max_B_dim)
        all_bettis = get_vectorized_bettis(eirene_results, max_B_dim)

    end
end

@testset "BettiCurves.jl -> area under betti curves" begin

    let sample_distance_matrix1 = [0  1  25 4  5  9  13 17;
                                  1  0  2  26 6  10 14 18;
                                  25 2  0  3  7  11 15 19;
                                  4  26 3  0  8  12 16 20;
                                  5  6  7  8  0  21 27 24;
                                  9  10 11 12 21 0  22 28;
                                  13 14 15 16 27 22 0  23;
                                  17 18 19 20 24 28 23 0 ],
          sample_distance_matrix2 = [1  1  41 4  5  9  13 17 25 33;
                                     1  1  2  42 6  10 14 18 26 34;
                                     41 2  1  3  7  11 15 19 27 35;
                                     4  42 3  1  8  12 16 20 28 36;
                                     5  6  7  8  1  21 43 24 29 37;
                                     9  10 11 12 21 1  22 44 30 38;
                                     13 14 15 16 43 22 1  23 31 39;
                                     17 18 19 20 24 44 23 1  32 40;
                                     25 26 27 28 29 30 31 32 1  45;
                                     33 34 35 36 37 38 39 40 45 1;],
            max_B_dim = 5,
            min_B_dim = 1
    #==
    checks if the size is anyhow changed during proces;
    checks is the values Array{Matrix} and reshaped matrix are the same

    ==#
        for min_B_dim in [1, 2, 3, 4, 5]
            eirene_results1 =
                eirene(sample_distance_matrix1, model = "vr", maxdim = max_B_dim)
            eirene_results2 =
                eirene(sample_distance_matrix2, model = "vr", maxdim = max_B_dim)

            bettis_collection = [
                get_bettis(eirene_results1, max_B_dim),
                get_bettis(eirene_results2, max_B_dim),
            ]

            for bettis_col in bettis_collection
                total_vecs = length(bettis_col)
                vec_len, vec_width = size(bettis_col[1])

                reshaped_betti = TopologyPreprocessing.vectorize_bettis(bettis_col)
                @test vec_len .== size(reshaped_betti, 1)
                @test total_vecs .== size(reshaped_betti, 2)
                for k = 1:total_vecs
                    @test reshaped_betti[:, k] == bettis_col[k][:, 2]
                end
            end
        end
        #==
        checks if get vectorized bettis has same values as get_bettis
        ==#
        for min_B_dim in [1, 2, 3, 4, 5]
            eirene_results1 =
                eirene(sample_distance_matrix1, model = "vr", maxdim = max_B_dim)
            eirene_results2 =
                eirene(sample_distance_matrix2, model = "vr", maxdim = max_B_dim)

            bettis_collection = [
                get_bettis(eirene_results1, max_B_dim),
                get_bettis(eirene_results2, max_B_dim),
            ]
            vec_bett_collection = [
                get_vectorized_bettis(eirene_results1, max_B_dim),
                get_vectorized_bettis(eirene_results2, max_B_dim),
            ]

            for index = 1:length(bettis_collection)
                bettis_col = bettis_collection[index]
                vec_bettis_col = vec_bett_collection[index]

                total_vecs = length(bettis_col)
                for k = 1:total_vecs
                    @test vec_bettis_col[:, k] == bettis_col[k][:, 2]
                end
            end
        end
    end
end


text:
@testset "Connlu utils" begin

  filename = "test/fixtures/connlu_test.txt"
  correct_sentences = [
    "Attribution-ShareAlike 2.0.",
    "Any use of the work other than as authorized under this license or copyright law is prohibited.",
    "UNLESS OTHERWISE AGREED TO BY THE PARTIES IN WRITING, LICENSOR OFFERS THE WORK AS-IS AND MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE MATERIALS, EXPRESS, IMPLIED, STATUTORY OR OTHERWISE, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF TITLE, MERCHANTIBILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, ACCURACY, OR THE PRESENCE OF ABSENCE OF ERRORS, WHETHER OR not DISCOVERABLE.",
  ]
  load_result = DependencyParser.DependencyParsing.load_connlu_file(filename)

  foreach(correct_sentences) do correct_sentence
    @test correct_sentence in map(connlu_s -> connlu_s.string_doc.text, load_result)
  end
end

ext: jl

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using ModelingToolkit, SparseArrays, Test, GalacticOptim, GalacticOptimJL

@variables x y
@parameters a b
loss = (a - x)^2 + b * (y - x^2)^2
sys1 = OptimizationSystem(loss, [x, y], [a, b], name = :sys1)
sys2 = OptimizationSystem(loss, [x, y], [a, b], name = :sys2)

@variables z
@parameters β
loss2 = sys1.x - sys2.y + z * β
combinedsys = OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2],
                                 name = :combinedsys)

equations(combinedsys)
states(combinedsys)
parameters(combinedsys)

calculate_gradient(combinedsys)
calculate_hessian(combinedsys)
generate_function(combinedsys)
generate_gradient(combinedsys)
generate_hessian(combinedsys)
ModelingToolkit.hessian_sparsity(combinedsys)

u0 = [sys1.x => 1.0
      sys1.y => 2.0
      sys2.x => 3.0
      sys2.y => 4.0
      z => 5.0]
p = [sys1.a => 6.0
     sys1.b => 7.0
     sys2.a => 8.0
     sys2.b => 9.0
     β => 10.0]

prob = OptimizationProblem(combinedsys, u0, p, grad = true)
sol = solve(prob, NelderMead())
@test sol.minimum < -1e5

prob2 = remake(prob, u0 = sol.minimizer)
sol = solve(prob, BFGS(initial_stepnorm = 0.0001), allow_f_increases = true)
@test sol.minimum < -1e8
sol = solve(prob2, BFGS(initial_stepnorm = 0.0001), allow_f_increases = true)
@test sol.minimum < -1e9

rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p = [1.0, 100.0]

f = OptimizationFunction(rosenbrock, GalacticOptim.AutoModelingToolkit())
prob = OptimizationProblem(f, x0, _p)
sol = solve(prob, Newton())

# issue #819
@testset "Combined system name collisions" begin
    sys2 = OptimizationSystem(loss, [x, y], [a, b], name = :sys1)
    @test_throws ArgumentError OptimizationSystem(loss2, [z], [β], systems = [sys1, sys2])
end

# observed variable handling
@variables OBS
@named sys2 = OptimizationSystem(loss, [x, y], [a, b]; observed = [OBS ~ x + y])
OBS2 = OBS
@test isequal(OBS2, @nonamespace sys2.OBS)
@unpack OBS = sys2
@test isequal(OBS2, OBS)


text:
module StokesNitsche

# Testing Stokes problem with Dirichlet BCs for the velocity
# and pressure prescribed in a single point

using Test
using Gridap
import Gridap: ∇

#
const T = VectorValue{2,Float64}

u(x) = VectorValue(x[1]*x[1], x[2])
∇u(x) = TensorValue(2*x[1],0.0,0.0,1.0)
Δu(x) = VectorValue(2.0,0.0)

p(x) = x[1] - x[2]
∇p(x) = VectorValue(1.0,-1.0)

∇(::typeof(u)) = ∇u

n = 2
order = 2

# function run_test(n,order)

f(x) = - Δu(x) + ∇p(x)
g(x) = tr(∇u(x)) # i.e. div(u(x))
ud(x) = u(x)

# Discrete model
L = 1.0
limits = (0.0, L, 0.0, L)
ncellx = n
model = CartesianDiscreteModel(limits, (ncellx,ncellx))

h = L / ncellx

γ = order*(order+1)
γ0 = 1.0/10.0

# Construct the FEspace 1
# labels = FaceLabels(model)
V = TestFESpace(
reffe=:Lagrangian,
conformity=:H1,
valuetype=VectorValue{2,Float64},
model=model,
order=order)

# Construct the FEspace 2
Q = TestFESpace(
reffe=:PLagrangian,
conformity=:L2,
valuetype=Float64,
model=model,
order=order-1,
constraint=:zeromean)

# Define test and trial
Y = MultiFieldFESpace([V,Q])
# Y = [V,Q]

U = TrialFESpace(V)
P = TrialFESpace(Q)

X = MultiFieldFESpace([U,P])
# X = [U,P]

# Define integration mesh and quadrature for volume
trian = get_triangulation(model)
quad = CellQuadrature(trian,2*order)

btrian = BoundaryTriangulation(model)
bquad = CellQuadrature(btrian,2*order)
nb = get_normal_vector(btrian)

# Weak form
function A_Ω(x,y)
  u, p = x
  v, q = y
  ∇(v)⊙∇(u) - q*(∇⋅u) - (∇⋅v)*p
end

function B_Ω(y)
  v, q = y
  v⋅f - q*g
end

function A_∂Ω(x,y)
  u, p = x
  v, q = y
  (γ/h)*v⋅u - v⋅(nb⋅∇(u)) - (nb⋅∇(v))⋅u + (p*nb)⋅v + (q*nb)⋅u
end

function B_∂Ω(y)
  v, q = y
  (γ/h)*v⋅u - (nb⋅∇(v))⋅u + (q*nb)⋅u
end

t_Ω = AffineFETerm(A_Ω,B_Ω,trian,quad)

# t_∂Ω = FESource(B_∂Ω,btrian,bquad)
t_∂Ω = AffineFETerm(A_∂Ω,B_∂Ω,btrian,bquad)

# Dummy term with 0 facets to check degenerated case
trian0 = BoundaryTriangulation(model,fill(false,Gridap.ReferenceFEs.num_facets(model)))
quad0 = CellQuadrature(trian0,2*order)
s0(x) = VectorValue(2.0*x[1],2.0)
function L0(y)
  v,q = y
  s0⋅v
end
t_0 = FESource(L0,trian0,quad0)

# Define the FEOperator
op = AffineFEOperator(X,Y,t_Ω,t_∂Ω,t_0)
# op = LinearFEOperator(Yh,Xh,t_Ω)

# Solve!
xh = solve(op)
uh, ph = xh

# Define exact solution and error
eu = u - uh

ep = p - ph

# writevtk(trian,"trian",cellfields=["uh"=>uh,"ph"=>ph, "eu"=>eu, "ep"=>ep])

# Define norms to measure the error
l2(u) = u⊙u
h1(u) = ∇(u)⊙∇(u) + l2(u)

# Compute errors
eul2 = sqrt(sum( integrate(l2(eu),trian,quad) ))
euh1 = sqrt(sum( integrate(h1(eu),trian,quad) ))

epl2 = sqrt(sum( integrate(l2(ep),trian,quad) ))

@test eul2 < 1.e-8
@test euh1 < 1.e-8

@test epl2 < 1.e-8

end


text:
module AiyagariContinuousTime
import Distances: chebyshev
using Gensys
using Calculus

##############################################################################
##
## Load files
##
##############################################################################
include("aiyagari.jl")
include("dynamicaiyagari.jl")
##############################################################################
##
## Exported methods and types
##
##############################################################################
export solve,
AiyagariProblem,
AiyagariArrays,
AiyagariSolution,
DynamicAiyagariSolution,
simulate
end


text:
module AutoregressiveNodeTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules

@testset "AutoregressiveNode" begin

    @testset "ARTransitionMatrix" begin

        rng = MersenneTwister(1233)

        for γ in randn(rng, 3), order in 2:4
            transition = ReactiveMP.ARTransitionMatrix(order, γ)
            matrix      = rand(rng, order, order)
            ftransition = zeros(order, order)
            ftransition[1] = inv(γ)

            @test broadcast(+, matrix, transition) == (matrix + ftransition)
            @test_throws DimensionMismatch broadcast(+, zeros(order + 1, order + 1), transition)

            @test ReactiveMP.add_transition(matrix, transition) == (matrix + ftransition)
            @test_throws DimensionMismatch ReactiveMP.add_transition(zeros(order + 1, order + 1), transition)

            cmatrix = copy(matrix)
            broadcast!(+, cmatrix, transition)
            @test cmatrix == (matrix + ftransition)

            cmatrix = copy(matrix)
            ReactiveMP.add_transition!(cmatrix, transition)
            @test cmatrix == (matrix + ftransition)
        end

    end
end

end


text:
const stdlibs = isdefined(Pkg.Types, :stdlib) ? Pkg.Types.stdlib : Pkg.Types.stdlibs
const BUILTIN_PKGS = Dict(v=>string(k) for (k, v) in stdlibs())


text:
export G

"""

    G(p::Integer, q::Integer, r::Integer)

A convenient method to instantiate a Geometric Algebra with `p` positive, `q` negative and `r` zero dimensions.

`q` and `r` defaults to zero if omitted.

The basis will be named with prefix `e` and indexed starting from 1.

For example:

```julia
# Basic
Hyper   = G(1)       # Hyperbolic numbers.
ℂ       = G(0,1)     # Complex numbers.
Dual    = G(0,0,1)   # Dual numbers.
ℍ       = G(0,2)     # Quaternions.

# Clifford
Cl2 =       G(2)     # Clifford algebra for 2D vector space.
Cl3 =       G(3)     # Clifford algebra for 3D vector space.
Spacetime = G(1,3)   # Clifford algebra for timespace vectors.

# Geometric
PGA2D = G(2,0,1)     # Projective Euclidean 2D plane. (dual)
PGA3D = G(3,0,1)     # Projective Euclidean 3D space. (dual)
CGA2D = G(3,1)       # conformal 2D space.
CGA3D = G(4,1)       # Conformal 3D space.

# High-Dimensional GA
DCGA3D = G(6,2)      # Double Conformal 3D Space.
TCGA3D = G(9,3)      # Tripple Conformal 3D Space.
DCGSTA = G(4,8)      # Double Conformal Geometric Space Time Algebra.
QCGA   = G(9,6)      # Quadric Conformal Geometric Algebra.
```

To instantiate a Geometric Algebra with more parameters, use `galgebra.ga.Ga` instead.

For example:

```julia
import SymPy: sympy
using GAlgebra

Ga = galgebra.ga.Ga

g3d = Ga("e*x|y|z")

(r, th, phi) = coords = sympy.symbols("r theta phi")
s3d = Ga("e_r e_theta e_phi", g=[1 0 0; 0 r^2 0; 0 0 r^2 * sympy.sin(th)^2], coords=coords, norm=true)
(er, eth, ephi) = s3d.mv()
```

Please also consult the documentation of [GAlgebra](https://github.com/pygae/galgebra).

"""
function G(p::Integer, q::Integer, r::Integer)
    total = p + q + r
    basis = "e*" * join([string(i) for i in 1:total], "|")
    metric = [fill(1, p) ; fill(-1, q) ; fill(0, r)]
    galgebra.ga.Ga(basis, g=metric)
end

G(p::Integer, q::Integer) = G(p, q, 0)

G(p::Integer) = G(p, 0, 0)


text:
export factors
export primefactors
export totient

"""
Return sorted list of unique factors for a number n.
"""
function factors(n::Integer)::Array{Integer}
    factors = Set{Integer}(1)
    for i in 2:Integer(floor(sqrt(n)))
        if n % i == 0
            push!(factors, i)
            push!(factors, div(n, i))
        end
    end
    return sort(collect(factors))
end

"""
Return sorted list of prime factors for number n.
"""
function primefactors(n::Integer)::Array{Integer}
    factors = ones(Integer, 0)

    # Store all of the factors of 2.
    while n % 2 == 0
        push!(factors, 2)
        n = fld(n, 2)
    end

    # Store all other prime factors, increasing to sqrt(n).
    for i in 3:2:Integer(floor(sqrt(n)))
        while n % i == 0
            push!(factors, i)
            n = fld(n, i)
        end
    end

    # Store last prime factor.
    if n > 1
        push!(factors, n)
    end

    return sort(factors)
end

"""
Euler's totient function returns the number of positive integers up to n that
are releatively prime to n.  Pull out the prime factors f and take a produce of
n and all (1 - 1/f).
"""
function totient(n::Integer)::Integer
    factors = primefactors(n)
    last_prime = 0
    result = Float64(n)
    for f in factors
        if f == last_prime
            continue
        else
            result *= (1 - 1 / f)
            last_prime = f
        end
    end
    return Integer(round(result))
end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# Autogenerated wrapper script for Qt5QuickControls2_jll for x86_64-w64-mingw32-cxx03
export libqt5quickcontrols2, libqt5quicktemplates2, quickcontrols2

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5QuickControls2")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "Qt5QuickControls2.dll")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "Qt5QuickTemplates2.dll")
JLLWrappers.@declare_file_product(quickcontrols2)
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "bin\\Qt5QuickControls2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "bin\\Qt5QuickTemplates2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        quickcontrols2,
        "qml\\QtQuick\\Controls.2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, Web3

connection = Web3Connection("https://localhost:8545")
abi = """[{"constant":false,"inputs":[{"name":"i","type":"int32"}],"name":"add","outputs":[],"payable":false,"stateMutability":"nonpayable","type":"function"}]"""

Web3.hash(hashes::Dict, str::String) = hashes[str]

hashes = Dict([
    "add(int32)" => hex2bytes("57b8a50f")
])

contract = "0x694b0b0853c64aa51ebce186dd27bd676486cb43"

txn = Dict{String,Any}(
    "gasPrice" => "0x4a817c800",
    "r" => "0x9a8d626b553ac295eb6c0cf692e5d88b97c073d40e1093122881073f18180fda",
    "blockNumber" => "0x7",
    "value" => "0x0",
    "gas" => "0x6691b7",
    "s" => "0x6acc9e3f3f04edf65a25fd63681af90fe197c4a16978c2646a3d368573fdf69f",
    "v" => "0x1c",
    "hash" => "0x38f537910e76108cb65807939ca3ed72858e2f8ee4520b58e1c6069a70e9384a",
    "transactionIndex" => "0x0",
    "input" => "0x57b8a50f0000000000000000000000000000000000000000000000000000000000000003",
    "blockHash" => "0x09f1aebdc61f76fc906fd95176f72960f04c72c5da28efd620df60b2169a2829",
    "to" => "0x694b0b0853c64aa51ebce186dd27bd676486cb43",
    "from" => "0xdb810d2557b2e9be3d86ec02e4ea6b4aaa93b7c7",
    "nonce" => "0x6"
)

# install contract
@test readABI(connection, contract, IOBuffer(abi)) != nothing

# test basic function call
decoded = decodefunctioncall(IOBuffer(hex2bytes(txn["input"][3:end])), contracts[txn["to"]])
data = encodefunctioncall(IOBuffer(UInt8[], write=true), decoded.decl, decoded.inputs)
@test txn["input"] == "0x" * bytes2hex(data)


text:
"""
    1D1F1V: flux_kcu!(fw, ff, wL, fL, wR, fR, u, ω, inK, γ, visRef, visIdx, Pr, dt)
    1D2F1V: flux_kcu!(fw, fh, fb, wL, hL, bL, wR, hR, bR, u, ω, inK, γ, visRef, visIdx, Pr, dt)
    1D4F1V: flux_kcu!(fw, fh0, fh1, fh2, fh3, wL, h0L, h1L, h2L, h3L, wR, h0R, h1R, h2R, h3R, u, ω, inK, γ, visRef, visIdx, Pr, dt)
    2D1F2V: flux_kcu!(fw, ff, wL, fL, wR, fR, u, v, ω, inK, γ, visRef, visIdx, Pr, dt, len)
    2D2F2V: flux_kcu!(fw, fh, fb, wL, hL, bL, wR, hR, bR, u, v, ω, inK, γ, visRef, visIdx, Pr, dt, len)
    2D3F2V: flux_kcu!(fw, fh0, fh1, fh2, wL, h0L, h1L, h2L, wR, h0R, h1R, h2R, u, v, ω, inK, γ, visRef, visIdx, Pr, dt, len)

Kinetic central-upwind (KCU) method

- @args: particle distribution functions and their slopes at left/right sides of interface
- @args: particle velocity quadrature points and weights
- @args: time step and cell size

"""
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,1},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,1},
    B<:AbstractArray{<:AbstractFloat,1},
} # 1D1F1V

    # upwind reconstruction
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mxi1, 0, 0)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mxi2, 0, 0)

    w = similar(wL)
    @. w = primL[1] * Muv1 + primR[1] * Muv2

    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g = maxwellian(u, prim)

    fw[1] += Mt[2] * sum(ω .* u .* f)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* f)
    fw[3] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* f))

    @. ff = Mt[1] * u * g + Mt[2] * u * f

    return nothing

end

#--- mixture ---#
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:AbstractFloat,2},
    Z<:AbstractArray{<:Real,2},
    A<:AbstractArray{<:AbstractFloat,2},
    B<:AbstractArray{<:AbstractFloat,2},
}

    # upwind reconstruction
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mxi1, 0, 0)
    Muv2 = mixture_moments_conserve(MuR2, Mxi2, 0, 0)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)

    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    #@. tau +=
    #    abs(cellL.prim[1, :] / cellL.prim[end, :] - cellR.prim[1, :] / cellR.prim[end, :]) /
    #    (cellL.prim[1, :] / cellL.prim[end, :] + cellR.prim[1, :] / cellR.prim[end, :]) *
    #    dt *
    #    2.0
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mxi, 1, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    M = mixture_maxwellian(u, prim)

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* f[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* f[:, j])
        fw[3, j] += Mt[2, j] * 0.5 * sum(ω[:, j] .* u[:, j] .^ 3 .* f[:, j])

        @. ff[:, j] = Mt[1, j] * u[:, j] * M[:, j] + Mt[2, j] * u[:, j] * f[:, j]
    end

    return nothing

end

# ------------------------------------------------------------
# 1D1F3V
# ------------------------------------------------------------
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    uVelo::B,
    vVelo::B,
    wVelo::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,3},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,3},
    B<:AbstractArray{<:AbstractFloat,3},
} # 1D1F1V

    # upwind reconstruction
    δ = heaviside.(uVelo)
    f = @. fL * δ + fR * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    w = similar(wL)
    @. w = primL[1] * Muv1 + primR[1] * Muv2

    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g = maxwellian(uVelo, vVelo, wVelo, prim)

    fw[1] += Mt[2] * sum(ω .* uVelo .* f)
    fw[2] += Mt[2] * sum(ω .* uVelo .^ 2 .* f)
    fw[3] += Mt[2] * sum(ω .* uVelo .* vVelo .* f)
    fw[4] += Mt[2] * sum(ω .* uVelo .* wVelo .* f)
    fw[5] += Mt[2] * 0.5 * sum(ω .* uVelo .* (uVelo .^ 2 .+ vVelo .^ 2 .+ wVelo .^ 2) .* f)

    @. ff = Mt[1] * uVelo * g + Mt[2] * uVelo * f

    return nothing

end

# ------------------------------------------------------------
# 1D2F1V
# ------------------------------------------------------------
function flux_kcu!(
    fw::X,
    fh::Y,
    fb::Y,
    wL::Z,
    hL::A,
    bL::A,
    wR::Z,
    hR::A,
    bR::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,1},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,1},
    B<:AbstractArray{<:AbstractFloat,1},
}

    # upwind reconstruction
    δ = heaviside.(u)
    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mxi1, 0, 0)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mxi2, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)

    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    Mh = maxwellian(u, prim)
    Mb = Mh .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[2] * sum(ω .* u .* h)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h)
    fw[3] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* h) + sum(ω .* u .* b))

    @. fh = Mt[1] * u * Mh + Mt[2] * u * h
    @. fb = Mt[1] * u * Mb + Mt[2] * u * b

    return nothing

end

#--- mixture ---#
function flux_kcu!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    ω::T5,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
) where {
    T1<:AbstractArray{<:AbstractFloat,2},
    T2<:AbstractArray{<:AbstractFloat,2},
    T3<:AbstractArray{<:Real,2},
    T4<:AbstractArray{<:AbstractFloat,2},
    T5<:AbstractArray{<:AbstractFloat,2},
}

    # upwind reconstruction
    δ = heaviside.(u)

    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mxi1, 0, 0)
    Muv2 = mixture_moments_conserve(MuR2, Mxi2, 0, 0)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)

    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    #@. tau +=
    #    abs(cellL.prim[1, :] / cellL.prim[end, :] - cellR.prim[1, :] / cellR.prim[end, :]) /
    #    (cellL.prim[1, :] / cellL.prim[end, :] + cellR.prim[1, :] / cellR.prim[end, :]) *
    #    dt *
    #    2.0
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mxi, 1, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    MH = mixture_maxwellian(u, prim)
    MB = similar(MH)
    for j in axes(MB, 2)
        MB[:, j] .= MH[:, j] .* inK ./ (2.0 * prim[end, j])
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* h[:, j])
        fw[3, j] +=
            Mt[2, j] *
            0.5 *
            (sum(ω[:, j] .* u[:, j] .^ 3 .* h[:, j]) + sum(ω[:, j] .* u[:, j] .* b[:, j]))

        @. fh[:, j] = Mt[1, j] * u[:, j] * MH[:, j] + Mt[2, j] * u[:, j] * h[:, j]
        @. fb[:, j] = Mt[1, j] * u[:, j] * MB[:, j] + Mt[2, j] * u[:, j] * b[:, j]
    end

    return nothing

end

# ------------------------------------------------------------
# 1D4F1V
# ------------------------------------------------------------
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    h3L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    h3R::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,1},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,1},
    B<:AbstractArray{<:AbstractFloat,1},
}

    # upwind reconstruction
    δ = heaviside.(u)

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)
    h3 = @. h3L * δ + h3R * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)

    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g0 = maxwellian(u, prim)
    g1 = Mv[1] .* g0
    g2 = Mw[1] .* g0
    g3 = (Mv[2] + Mw[2]) .* g0

    fw[1] += Mt[2] * sum(ω .* u .* h0)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h0)
    fw[3] += Mt[2] * sum(ω .* u .* h1)
    fw[4] += Mt[2] * sum(ω .* u .* h2)
    fw[5] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* h0) + sum(ω .* u .* h3))

    @. fh0 = Mt[1] * u * g0 + Mt[2] * u * h0
    @. fh1 = Mt[1] * u * g1 + Mt[2] * u * h1
    @. fh2 = Mt[1] * u * g2 + Mt[2] * u * h2
    @. fh3 = Mt[1] * u * g3 + Mt[2] * u * h3

    return nothing

end

#--- mixture ---#
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    h3L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    h3R::A,
    u::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    isMHD = false::Bool,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:AbstractFloat,2},
    Z<:AbstractArray{<:Real,2},
    A<:AbstractArray{<:AbstractFloat,2},
    B<:AbstractArray{<:AbstractFloat,2},
}

    #--- upwind reconstruction ---#
    δ = heaviside.(u)

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)
    h3 = @. h3L * δ + h3R * (1.0 - δ)

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv2 = mixture_moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)

    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    #@. tau +=
    #    abs(primL[1, :] / primL[end, :] - primR[1, :] / primR[end, :]) /
    #    (primL[1, :] / primL[end, :] + primR[1, :] / primR[end, :]) *
    #    dt *
    #    5.0
    if !isMHD
        prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)
    end

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mv, Mw, 1, 0, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    g0 = mixture_maxwellian(u, prim)

    g1 = similar(h0)
    g2 = similar(h0)
    g3 = similar(h0)
    for j in axes(g0, 2)
        g1[:, j] .= Mv[1, j] .* g0[:, j]
        g2[:, j] .= Mw[1, j] .* g0[:, j]
        g3[:, j] .= (Mv[2, j] + Mw[2, j]) .* g0[:, j]
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h0[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* h0[:, j])
        fw[3, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h1[:, j])
        fw[4, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h2[:, j])
        fw[5, j] +=
            Mt[2, j] *
            0.5 *
            (
                sum(ω[:, j] .* u[:, j] .^ 3 .* h0[:, j]) +
                sum(ω[:, j] .* u[:, j] .* h3[:, j])
            )

        @. fh0[:, j] = Mt[1, j] * u[:, j] * g0[:, j] + Mt[2, j] * u[:, j] * h0[:, j]
        @. fh1[:, j] = Mt[1, j] * u[:, j] * g1[:, j] + Mt[2, j] * u[:, j] * h1[:, j]
        @. fh2[:, j] = Mt[1, j] * u[:, j] * g2[:, j] + Mt[2, j] * u[:, j] * h2[:, j]
        @. fh3[:, j] = Mt[1, j] * u[:, j] * g3[:, j] + Mt[2, j] * u[:, j] * h3[:, j]
    end

    return nothing

end

# ------------------------------------------------------------
# 2D1F2V
# ------------------------------------------------------------
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,2},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,2},
    B<:AbstractArray{<:AbstractFloat,2},
}

    # reconstruct initial distribution
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv * len

    ## flux from f0
    g = maxwellian(u, v, prim)

    fw[1] += Mt[2] * sum(ω .* u .* f) * len
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* f) * len
    fw[3] += Mt[2] * sum(ω .* v .* u .* f) * len
    fw[4] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* f)) * len

    @. ff = (Mt[1] * u * g + Mt[2] * u * f) * len

    return nothing

end

# ------------------------------------------------------------
# 2D2F2V
# ------------------------------------------------------------
function flux_kcu!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {
    T1<:AbstractArray{<:AbstractFloat,1},
    T2<:AbstractArray{<:AbstractFloat,2},
    T3<:AbstractArray{<:Real,1},
    T4<:AbstractArray{<:AbstractFloat,2},
    T5<:AbstractArray{<:AbstractFloat,2},
}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv * len

    ## flux from f0
    H = maxwellian(u, v, prim)
    B = H .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[2] * sum(ω .* u .* h) * len
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h) * len
    fw[3] += Mt[2] * sum(ω .* v .* u .* h) * len
    fw[4] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h) + sum(ω .* u .* b)) * len

    @. fh = (Mt[1] * u * H + Mt[2] * u * h) * len
    @. fb = (Mt[1] * u * B + Mt[2] * u * b) * len

    return nothing

end

# ------------------------------------------------------------
# 2D3F2V
# ------------------------------------------------------------
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:AbstractFloat,2},
    Z<:AbstractArray{<:Real,1},
    A<:AbstractArray{<:AbstractFloat,2},
    B<:AbstractArray{<:AbstractFloat,2},
}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    H0 = maxwellian(u, v, prim)
    H1 = H0 .* prim[4]
    H2 = H0 .* (prim[4]^2 + 1.0 / (2.0 * prim[5]))

    fw[1] += Mt[2] * sum(ω .* u .* h0)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h0)
    fw[3] += Mt[2] * sum(ω .* v .* u .* h0)
    fw[4] += Mt[2] * sum(ω .* u .* h1)
    fw[5] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h0) + sum(ω .* u .* h2))

    @. fw *= len
    @. fh0 = (Mt[1] * u * H0 + Mt[2] * u * h0) * len
    @. fh1 = (Mt[1] * u * H1 + Mt[2] * u * h1) * len
    @. fh2 = (Mt[1] * u * H2 + Mt[2] * u * h2) * len

    return nothing

end

#--- mixture ---#
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    len,
    isMHD = false::Bool,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:AbstractFloat,3},
    Z<:AbstractArray{<:Real,2},
    A<:AbstractArray{<:AbstractFloat,3},
    B<:AbstractArray{<:AbstractFloat,3},
}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv2 = mixture_moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)
    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    #@. tau +=
    #    abs(primL[1, :] / primL[end, :] - primR[1, :] / primR[end, :]) /
    #    (primL[1, :] / primL[end, :] + primR[1, :] / primR[end, :]) *
    #    dt *
    #    2.0
    if !isMHD
        prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)
    end

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)

    ## flux from M0
    Muv = mixture_moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    H0 = mixture_maxwellian(u, v, prim)
    H1 = similar(H0)
    H2 = similar(H0)
    for j in axes(H0, 3)
        H1[:, :, j] = H0[:, :, j] .* prim[4, j]
        H2[:, :, j] .= H0[:, :, j] .* (prim[4, j]^2 + 1.0 / (2.0 * prim[5, j]))
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .* h0[:, :, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* h0[:, :, j])
        fw[3, j] += Mt[2, j] * sum(ω[:, :, j] .* v[:, :, j] .* u[:, :, j] .* h0[:, :, j])
        fw[4, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .* h1[:, :, j])
        fw[5, j] +=
            Mt[2, j] *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    h0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .* h2[:, :, j])
            )

        @. fh0[:, :, j] =
            Mt[1, j] * u[:, :, j] * H0[:, :, j] + Mt[2, j] * u[:, :, j] * h0[:, :, j]
        @. fh1[:, :, j] =
            Mt[1, j] * u[:, :, j] * H1[:, :, j] + Mt[2, j] * u[:, :, j] * h1[:, :, j]
        @. fh2[:, :, j] =
            Mt[1, j] * u[:, :, j] * H2[:, :, j] + Mt[2, j] * u[:, :, j] * h2[:, :, j]
    end

    @. fw *= len
    @. fh0 *= len
    @. fh1 *= len
    @. fh2 *= len

    return nothing

end


text:
@testset "Inference" begin
f = readc3d(joinpath(datadir, "sample01", "Eb015pr.c3d"))

@test @inferred(f.groups[:ANALOG][Vector{Int16}, :OFFSET]) == fill(Int16(2048), 32)
@test @inferred(f.groups[:ANALOG][Int16, :USED]) == 16

# Test implicit conversion
@test @inferred(f.groups[:ANALOG][Int, :USED]) == 16
@test @inferred(f.groups[:ANALOG][Vector{Int}, :OFFSET]) == fill(2048, 32)
end


text:
function flux(u, q, w, v, x, t, mu, eta)
    z = x[1];
    r = x[2];
    f = mu[1]*r*q;
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    z = x[1];
    r = x[2];
    s = sin(r)/exp(z);
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    z = x[1];
    r = x[2];
    uexact = exp(-z)*cos(r);
    ub = [u uexact uexact uexact];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    tm = f[1]*n[1] + f[2]*n[2] + tau[1]*(u[1]-uhat[1]);
    fb = [0.0 tm tm tm];
    return fb;
end
function initu(x, mu, eta)
    u0 = 0.0;
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function pymatgen_lattice(lattice::AbstractArray)
    # Notice: Pymatgen uses rows as lattice vectors, so we unpeel
    # our lattice column by column. The default unit in pymatgen is Ångström
    Lattice = pyimport("pymatgen.core.lattice").Lattice
    bohr_to_A = 1 / austrip(1u"Å")
    Lattice([Array(bohr_to_A .* lattice[:, i]) for i in 1:3])
end
pymatgen_lattice(model::Model) = pymatgen_lattice(model.lattice)

function pymatgen_structure(model_or_lattice, atoms, positions)
    Structure = pyimport("pymatgen.core.structure").Structure
    Structure(pymatgen_lattice(model_or_lattice),
              charge_nuclear.(atoms),
              Vector{Float64}.(positions))
end
pymatgen_structure(model::Model) = pymatgen_structure(model, model.atoms, model.positions)

function load_lattice_pymatgen(pyobj::PyObject)
    Structure = pyimport("pymatgen.core.structure").Structure
    Lattice   = pyimport("pymatgen.core.lattice").Lattice

    if pyisinstance(pyobj, Structure)
        load_lattice_pymatgen(pyobj.lattice)
    elseif pyisinstance(pyobj, Lattice)
        lattice = zeros(3, 3)
        for i in 1:3, j in 1:3
            lattice[i, j] = austrip(get(get(pyobj.matrix, j-1), i-1) * u"Å")
        end
        Mat3(lattice)
    else
        error("load_lattice_pymatgen not implemented for python type $pyobj")
    end
end

function load_atoms_pymatgen(pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("pymatgen.core.structure").Structure)
    [ElementCoulomb(spec.number) for spec in pyobj.species]
end

function load_positions_pymatgen(pyobj::PyObject)
    @assert pyisinstance(pyobj, pyimport("pymatgen.core.structure").Structure)
    [Vec3{Float64}(site.frac_coords) for site in pyobj.sites]
end


text:
using LinearAlgebra

"""
Performs a random hill climb solution for a discretized search space.

    Args:
      fn (Function): Objective function
      x0 (Array{Any}): Starting point
      h (Float64, optional): Discretization
      max_failed (Int64, optional): Max number of failed moves before exiting.
        Defaults to the dimensionality of the design space.
      verbose (Bool, optional): Optionally print results

    Returns:
      tuple: (x, fx, f_calls) where x is the minimizer, fx is the minimum, and
      f_calls is the number of function calls made.
"""
function random_hill(fn::Function, x0::Array{Float64};
                     h=.1, max_failed=nothing, verbose=false)
    x = copy(x0)
    fx = fn(x)
    n = length(x)
    failed = 0
    f_calls = 1

    # use an I matrix to randomly determine direction
    e_hat = Matrix(1I, n, n)

    if isnothing(max_failed)
        max_failed = n
    end

    while failed < max_failed
        d = e_hat[:, rand(1:n)] # randomly choose basis vec
        d *= rand([-1, 1]) # randomly choose pos/neg

        xnew = x + h*d
        fnew = fn(xnew)
        f_calls += 1

        if fnew < fx
            fx = fnew
            x = xnew
            failed = 0
        else
            failed += 1
        end
    end

    if verbose
        println("Iterations: $(f_calls)\nx*: $(x)\nf*: $(fx)\n")
    end

    return x, fx, f_calls
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "CSPICE"
version = v"66.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource(
        "http://naif.jpl.nasa.gov/pub/naif/toolkit//C/PC_Linux_GCC_64bit/packages/cspice.tar.Z",
        "93cd4fbce5818f8b7fecf3914c5756b8d41fd5bdaaeac1f4037b5a5410bc4768",
    ),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd cspice
cp -r $WORKSPACE/srcdir/cmake .
mv cmake/CMakeLists.txt .
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/dskx02.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/subpnt.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/inquire.patch"
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/license/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcspice", :libcspice),
    LibraryProduct("libcsupport", :libcsupport),
    ExecutableProduct("brief", :brief),
    ExecutableProduct("chronos", :chronos),
    ExecutableProduct("ckbrief", :ckbrief),
    ExecutableProduct("commnt", :commnt),
    ExecutableProduct("dskbrief", :dskbrief),
    ExecutableProduct("dskexp", :dskexp),
    ExecutableProduct("frmdiff", :frmdiff),
    ExecutableProduct("inspekt", :inspekt),
    ExecutableProduct("mkdsk", :mkdsk),
    ExecutableProduct("mkspk", :mkspk),
    ExecutableProduct("msopck", :msopck),
    ExecutableProduct("spacit", :spacit),
    ExecutableProduct("spkdiff", :spkdiff),
    ExecutableProduct("spkmerge", :spkmerge),
    ExecutableProduct("tobin", :tobin),
    ExecutableProduct("toxfr", :toxfr),
    ExecutableProduct("version", :version),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function _pthread_macro(pool, log, args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in thread macro"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to thread macro"))
    end
    if ex.head === :for
        if ex.args[1] isa Expr && ex.args[1].head === :(=)
            index = ex.args[1].args[1]
            range = ex.args[1].args[2]
            body = ex.args[2]
            return quote
                pool = twith($pool) do p
                    tforeach(p, $(esc(range))) do $(esc(index))
                        for _tmp_var_dontusename in [0] # To allow `continue` to function in body
                            $(esc(body))
                        end
                    end
                end
                $(esc(log)) ? pool : nothing
            end
        else
            throw(ArgumentError("nested outer loops are not currently supported by thread macro"))
        end
    else
        throw(ArgumentError("unrecognized argument tothread macro"))
    end
end

"""
    @bthreads

Mimics `Base.Threads.@threads, but keeps the iterated tasks off if the primary
thread.`

# Example
```julia
julia> @bthreads for x in 1:8
         println((x, Threads.threadid()))
       end
(1, 2)
(6, 4)
(3, 3)
(7, 4)
(4, 3)
(8, 4)
(5, 3)
(2, 2)
```
Note that execution order is not guaranteed, but the primary thread does not
show up on any of the jobs.
"""
macro bthreads(args...)
    return _pthread_macro(:(StaticPool(2)), false, args...)
    nothing
end

"""
    @qthreads

Mimics `Base.Threads.@threads`, but uses a task queueing strategy, only starting
a new task when an previous one (on any thread) has completed.  This can provide
performance advantages when the iterated tasks are very nonuniform in length.
The primary thread is used.  To prevent usage of the primary thread, see
[`@qbthreads`](@ref).

# Example
```julia
julia> @qthreads for x in 1:8
         println((x, Threads.threadid()))
       end
(2, 4)
(3, 3)
(4, 2)
(5, 4)
(6, 3)
(7, 2)
(8, 4)
(1, 1)
```
Note that execution order is not guaranteed and the primary thread is used.
"""
macro qthreads(args...)
    return _pthread_macro(:(QueuePool(1)), false, args...)
end

"""
    @qbthreads

Mimics `Base.Threads.@threads`, but uses a task queueing strategy, only starting
a new task when an previous one (on any thread) has completed.  This can provide
performance advantages when the iterated tasks are very nonuniform in length.
The primary thread is not used.  To allow usage of the primary thread, see
[`@qthreads`](@ref).

# Example
```julia
julia> @qbthreads for x in 1:8
         println((x, Threads.threadid()))
       end
(2, 4)
(3, 2)
(1, 3)
(4, 4)
(5, 2)
(6, 3)
(7, 4)
(8, 2)
```
Note that execution order is not guaranteed, but the primary thread does not
show up on any of the jobs.
"""
macro qbthreads(args...)
    return _pthread_macro(:(QueuePool(2)), false, args...)
end

"""
    @logthreads -> pool

Mimics `Base.Threads.@threads`.  Returns a logged pool that can be analyzed with
the logging functions and `plot`ted.

# Example
```julia
julia> pool = @logthreads for x in 1:8
         println((x, Threads.threadid()))
       end;
(1, 1)
(5, 3)
(7, 4)
(2, 1)
(6, 3)
(8, 4)
(3, 2)
(4, 2)

julia> plot(pool)
```
Note that execution order is not guaranteed and the primary thread is used.
"""
macro logthreads(args...)
    return _pthread_macro(:(LoggedStaticPool(1)), true, args...)
end

"""
    @logbthreads -> pool

Mimics `Base.Threads.@threads, but keeps the iterated tasks off if the primary
thread.`  Returns a logged pool that can be analyzed with the logging functions
and `plot`ted.

# Example
```julia
julia> pool = @logbthreads for x in 1:8
         println((x, Threads.threadid()))
       end;
(3, 4)
(2, 3)
(1, 2)
(4, 4)
(5, 3)
(6, 2)
(8, 3)
(7, 4)

julia> plot(pool)
```
Note that execution order is not guaranteed, but the primary thread does not
show up on any of the jobs.
"""
macro logbthreads(args...)
    return _pthread_macro(:(LoggedStaticPool(2)), true, args...)
end

"""
    @logqthreads -> pool

Mimics `Base.Threads.@threads`, but uses a task queueing strategy, only starting
a new task when an previous one (on any thread) has completed.  Returns a logged
pool that can be analyzed with the logging functions and `plot`ted. The primary
thread is used.  To prevent usage of the primary thread, see
[`@logqbthreads`](@ref).

# Example
```julia
julia> pool = @logqthreads for x in 1:8
         println((x, Threads.threadid()))
       end;
(1, 1)
(3, 2)
(7, 4)
(5, 3)
(2, 1)
(8, 4)
(6, 3)
(4, 2)

julia> plot(pool)
```
Note that execution order is not guaranteed and the primary thread is used.
"""
macro logqthreads(args...)
    return _pthread_macro(:(LoggedQueuePool(1)), true, args...)
end

"""
    @logqbthreads -> pool

Mimics `Base.Threads.@threads`, but uses a task queueing strategy, only starting
a new task when an previous one (on any thread) has completed.  Returns a logged
pool that can be analyzed with the logging functions and `plot`ted. The primary
thread is not used.  To allow usage of the primary thread, see
[`@logqthreads`](@ref).

# Example
```julia
julia> pool = @logqbthreads for x in 1:8
         println((x, Threads.threadid()))
       end;
(2, 3)
(1, 4)
(3, 2)
(4, 3)
(5, 4)
(6, 2)
(7, 3)
(8, 4)

julia> plot(pool)
```
Note that execution order is not guaranteed, but the primary thread does not
show up on any of the jobs.
"""
macro logqbthreads(args...)
    return _pthread_macro(:(LoggedQueuePool(2)), true, args...)
end

"""
    @tspawnat tid -> task

Mimics `Base.Threads.@spawn`, but assigns the task to thread `tid`.

# Example
```julia
julia> t = @tspawnat 4 Threads.threadid()
Task (runnable) @0x0000000010743c70

julia> fetch(t)
4
```
"""
macro tspawnat(thrdid, expr)
    if VERSION >= v"1.4"
        letargs = Base._lift_one_interp!(expr)

        thunk = esc(:(()->($expr)))
        var = esc(Base.sync_varname)
        tid = esc(thrdid)
        quote
            if $tid < 1 || $tid > Threads.nthreads()
                throw(AssertionError("@tspawnat thread assignment ($($tid)) must be between 1 and Threads.nthreads() (1:$(Threads.nthreads()))"))
            end
            let $(letargs...)
                local task = Task($thunk)
                task.sticky = false
                ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, $tid-1)
                if $(Expr(:islocal, var))
                    put!($var, task)
                end
                schedule(task)
                task
            end
        end
    else
        thunk = esc(:(()->($expr)))
        var = esc(Base.sync_varname)
        tid = esc(thrdid)
        quote
            if $tid < 1 || $tid > Threads.nthreads()
                throw(AssertionError("@tspawnat thread assignment ($($tid)) must be between 1 and Threads.nthreads() (1:$(Threads.nthreads()))"))
            end
            local task = Task($thunk)
            task.sticky = false
            ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, $tid-1)
            if $(Expr(:isdefined, var))
                push!($var, task)
            end
            schedule(task)
            task
        end
    end
end


text:
"""
  coefs_table(mf,lo, up, betass, std_errors)

  This function takes a `ModelFrame` object `mf`, numeric arguments `lo`, `up`, `betass` estimates and `std_errors` (beta's
  standard errors, returns a table with model estimates, 95% CI, and p-values.
"""

coefs_table(mf, lo, up, betass, std_errors)

function coefs_table(; mf = throw(ArgumentError("mf is missing")),
   lo = throw(ArgumentError("lo is missing")),
   up = throw(ArgumentError("up is missing")),
   betass = throw(ArgumentError("betass is missing")),
   std_errors = throw(ArgumentError("std_errors is missing")))
   zz = betass ./ std_errors
   result = (StatsBase.CoefTable(hcat(round(betass,4),
              lo,
              up,
              round(std_errors,4),
              round(zz,4),
              2.0 * Distributions.ccdf(Distributions.Normal(),
              abs.(zz))),
              ["Estimate","2.5%","97.5%","Std.Error","t value", "Pr(>|t|)"],
              ["$i" for i = coefnames(mf)], 4))
  return(result)
end


text:
#
#
# Proximal Point Problem and Options
#
#
export ProximalProblem
export CyclicProximalPointOptions, DouglasRachfordOptions
export get_cost, getProximalMap
export DebugProximalParameter
export RecordProximalParameter

@doc raw"""
    ProximalProblem <: Problem
specify a problem for solvers based on the evaluation of proximal map(s).

# Fields
* `M`            - a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) $\mathcal M$
* `cost` - a function $F\colon\mathcal M\to\mathbb R$ to
  minimize
* `proxes` - proximal maps $\operatorname{prox}_{\lambda\varphi}\colon\mathcal M\to\mathcal M$
  as functions (λ,x) -> y, i.e. the prox parameter λ also belongs to the signature of the proximal map.
* `number_of_proxes` - (length(proxes)) number of proxmal Maps,
  e.g. if one of the maps is a compined one such that the proximal Maps
  functions return more than one entry per function

# See also
[`cyclic_proximal_point`](@ref), [`get_cost`](@ref), [`getProximalMap`](@ref)
"""
mutable struct ProximalProblem{mT <: Manifold} <: Problem
  M::mT
  cost::Function
  proxes::Array{Function,N} where N
  number_of_proxes::Array{Int,1}
  ProximalProblem(M::mT, cF::Function, proxMaps::Array{Function,1}) where {mT <: Manifold}= new{mT}(M,cF,proxMaps,ones(length(proxMaps)))
  ProximalProblem(M::mT, cF::Function, proxMaps::Array{Function,1}, nOP::Array{Int,1}) where {mT <: Manifold} =
    length(nOP) != length(proxMaps) ? throw(ErrorException("The number_of_proxes ($(nOP)) has to be the same length as the number of Proxes ($(length(proxMaps)).")) :
    new{mT}(M,cF,proxMaps,nOP)
end
@doc raw"""
    getProximalMap(p,λ,x,i)

evaluate the `i`th proximal map of `ProximalProblem p` at the point `x` of `p.M` with parameter `λ`$>0$.
"""
function getProximalMap(p::P,λ,x,i) where {P <: ProximalProblem{M} where M <: Manifold}
    if i>length(p.proxes)
        throw( ErrorException("the $(i)th entry does not exists, only $(length(p.proxes)) available.") )
    end
    return p.proxes[i](λ,x);
end
#
#
# Proximal based Options
#
#
"""
    CyclicProximalPointOptions <: Options

stores options for the [`cyclic_proximal_point`](@ref) algorithm. These are the

# Fields
* `x0` – an point to start
* `stopping_criterion` – a function `@(iter,x,xnew,λ_k)` based on the current
    `iter`, `x` and `xnew` as well as the current value of `λ`.
* `λ` – (@(iter) -> 1/iter) a function for the values of λ_k per iteration/cycle
* `evaluationOrder` – ([`LinearEvalOrder`](@ref)`()`) how to cycle through the proximal maps.
    Other values are [`RandomEvalOrder`](@ref)`()` that takes a new random order each
    iteration, and [`FixedRandomEvalOrder`](@ref)`()` that fixes a random cycle for all iterations.

# See also
[`cyclic_proximal_point`](@ref)
"""
mutable struct CyclicProximalPointOptions <: Options
    x
    stop::StoppingCriterion
    λ::Function
    orderType::EvalOrder
    order::Array{Int,1}
end
function CyclicProximalPointOptions(
    x,
    s::StoppingCriterion,
    λ::Function=(iter)-> 1.0/iter,
    o::EvalOrder=LinearEvalOrder()
    )
    return CyclicProximalPointOptions(x,s,λ,o,[])
end
@doc raw"""
    DouglasRachfordOptions <: Options

Store all options required for the DouglasRachford algorithm,

# Fields
* `x` - the current iterate (result) For the parallel Douglas-Rachford, this is
  not a value from the `PowerManifold` manifold but the mean.
* `s` – the last result of the double reflection at the proxes relaxed by `α`.
* `λ` – (`(iter)->1.0`) function to provide the value for the proximal parameter
  during the calls
* `α` – (`(iter)->0.9`) relaxation of the step from old to new iterate, i.e.
  $x^{(k+1)} = g(α(k); x^{(k)}, t^{(k)})$, where $t^{(k)}$ is the result
  of the double reflection involved in the DR algorithm
* `R` – ([`reflect`](@ref)) method employed in the iteration to perform the reflection of `x` at
  the prox `p`.
* `stop` – ([`StopAfterIteration`](@ref)`(300)`) a [`StoppingCriterion`](@ref)
* `parallel` – (`false`) inducate whether we are running a pallel Douglas-Rachford
  or not.
"""
mutable struct DouglasRachfordOptions <: Options
    x
    s
    λ::Function
    α::Function
    R::Function
    stop::StoppingCriterion
    parallel::Bool
    function DouglasRachfordOptions(
        x,
        λ::Function=(iter)->1.0,
        α::Function=(iter)->0.9,
        R=reflect,
        stop::StoppingCriterion = StopAfterIteration(300),
        parallel=false
    )
    return new(x,x,λ,α,R,stop,parallel)
    end
end
#
# Debug
#
# overwrite defaults, since we store the result in the mean field
#
# Debug the Cyclic Proximal point parameter
#
@doc raw"""
    DebugProximalParameter <: DebugAction

print the current iterates proximal point algorithm parameter given by
[`Options`](@ref)s `o.λ`.
"""
mutable struct DebugProximalParameter <: DebugAction
    print::Function
    prefix::String
    DebugProximalParameter(long::Bool=false,print::Function=print) = new(print, long ? "Proximal Map Parameter λ(i):" : "λ:" )
end
(d::DebugProximalParameter)(p::ProximalProblem,o::DouglasRachfordOptions,i::Int) = d.print((i>0) ? d.prefix*string(o.λ(i)) : "")
(d::DebugProximalParameter)(p::ProximalProblem,o::CyclicProximalPointOptions,i::Int) = d.print((i>0) ? d.prefix*string(o.λ(i)) : "")

#
# Record
@doc raw"""
    RecordProximalParameter <: RecordAction

recoed the current iterates proximal point algorithm parameter given by in
[`Options`](@ref)s `o.λ`.
"""
mutable struct RecordProximalParameter <: RecordAction
    recordedValues::Array{Float64,1}
    RecordProximalParameter() = new(Array{Float64,1}())
end
(r::RecordProximalParameter)(p::P,o::O,i::Int) where {P <: ProximalProblem, O <: CyclicProximalPointOptions} = record_or_eset!(r, o.λ(i), i)
(r::RecordProximalParameter)(p::P,o::O,i::Int) where {P <: ProximalProblem, O <: DouglasRachfordOptions} = record_or_eset!(r, o.λ(i), i)


text:
# using Memoization
using Printf

function run_euler(p)
    # GC.gc(false)
    t = @timed p()
    # GC.gc(true)
  @printf "%s: %d: %2.8fs\n" p t.value t.time
end

# @memoize function fib(n)
#    n < 3 ? 1 : fib(n-1) + fib(n-2)
# end

function fib_nonmem(n)
    n < 3 ? 1 : fib(n-1) + fib(n-2)
end

# For Fibonacci with BigInt
# Note: @memoize is too slow, probably because of BigInts
# (or that I have not used @memoize correctly)
fibmemBigInt_dict = Dict{Int,BigInt}()
function fibmemBigInt(n)
    get!(fibmemBigInt_dict, n) do
        n < 3 ? 1 : fibmemBigInt(BigInt(n-1)) + fibmemBigInt(BigInt(n-2))
    end
end

function isPrime(n::Int)
    if n < 2
        return false
    end
    if n === 2 || n === 3
        return true
    end

    if n % 2 === 0
        return false
    end

    m = 1+ceil(Int,sqrt(n))
    for i in 3:2:m
        if n % i === 0
            return false;
        end
    end
    return true;
end

# @memoize function isPrimeCached(n)
#    return isPrime(n)
# end

# return all primes <= n
function primes(n::Int)
    primes = [2]
    for i in 3:2:n
        if isPrime(i)
            push!(primes,i)
        end
    end
    return primes;
end

#
# returns the prime divisors of n.
#

function prime_divisors(n::Int)
    if n === 1
        return []
    end
    divisors = []
    while (n > 1)
        m = round(Int,ceil(sqrt(n)))
        for i in 2:m
            if n % i == 0 && isPrime(i)
                push!(divisors,i)
                n /= i
                continue
            end
        end
    end
    return divisors
end

# Find the primes <= n
function sieve(n::Int)
    a = Array(range(1,length=n))
    a[1] = 0;
    a[2] = 2;
    for i in 2:ceil(Int,sqrt(n))
        for j in i*i:i:n
            a[j] = 0
        end
    end
    return filter(i->a[i]>0,1:n);

end

# Prime factors of n
function factors(n::Int)
    if n == 1
        return [1]
    end
    f = []
    while n % 2 == 0
        push!(f,2)
        n /= 2
    end
    t = 3
    while n > 1 && t < ceil(sqrt(n))
        while n % t == 0
            push!(f,t)
            n /= t
        end
        t += 2
    end
    if n > 1
        push!(f,round(Int,n))
    end
    return f
end

# Is n a palindromic number?
function palindromic_number(n::Int)
    s = string(n)
    len = length(s)
    for i in 1:ceil(Int,len/2)
        if s[i] != s[len-i+1]
            return false
        end
    end
    return true
end

# is the list a palindrom?
function palindromic_list(s)
    # return s == reverse(s)
    len = length(s)
    for i in 1:ceil(Int,len/2)
        if s[i] != s[len-i+1]
            return false
        end
    end
    return true
end

# All divisors of a number, except 1 and n
function all_divisors(n::Int)
    divisors = [];
    m = round(Int,n / 2)
    for i in 2:m
        if n % i == 0
            push!(divisors,i)
        end
    end
    return divisors;

end

# // All divisors of n, including 1 and n
function all_divisors2(n::Int)
    return hcat([1], all_divisors(n),[n])
end

# All divisors of n, including 1 (but not n)
function all_divisors3(n::Int)
    divisors = [1]
    m = round(Int,n / 2)
    for i in 2:m
        if n % i == 0
            push!(divisors,i)
        end
    end
    return divisors
end

# convert a list of elements in a to a hash with
#   { element: number of occurrences of element, ...}
function collect_list(a)
    m = Dict()
    a.|>(e->
        if !haskey(m,e)
            m[e] = 1
        else
            m[e] += 1
        end

    )
    return m
end

# next permutation of array p
function next_permutation(p)
    perm = p
    n = length(perm)
    k = n - 1
    while perm[k] > perm[k+1] && k >= 0
        k = k - 1
    end
    if k > 0
        j = n
        while perm[k] > perm[j]
            j = j - 1
        end
        perm[k],perm[j] = perm[j], perm[k]
        r = n
        s = k + 1
        while r > s
            perm[r],perm[s] = perm[s],perm[r]
            r = r - 1
            s = s + 1
        end
    end
    return perm
end

#
# Generate all permutations
# This is much faster than all_permutations2
#
function all_permutations(s)
    perms = []
    push!(perms,s)
    ss = copy(s)
    p = next_permutation(copy(s));;
    push!(perms,p)
    while true
        p = next_permutation(p)
        # Note: I have to copy p here
        push!(perms,copy(p))
        if p == reverse(s)
            break
        end
    end
    return perms;
end

# Inspired by an answer from
# http://stackoverflow.com/questions/4240080/generating-all-permutations-of-a-given-string
# This is too slow. Use all_permutations instead
function all_permutations2(s)
    perms = [];
    permutation_tmp([], s, perms);
    return perms;
end

function permutation_tmp(prefix,s, perms)
    n = length(s)
    if n == 0
        push!(perms,prefix)
    else
        for i in 1:n
            permutation_tmp(vcat(prefix,s[i]), vcat(s[1:i-1],s[i+1:n]), perms)
        end
    end
end

# Convert a number to its (integer) digits
# Note: this should be replaced with
#  digits(n)
function num_to_list(n::Int)
    # return split(string(n),"").|>i->parse(Int,i)
    return digits(n)
end

# length of a number (i.e. the number of digits)
function nlen(n::Int)
    return floor(Int,log10(n))+1
end

# check is s is a pandigital number  1..9
function is_pandigital(s)
    return length(s) === 9 && !occursin("0",s) && length(Set(s)) == 9
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1Event <: SwaggerModel
    action::Any # spec type: Union{ Nothing, String } # spec name: action
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count
    eventTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1MicroTime } # spec name: eventTime
    firstTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: firstTimestamp
    involvedObject::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: involvedObject
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    lastTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastTimestamp
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    related::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: related
    reportingComponent::Any # spec type: Union{ Nothing, String } # spec name: reportingComponent
    reportingInstance::Any # spec type: Union{ Nothing, String } # spec name: reportingInstance
    series::Any # spec type: Union{ Nothing, IoK8sApiCoreV1EventSeries } # spec name: series
    source::Any # spec type: Union{ Nothing, IoK8sApiCoreV1EventSource } # spec name: source
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiCoreV1Event(;action=nothing, apiVersion=nothing, count=nothing, eventTime=nothing, firstTimestamp=nothing, involvedObject=nothing, kind=nothing, lastTimestamp=nothing, message=nothing, metadata=nothing, reason=nothing, related=nothing, reportingComponent=nothing, reportingInstance=nothing, series=nothing, source=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Event, Symbol("action"), action)
        setfield!(o, Symbol("action"), action)
        validate_property(IoK8sApiCoreV1Event, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1Event, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        validate_property(IoK8sApiCoreV1Event, Symbol("eventTime"), eventTime)
        setfield!(o, Symbol("eventTime"), eventTime)
        validate_property(IoK8sApiCoreV1Event, Symbol("firstTimestamp"), firstTimestamp)
        setfield!(o, Symbol("firstTimestamp"), firstTimestamp)
        validate_property(IoK8sApiCoreV1Event, Symbol("involvedObject"), involvedObject)
        setfield!(o, Symbol("involvedObject"), involvedObject)
        validate_property(IoK8sApiCoreV1Event, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1Event, Symbol("lastTimestamp"), lastTimestamp)
        setfield!(o, Symbol("lastTimestamp"), lastTimestamp)
        validate_property(IoK8sApiCoreV1Event, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1Event, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1Event, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiCoreV1Event, Symbol("related"), related)
        setfield!(o, Symbol("related"), related)
        validate_property(IoK8sApiCoreV1Event, Symbol("reportingComponent"), reportingComponent)
        setfield!(o, Symbol("reportingComponent"), reportingComponent)
        validate_property(IoK8sApiCoreV1Event, Symbol("reportingInstance"), reportingInstance)
        setfield!(o, Symbol("reportingInstance"), reportingInstance)
        validate_property(IoK8sApiCoreV1Event, Symbol("series"), series)
        setfield!(o, Symbol("series"), series)
        validate_property(IoK8sApiCoreV1Event, Symbol("source"), source)
        setfield!(o, Symbol("source"), source)
        validate_property(IoK8sApiCoreV1Event, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiCoreV1Event

const _property_map_IoK8sApiCoreV1Event = Dict{Symbol,Symbol}(Symbol("action")=>Symbol("action"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("count")=>Symbol("count"), Symbol("eventTime")=>Symbol("eventTime"), Symbol("firstTimestamp")=>Symbol("firstTimestamp"), Symbol("involvedObject")=>Symbol("involvedObject"), Symbol("kind")=>Symbol("kind"), Symbol("lastTimestamp")=>Symbol("lastTimestamp"), Symbol("message")=>Symbol("message"), Symbol("metadata")=>Symbol("metadata"), Symbol("reason")=>Symbol("reason"), Symbol("related")=>Symbol("related"), Symbol("reportingComponent")=>Symbol("reportingComponent"), Symbol("reportingInstance")=>Symbol("reportingInstance"), Symbol("series")=>Symbol("series"), Symbol("source")=>Symbol("source"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiCoreV1Event = Dict{Symbol,String}(Symbol("action")=>"String", Symbol("apiVersion")=>"String", Symbol("count")=>"Int32", Symbol("eventTime")=>"IoK8sApimachineryPkgApisMetaV1MicroTime", Symbol("firstTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("involvedObject")=>"IoK8sApiCoreV1ObjectReference", Symbol("kind")=>"String", Symbol("lastTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("reason")=>"String", Symbol("related")=>"IoK8sApiCoreV1ObjectReference", Symbol("reportingComponent")=>"String", Symbol("reportingInstance")=>"String", Symbol("series")=>"IoK8sApiCoreV1EventSeries", Symbol("source")=>"IoK8sApiCoreV1EventSource", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1Event }) = collect(keys(_property_map_IoK8sApiCoreV1Event))
Swagger.property_type(::Type{ IoK8sApiCoreV1Event }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1Event[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Event }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Event[property_name]

function check_required(o::IoK8sApiCoreV1Event)
    (getproperty(o, Symbol("involvedObject")) === nothing) && (return false)
    (getproperty(o, Symbol("metadata")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Event }, name::Symbol, val)
end


text:
function plotFunc_timecourse(Sim::Module)
    rc("figure",figsize = (12,8));
    rc("font",size = 18);
    rc("axes",linewidth = 1);
    rc("lines",linewidth = 2);

    subplots_adjust(wspace=0.5, hspace=0.5);

    subplot(2,2,1);
    plot(Sim.t,Sim.CBMa./maximum(Sim.CBMa).*100,"k");
    xlim(-0.5,9.5);
    xticks([0,3,6,9]);
    ylim(0,120);
    yticks([0,20,40,60,80,100]);
    xlabel("Time (min)");
    ylabel("CP668 activity \n (relative % max.)");

    subplot(2,2,2);
    plot(Sim.t,Sim.TAK1a./maximum(Sim.TAK1a).*100,"k");
    xlim(-0.5,9.5);
    xticks([0,3,6,9]);
    ylim(0,120);
    yticks([0,20,40,60,80,100]);
    xlabel("Time (min)");
    ylabel("TAK1 activity \n (relative % max.)");

    subplot(2,2,3);
    plot(Sim.t,Sim.IKKa./maximum(Sim.IKKa).*100,"k");
    xlim(-0.5,9.5);
    xticks([0,3,6,9]);
    ylim(0,120);
    yticks([0,20,40,60,80,100]);
    xlabel("Time (min)");
    ylabel("IKK activity \n (relative % max.)");

    subplot(2,2,4);
    plot(Sim.t,Sim.NFkBa./maximum(Sim.NFkBa).*100,"k");
    xlim(0,360);
    xticks([0,120,240,360]);
    ylim(0,120);
    yticks([0,20,40,60,80,100]);
    xlabel("Time (min)");
    ylabel("NFκB activity \n (relative % max.)");

    show();
end


text:
"""
parse_tuple(ps, ret)

`ret` is followed by a comma so tries to parse the rest of the
tuple.
"""
function parse_tuple end

iserrorwrapped(x, head) = x.head === head || (x.head === :errortoken && length(x.args) > 0 && iserrorwrapped(x.args[1], head))

function parse_tuple(ps::ParseState, ret::EXPR)
    op = EXPR(next(ps))
    if istuple(ret) && !(length(ret.trivia) > 1 && iserrorwrapped(last(ret.trivia), :RPAREN))
        if (isassignmentop(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
            pushtotrivia!(ret, op)
        elseif closer(ps)
            @static if VERSION > v"1.1-"
                pushtotrivia!(ret, mErrorToken(ps, op, Unknown))
            else
                pushtotrivia!(ret, op)
            end
        else
            nextarg = @closer ps :tuple parse_expression(ps)
            if !(is_lparen(first(ret.trivia)))
                pushtotrivia!(ret, op)
                push!(ret, nextarg)
            else
                ret = EXPR(:tuple, EXPR[ret, nextarg], EXPR[op])
            end
        end
    else
        if (isassignmentop(ps.nt) && kindof(ps.nt) != Tokens.APPROX)
            ret = EXPR(:tuple, EXPR[ret], EXPR[op])
        elseif closer(ps)
            @static if VERSION > v"1.1-"
                ret = mErrorToken(ps, EXPR(:tuple, EXPR[ret], EXPR[op]), Unknown)
            else
                ret = EXPR(:tuple, EXPR[ret], EXPR[op])
            end
        else
            nextarg = @closer ps :tuple parse_expression(ps)
            ret = EXPR(:tuple, EXPR[ret, nextarg], EXPR[op])
        end
    end
    return ret
end

"""
    parse_array(ps)
Having hit '[' return either:
+ A vect
+ A vcat
+ A comprehension
+ An array (vcat of hcats)
"""
function parse_array(ps::ParseState, isref = false)
    args = EXPR[]
    trivia = EXPR[EXPR(ps)]
    if kindof(ps.nt) === Tokens.RSQUARE
        push!(trivia, accept_rsquare(ps))
        ret = EXPR(:vect, args, trivia)
    else
        first_arg = @nocloser ps :newline @closesquare ps  @closer ps :insquare @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)
        if isref && _do_kw_convert(ps, first_arg)
            first_arg = _kw_convert(first_arg)
        end

        if kindof(ps.nt) === Tokens.RSQUARE
            if headof(first_arg) === :generator || headof(first_arg) === :flatten
                push!(trivia, accept_rsquare(ps))
                if isbinarycall(first_arg.args[1]) && is_pairarrow(first_arg.args[1].args[2])
                    return EXPR(:dict_comprehension, EXPR[first_arg], trivia)
                else
                    return EXPR(:comprehension, EXPR[first_arg], trivia)
                end
            elseif kindof(ps.ws) == SemiColonWS
                push!(args, first_arg)
                push!(trivia, accept_rsquare(ps))
                return EXPR(:vcat, args, trivia)
            else
                push!(args, first_arg)
                push!(trivia, accept_rsquare(ps))
                ret = EXPR(:vect, args, trivia)
            end
        elseif iscomma(ps.nt)
            push!(args, first_arg)
            push!(trivia, accept_comma(ps))
            @closesquare ps parse_comma_sep(ps, args, trivia, isref, insert_params_at = 1)
            push!(trivia, accept_rsquare(ps))
            return EXPR(:vect, args, trivia)
        elseif kindof(ps.ws) == WS || kindof(ps.ws) == SemiColonWS || kindof(ps.ws) == NewLineWS
            ps.closer.inref = false
            args1 = EXPR[first_arg]
            prevpos = position(ps)
            while kindof(ps.nt) !== Tokens.RSQUARE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                a = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                push!(args1, a)

                prevpos = loop_check(ps, prevpos)
            end
            if kindof(ps.nt) === Tokens.RSQUARE && kindof(ps.ws) != SemiColonWS
                push!(trivia, accept_rsquare(ps))
                if length(args1) == 1
                    return EXPR(:vcat, args1, trivia)
                else
                    return EXPR(:hcat, args1, trivia)
                end
            else
                if length(args1) == 1
                    first_row = args1[1]
                else
                    first_row = EXPR(:row, args1, nothing)
                end
                push!(args, first_row)
                prevpos = position(ps)
                while kindof(ps.nt) !== Tokens.RSQUARE && kindof(ps.nt) !== Tokens.ENDMARKER
                    first_arg = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                    args2 = EXPR[first_arg]
                    prevpos1 = position(ps)
                    while kindof(ps.nt) !== Tokens.RSQUARE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                        a = @closesquare ps @closer ps :ws @closer ps :wsop parse_expression(ps)

                        push!(args2, a)
                        prevpos1 = loop_check(ps, prevpos1)
                    end
                    # if only one entry dont use :row
                    if length(args2) == 1
                        push!(args, args2[1])
                    else
                        push!(args, EXPR(:row, args2, nothing))
                    end
                    prevpos = loop_check(ps, prevpos)
                end
                push!(trivia, accept_rsquare(ps))
                return EXPR(:vcat, args, trivia)
            end
        else
            push!(args, first_arg)
            push!(trivia, accept_rsquare(ps))
            ret = EXPR(:vect, args, trivia)
        end
    end
    return ret
end

function parse_array_row()

end

"""
    parse_ref(ps, ret)

Handles cases where an expression - `ret` - is followed by
`[`. Parses the following bracketed expression and modifies it's
`.head` appropriately.
"""
function parse_ref(ps::ParseState, ret::EXPR)
    next(ps)
    ref = @closer ps :inref @nocloser ps :inwhere parse_array(ps, true)
    if headof(ref) === :vect
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:ref, args, ref.trivia)
    elseif headof(ref) === :hcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_hcat, args, ref.trivia)
    elseif headof(ref) === :vcat
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_vcat, args, ref.trivia)
    else
        args = EXPR[ret]
        for a in ref.args
            push!(args, a)
        end
        return EXPR(:typed_comprehension, args, ref.trivia)
    end
end

"""
parse_curly(ps, ret)

Parses the juxtaposition of `ret` with an opening brace. Parses a comma
seperated list.
"""
function parse_curly(ps::ParseState, ret::EXPR)
    args = EXPR[ret]
    trivia = EXPR[EXPR(next(ps))]
    parse_comma_sep(ps, args, trivia, true)
    accept_rbrace(ps, trivia)
    return EXPR(:curly, args, trivia)
end

function parse_braces(ps::ParseState)
    return @default ps @nocloser ps :inwhere parse_barray(ps)
end

function parse_barray(ps::ParseState)
    args = EXPR[]
    trivia = EXPR[EXPR(ps)]

    if kindof(ps.nt) === Tokens.RBRACE
        accept_rbrace(ps, trivia)
        ret = EXPR(:braces, args, trivia)
    else
        first_arg = @nocloser ps :newline @closebrace ps  @closer ps :ws @closer ps :wsop @closer ps :comma parse_expression(ps)
        if kindof(ps.nt) === Tokens.RBRACE
            push!(args, first_arg)
            if kindof(ps.ws) == SemiColonWS
                pushfirst!(args, EXPR(:parameters, EXPR[], nothing))
            end
            accept_rbrace(ps, trivia)
            ret = EXPR(:braces, args, trivia)
        elseif iscomma(ps.nt)
            push!(args, first_arg)
            accept_comma(ps, trivia)
            @closebrace ps parse_comma_sep(ps, args, trivia, true, insert_params_at = 1)
            accept_rbrace(ps, trivia)
            return EXPR(:braces, args, trivia)
        elseif kindof(ps.ws) == NewLineWS
            ret = EXPR(:bracescat, args, trivia)
            push!(ret, first_arg)
            prevpos = position(ps)
            while kindof(ps.nt) != Tokens.RBRACE && kindof(ps.nt) !== Tokens.ENDMARKER
                a = @closebrace ps  parse_expression(ps)
                push!(ret, a)
                prevpos = loop_check(ps, prevpos)
            end
            pushtotrivia!(ret, accept_rbrace(ps))
            return ret
        elseif kindof(ps.ws) == WS || kindof(ps.ws) == SemiColonWS
            first_row = EXPR(:row, EXPR[first_arg])
            prevpos = position(ps)
            while kindof(ps.nt) !== Tokens.RBRACE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                push!(first_row, a)
                prevpos = loop_check(ps, prevpos)
            end
            if kindof(ps.nt) === Tokens.RBRACE && kindof(ps.ws) != SemiColonWS
                if length(first_row.args) == 1
                    first_row = EXPR(:bracescat, first_row.args)
                end
                push!(args, first_row)
                push!(trivia, INSTANCE(next(ps)))
                return EXPR(:bracescat, args, trivia)
            else
                if length(first_row.args) == 1
                    first_row = first_row.args[1]
                else
                    first_row = EXPR(:row, first_row.args)
                end
                ret = EXPR(:bracescat, EXPR[first_row], trivia)
                prevpos = position(ps)
                while kindof(ps.nt) != Tokens.RBRACE
                    if kindof(ps.nt) === Tokens.ENDMARKER
                        break
                    end
                    first_arg = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                    push!(ret, EXPR(:row, EXPR[first_arg]))
                    prevpos1 = position(ps)
                    while kindof(ps.nt) !== Tokens.RBRACE && kindof(ps.ws) !== NewLineWS && kindof(ps.ws) !== SemiColonWS && kindof(ps.nt) !== Tokens.ENDMARKER
                        a = @closebrace ps @closer ps :ws @closer ps :wsop parse_expression(ps)
                        push!(last(ret.args), a)
                        prevpos1 = loop_check(ps, prevpos1)
                    end
                    # if only one entry dont use :row
                    if length(last(ret.args).args) == 1
                        ret.args[end] = setparent!(ret.args[end].args[1], ret)
                    end
                    update_span!(ret)
                    prevpos = loop_check(ps, prevpos)
                end
                pushtotrivia!(ret, accept_rbrace(ps))
                update_span!(ret)
                return ret
            end
        else
            ret = EXPR(:braces, args, trivia)
            push!(ret, first_arg)
            pushtotrivia!(ret, accept_rbrace(ps))
        end
    end
    return ret
end


text:
using Memoize, IterTools, LinearAlgebra

@memoize function mem_det(A::Array{Int,2})
    return round(Int, det(A))
end

"""
`istum(A)` determines if the integer matrix `A`
is totally unimodular. This may be called with an
optional second argument `istum(A,true)` that,
in case the matrix is not totally unimodular, will
report the first submatrix found whose determinant
is not in `{-1,0,1}`.

**Warning**: This is not an efficient function.
"""
function istum(A::Array{Int,2}, verbose::Bool = false)
    r, c = size(A)
    n = min(r, c)

    for k = 1:n
        rows = subsets(collect(1:r), k)
        cols = subsets(collect(1:c), k)
        for x in rows
            for y in cols
                B = A[x, y]
                d = mem_det(B)
                if abs(d) > 1
                    if verbose
                        println("Rows: ", x)
                        println("Cols: ", y)
                        println("$B has determinant $d")
                    end
                    return false
                end
            end
        end
    end
    return true
end


text:
ENV["MATLAB_HOME"] = "D:\\Matlab\\R2020b"
using MATLAB
using BSON

s = MSession()

X = BSON.load("data.bson")[:X]
X = mxarray(X)
@mput X

eval_string(s, "MakeVideo(X, 30, \"Video.gif\");")

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse