text:
@resumable function readcards(io::IO)::String
    iscontinueline, card = cardpart(readnoncommentline(io))
    nextcard = ""
    while ~eof(io)
        iscontinueline, nextcard = cardpart(readnoncommentline(io))
        if iscontinueline
            card = card * nextcard
        else
            @yield card
            card, nextcard = nextcard, card
        end
    end
    @yield card
    nothing # this is nessessary, why?
end

function readnoncommentline(io::IO)
    unwantedline = true
    line = ""
    while unwantedline && ~eof(io)
        line = readline(io)
        unwantedline = occursin(r"^\s*(?:\*)",line) || occursin(r"^\s*$",line)
    end
    return unwantedline ? "" : line
end

function cardpart(line)
    pos = 1
    m  = match(r"^\s*(\+)?(.*)"i,line,pos)
    iscontinueline = m.captures[1] == "+"
    pos = m.offsets[2]
    cardstart = pos
    while (m = match(r".*?([\"']).*?\1(.*)"i,line,pos)) != nothing
        pos = m.offsets[2]
    end
    m = match(r"(.)(?:(?:\$)|(?://))"i,line,pos)
    if m != nothing
        return (iscontinueline,line[cardstart:m.offset])
    else
        return (iscontinueline,line[cardstart:end])
    end
end


text:
using LinearAlgebra

# define the basis functions for mean functions
# return matrix X (n by p) such that X_ij = f_j(x_i)
function meanbasisMat(x, p)
    # for now we only consider constant case
    @assert p == 1
    n = size(x, 1)
    X = ones(n)
    return X
end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# We output a message after loading each package to work around a
# ten-minute timeout limit on Travis. Travis assumes the tests have hung
# if the interval between printing something to stdio exceeds ten minutes.
using Images
@info "Finished loading Images package."
using StaticArrays
@info "Finished loading StaticArrays package."

@testset "Visualization" begin
    @info "Running flow visualization test."
    flow = Array{SVector{2, Float64}, 2}(undef,(3, 2))

    flow[1] = SVector{2, Float64}(0.5, 0.5)
    flow[2] = SVector{2, Float64}(-0.4, -0.8)
    flow[3] = SVector{2, Float64}(-0.5, 0.3)
    flow[4] = SVector{2, Float64}(-0.4, -0.8)
    flow[5] = SVector{2, Float64}(-0.5, 0.3)
    flow[6] = SVector{2, Float64}(0.5, 0.5)

    @time hsv = visualize_flow(flow, ColorBased(), RasterConvention())

    @test hsv[1] == hsv[6]
    @test hsv[2] == hsv[4]
    @test hsv[3] == hsv[5]

end;


text:
# tests for kinetic_monte_carlo.jl
using MonteCarloX
using StatsBase
using Random

# for kldivergence
include("utils.jl")

function test_kmc_next(;verbose = false)
    pass = true

    list_rates = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
    sum_rates = sum(list_rates)
    weights = ProbabilityWeights(list_rates)
    tbin = 0.5 / sum_rates
    pdf_time(t) = sum_rates * exp(-(t + tbin) * sum_rates)
    cdf_time(t) = 1.0 - exp(-(t + tbin) * sum_rates)
    pdf_event(e) = list_rates[e] / sum_rates

    if verbose
        println("... check sum(weights)[$(sum(weights))] == sum_rates [$(sum_rates)]")
    end
    pass &= abs(sum_rates - sum(weights)) < 1e-10

    samples = 10^5
    times = zeros(samples)
    events = zeros(samples)
    alg = KineticMonteCarlo()
    rng = MersenneTwister(1000)
    for i = 1:samples
        time, event = next(alg, rng, weights)
        times[i] = time
        events[i] = event
    end

    t_max = maximum(times)
    hist_time = fit(Histogram, times, 0:tbin:t_max + tbin)
    cdf_time_meas = normalize(float(hist_time), mode = :probability)
    pdf_sum = 0.0
    for i in 1:length(cdf_time_meas.weights)
        pdf_sum += cdf_time_meas.weights[i]
        cdf_time_meas.weights[i] = pdf_sum
    end
    # since measurements only performed until time_max, we have to rescale
    # pdf_time(t) to integrate to 1 in the bounds (-infty, t_max). This is
    # achieved by deviding by cdf(t_max). This then percolates to cdf(t)
    cdf_time_renorm(t) = cdf_time(t) / cdf_time(cdf_time_meas.edges[1][end])
    kld_time = kldivergence(cdf_time_meas, cdf_time_renorm)
    pass &= abs(kld_time) < 0.01
    if verbose
        println("... kldivergence cdf_time to cdf_true = $(kld_time)")
    end

    hist_event = fit(Histogram, events, 1:(length(list_rates) + 1))
    pdf_event_meas = normalize(float(hist_event))
    kld_event = kldivergence(pdf_event_meas, pdf_event)
    if verbose
        println("... kldivergence pdf_event to pdf_true = $(kld_event)")
    end
    pass &= abs(kld_event) < 0.01

    return pass
end

# TODO
function test_kmc_event_handler(;verbose = false)
    pass = true

    return pass
end


text:
"""

```
mddf_frame!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)
```

Computes the MDDF for a single frame, modifies the data in the `R` (type `Result`) structure.

"""
function mddf_frame!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)

    # Simplify code by assigning some shortened names
    trajectory = framedata.trajectory
    volume_frame = framedata.volume_frame
    rdf_count_random_frame = framedata.rdf_count_random_frame
    md_count_random_frame = framedata.md_count_random_frame
    dc = framedata.dc
    dmin_mol = framedata.dmin_mol
    dref_mol = framedata.dref_mol
    x_solvent_random = framedata.x_solvent_random
    lc_solvent = framedata.lc_solvent
    solute = trajectory.solute
    solvent = trajectory.solvent
    x_solute = trajectory.x_solute
    x_solvent = trajectory.x_solvent

    # Reset counters for this frame
    reset!(volume_frame)
    @. rdf_count_random_frame = 0.0
    @. md_count_random_frame = 0.0

    # get pbc sides in this frame
    sides = getsides(trajectory, iframe)

    volume_frame.total = sides[1] * sides[2] * sides[3]
    R.volume.total = R.volume.total + volume_frame.total

    R.density.solute = R.density.solute + (solute.nmols / volume_frame.total)
    R.density.solvent = R.density.solvent + (solvent.nmols / volume_frame.total)

    # Add the box side information to the box structure, in this frame
    box = Box(options.lcell, sides, R.cutoff)

    # Will wrap everthing relative to the reference atom of the first molecule
    # and move everything such that that center is in the origin. This is important
    # to simplify the computation of cell indexes, as the minimum coordinates are
    # automatically -side/2 at each direction. The reason for wrapping to the solute
    # center and not to the origin is that if the solute is a protein, for example,
    # the protein gets centered, and hardly ever it will be necessary to wrap the
    # coordinates afterwards, accelerating the computation.
    solute_center = x_solute[R.irefatom]
    wrap!(x_solute, sides, solute_center)
    center_to_origin!(x_solute, solute_center)
    wrap!(x_solvent, sides, solute_center)
    center_to_origin!(x_solvent, solute_center)

    # Initialize linked cells
    initcells!(x_solvent, box, lc_solvent)

    # Check if the cutoff is not too large considering the periodic cell size
    if R.cutoff > sides[1] / 2.0 || R.cutoff > sides[2] / 2.0 || R.cutoff > sides[3] / 2.0
        error("""
        in MDDF: cutoff or dbulk > periodic_dimension/2 in frame: $iframe
                 max(cutoff,dbulk) = $(R.cutoff)
                 sides read from file = $(sides)
                 If sides are zero it is likely that the PBC information is not available
                 in the trajectory file.
        """)
    end

    local n_dmin_in_bulk
    n_solvent_in_bulk = 0.0
    for isolute = 1:solute.nmols

        # We need to do this one solute molecule at a time to avoid exploding the memory
        # requirements
        x_this_solute = viewmol(isolute, x_solute, solute)

        # Compute all distances between solute and solvent atoms which are smaller than the
        # cutoff (this is the most computationally expensive part), the distances are returned
        # in the dc structure
        cutoffdistances!(R.cutoff, x_this_solute, x_solvent, lc_solvent, box, dc)

        # For each solute molecule, update the counters (this is highly suboptimal, because
        # within updatecounters there are loops over solvent molecules, in such a way that
        # this will loop with cost nsolute*nsolvent. However, I cannot see an easy solution
        # at this point with acceptable memory requirements
        n_dmin_in_bulk, n_dref_in_bulk =
            updatecounters!(R, solute, solvent, dc, dmin_mol, dref_mol)
        n_solvent_in_bulk += n_dref_in_bulk
    end
    n_solvent_in_bulk = n_solvent_in_bulk / solute.nmols

    #
    # Computing the random-solvent distribution to compute the random minimum-distance count
    #
    bulk_range = (solvent.nmols-n_dmin_in_bulk+1):solvent.nmols
    for i = 1:options.n_random_samples

        # generate random solvent box, and store it in x_solvent_random
        for j = 1:solvent.nmols
            # Choose randomly one molecule from the bulk, if there are actually bulk molecules
            if n_dmin_in_bulk != 0
                jmol = dmin_mol[random(RNG, bulk_range)].jmol
            else
                jmol = random(RNG, 1:solvent.nmols)
            end
            # Indexes of this molecule in the x_solvent array
            x_ref = viewmol(jmol, x_solvent, solvent)
            # Indexes of the random molecule in random array
            x_rnd = viewmol(j, x_solvent_random, solvent)
            # Generate new random coordinates (translation and rotation) for this molecule
            random_move!(x_ref, R.irefatom, sides, x_rnd, RNG)
        end

        # wrap random solvent coordinates to box, with the center at the origin
        wrap!(x_solvent_random, sides)

        # Initialize linked cells
        initcells!(x_solvent_random, box, lc_solvent)

        # Choose randomly one solute molecule to be the solute in this sample
        i_rand_mol = random(RNG, 1:solute.nmols)
        x_this_solute = viewmol(i_rand_mol, x_solute, solute)

        # Compute all distances between solute and solvent atoms which are smaller than the
        # cutoff (this is the most computationally expensive part), the distances are returned
        # in the dc structure
        cutoffdistances!(R.cutoff, x_this_solute, x_solvent_random, lc_solvent, box, dc)

        # Update the counters of the random distribution
        updatecounters!(
            R,
            rdf_count_random_frame,
            md_count_random_frame,
            solvent,
            dc,
            dmin_mol,
            dref_mol,
        )

    end # random solvent sampling

    # Update counters with the data of this frame
    update_counters_frame!(
        R,
        rdf_count_random_frame,
        md_count_random_frame,
        volume_frame,
        solute,
        solvent,
        n_solvent_in_bulk,
    )

    nothing
end


text:
using StatsBase
using LinearAlgebra, Random, Test

@testset "Cronbach's Alpha" begin
    # basic vanilla test
    cov_X = [10 6 6 6;
             6 11 6 6;
             6 6 12 6;
             6 6 6 13]
    cronbach_X = cronbachalpha(cov_X)
    @test cronbach_X isa CronbachAlpha{Float64}
    @test cronbach_X.alpha ≈ 0.8135593220338981
    @test cronbach_X.dropped ≈
        [0.75, 0.7605633802816901, 0.7714285714285715, 0.782608695652174]

    # testing Rational
    cov_rational = cov_X .// 1
    cronbach_rational = cronbachalpha(cov_rational)
    @test cronbach_rational isa CronbachAlpha{Rational{Int}}
    @test cronbach_rational.alpha == 48 // 59
    @test cronbach_rational.dropped ==
        [3 // 4, 54 // 71, 27 // 35, 18 // 23]

    # testing BigFloat
    cov_bigfloat = BigFloat.(cov_X)
    cronbach_bigfloat = cronbachalpha(cov_bigfloat)
    @test cronbach_bigfloat isa CronbachAlpha{BigFloat}
    @test cronbach_bigfloat.alpha ≈ 0.8135593220338981
    @test cronbach_bigfloat.dropped ≈
        [0.75, 0.7605633802816901, 0.7714285714285715, 0.782608695652174]

    # testing corner cases
    @test_throws MethodError cronbachalpha([1.0, 2.0])
    cov_k2 = [10 6;
              6 11]
    cronbach_k2 = cronbachalpha(cov_k2)
    @test cronbach_k2.alpha ≈ 0.7272727272727273
    @test isempty(cronbach_k2.dropped)

    # testing when Matrix is not positive-definite
    cov_not_pos = [-1 1;
                   -1 1]
    @test_throws ArgumentError cronbachalpha(cov_not_pos)

    # testing with a zero
    cov_zero = [1 2;
                0 1]
    @test_throws ArgumentError cronbachalpha(cov_not_pos)

    # testing with one column
    cov_k1 = reshape([1, 2], 2, 1)
    @test_throws ArgumentError cronbachalpha(cov_k1)

    # testing with Missing
    cov_missing = [1 2;
                   missing 1]
    @test_throws MethodError cronbachalpha(cov_missing)

    # testing Base.show
    cronbach_X = cronbachalpha(cov_X)
    io = IOBuffer()
    show(io, cronbach_X)
    str = String(take!(io))
    @test str == """
        Cronbach's alpha for all items: 0.8136

        Cronbach's alpha if an item is dropped:
        item 1: 0.7500
        item 2: 0.7606
        item 3: 0.7714
        item 4: 0.7826
        """
    # for two columns
    io = IOBuffer()
    show(io, cronbach_k2)
    str = String(take!(io))
    @test str == "Cronbach's alpha for all items: 0.7273\n"

end # @testset "Cronbach's Alpha"


text:
function CartGrid(nx,ny,lx,ly,x0,y0,OrientFace,Boundary,Param)
Grid = GridStruct(;)
Grid.nBar=[ 0  1   0   1
           -1  0  -1   0];
Grid.nBar3=[ 0  1   0   1
           -1  0  -1   0
           0   0  0    0];
Grid.Dim=3;
Grid.Type="Quad";
Grid.Form="Planar";
Pert=0.0;
PertX=0.2;
PertY=0.2;
if strcmp(Boundary.WE,"Period") && strcmp(Boundary.BT,"Period")
  NumNodes=nx*ny;
elseif strcmp(Boundary.WE,"Period")
  NumNodes=nx*(ny+1);
elseif strcmp(Boundary.BT,"Period")
  NumNodes=(nx+1)*ny;
else
  NumNodes=(nx+1)*(ny+1);
end
Nodes = map(1:NumNodes) do i
  Node([0 0 0],0)
end
# Nodes[1:NumNodes]=Node([0 0 0],0);
NodeNumber=1;
dx=lx/nx;
dy=ly/ny;
y=y0;
P=zeros(3,nx+1,ny+1);
for iy=1:ny+1
  eta=(iy-1)/ny;
  x=x0;
  for ix=1:nx+1
    P[1,ix,iy]=x;
    P[2,ix,iy]=y;
    P[3,ix,iy]=0;
    x=x+dx;
  end
  y=y+dy;
end

y=y0;
for iy=1:ny+1
  x=x0;
  if iy==ny+1 && strcmp(Boundary.BT,"Period")
#   for ix=1:nx+1
#     if ix==nx+1 && strcmp(Boundary.WE,"Period")
#     else
#       Nodes[NodeNumber]=Node([x;y;0],NodeNumber);
#       NodeNumber=NodeNumber+1;
#     end
#     x=x+dx;
#   end
  else
    for ix=1:nx+1
      if ix==nx+1 && strcmp(Boundary.WE,"Period")
      else
        Nodes[NodeNumber]=Node([x;y;0],NodeNumber);
        NodeNumber=NodeNumber+1;
      end
      x=x+dx;
    end
  end
  y=y+dy;
end

Grid.Nodes=Nodes;
Grid.NumNodes=NumNodes;

if strcmp(Boundary.WE,"Period" ) && strcmp(Boundary.BT,"Period")
  NumEdges=2*nx*ny;
  NumEdgesX=nx*ny;
elseif strcmp(Boundary.WE,"Period")
  NumEdges=nx*ny+nx*(ny+1);
  NumEdgesX=nx*(ny+1);
else
  NumEdges=(nx+1)*ny+nx*(ny+1);
  #NumEdgesX=nx*ny;
end

Edges= map(1:NumEdges) do i
  Edge([1 2],Grid,0,0,"",0)
end
# Edges[1:NumEdges]=Edge([1 2],Grid,0,0,"",0);

EdgeNumber=1;
EdgeNumberX=1;
EdgeNumberY=1;
BC="";
if strcmp(Boundary.WE,"Period")
  N1=1;
  N2=nx+1;
else
  N1=1;
  N2=nx+2;
end
for iy=1:ny
  for ix=1:nx+1
    if ix==nx+1 && strcmp(Boundary.WE,"Period")
    else
      if iy==ny && strcmp(Boundary.BT,"Period")
        Edges[EdgeNumber]=Edge([N1 1+(ix-1)],Grid,EdgeNumber,EdgeNumber,"Y",EdgeNumberY);
        EdgeNumber=EdgeNumber+1;
        EdgeNumberY=EdgeNumberY+1;
        N1=N1+1;
        N2=N2+1;
      else
        Edges[EdgeNumber]=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber,"Y",EdgeNumberY);
        EdgeNumber=EdgeNumber+1;
        EdgeNumberY=EdgeNumberY+1;
        N1=N1+1;
        N2=N2+1;
      end
    end
  end
end
N1=1;
N2=2;
for iy=1:ny+1
  if iy==ny+1 && strcmp(Boundary.BT,"Period")
  else
    for ix=1:nx
      if ix==nx && strcmp(Boundary.WE,"Period")
        Edges[EdgeNumber]=Edge([N1 1+(iy-1)*nx],Grid,EdgeNumber,EdgeNumber,"X",EdgeNumberX);
        EdgeNumber=EdgeNumber+1;
        EdgeNumberX=EdgeNumberX+1;
        N1=N1+1;
        N2=N2+1;
      else
        Edges[EdgeNumber]=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber,"X",EdgeNumberX);
        EdgeNumber=EdgeNumber+1;
        EdgeNumberX=EdgeNumberX+1;
        N1=N1+1;
        N2=N2+1;
        if ix==nx
          N1=N1+1;
          N2=N2+1;
        end
      end
    end
  end
end

Grid.Edges=Edges;
Grid.NumEdges=NumEdges;

NumFaces=nx*ny;
Faces=map(1:NumFaces) do i
  (F, Grid) = Face([0 0 0 0],Grid,0,"x",[]);
  F
end
# Faces[1:NumFaces]=Face([0 0 0 0],Grid,0,"x",[]);
if strcmp(Boundary.WE,"Period")
  E1=nx*ny+1;
  E3=nx*ny+1+nx;
else
  E1=(nx+1)*ny+1;
  E3=(nx+1)*ny+1+nx;
end
E2=2;
E4=1;
FaceNumber=1;
Type="o";
for iy=1:ny
  if iy==ny && strcmp(Boundary.BT,"Period")
    for ix=1:nx
      if ix==nx && strcmp(Boundary.WE,"Period")
        (Faces[FaceNumber],Grid)=Face([E1 1+(iy-1)*nx NumEdgesX+1+(ix-1) E4],Grid,FaceNumber,Type,OrientFace,
          [P[:,ix,iy] P[:,ix+1,iy] P[:,ix+1,iy+1] P[:,ix,iy+1]]);
        FaceNumber=FaceNumber+1;
        E1=E1+1;
        E4=E4+1;
      else
        (Faces[FaceNumber],Grid)=Face([E1 E2 NumEdgesX+1+(ix-1) E4],Grid,FaceNumber,Type,OrientFace,
          [P[:,ix,iy] P[:,ix+1,iy] P[:,ix+1,iy+1] P[:,ix,iy+1]]);
        FaceNumber=FaceNumber+1;
        E1=E1+1;
        E2=E2+1;
        E4=E4+1;
      end
    end
  else
    for ix=1:nx
      if ix==nx && strcmp(Boundary.WE,"Period")
        (Faces[FaceNumber],Grid)=Face([E1 1+(iy-1)*nx E3 E4],Grid,FaceNumber,Type,OrientFace,
          [P[:,ix,iy] P[:,ix+1,iy] P[:,ix+1,iy+1] P[:,ix,iy+1]]);
        FaceNumber=FaceNumber+1;
        E1=E1+1;
        E3=E3+1;
      else
        (Faces[FaceNumber],Grid)=Face([E1 E2 E3 E4],Grid,FaceNumber,Type,OrientFace,
          [P[:,ix,iy] P[:,ix+1,iy] P[:,ix+1,iy+1] P[:,ix,iy+1]]);
        FaceNumber=FaceNumber+1;
        E1=E1+1;
        E2=E2+1;
        E3=E3+1;
        E4=E4+1;
      end
    end
    E2=E2+1;
    E4=E4+1;
  end
end
Grid.Faces=Faces;
Grid.NumFaces=NumFaces;
Grid=Orientation(Grid);
Grid=Renumbering(Grid);
Grid=FacesInNodes(Grid);
return Grid
end


text:
module ReactiveMPPermutationMatrixTest

using Test
using ReactiveMP
using ReactiveMP: getind, PT_X_P
using Random

using LinearAlgebra

@testset "Permutation Matrix" begin

    @testset "Creation" begin

        ind = shuffle(collect(1:100))
        P = PermutationMatrix(ind)
        @test P.ind == ind

        for k = 2:10
            P = PermutationMatrix(k)
            @test length(P.ind) == k
        end

        # check whether the first entry is one
        for k = 1:100
            P = PermutationMatrix(3)
            @test P[1,1] == 0
        end
        ind = 0
        for k = 1:1000
            P = PermutationMatrix(3; switch_first=false)
            if P[1,1] == 1
                ind += 1
            end
        end
        @test ind > 0

    end

    @testset "Get functions" begin

        ind = shuffle(collect(1:100))
        P = PermutationMatrix(ind)
        @test P.ind     == getind(P)

        ind = shuffle(collect(1:100))
        P = PermutationMatrix(ind)'
        @test sortperm(ind) == getind(P)

    end

    @testset "Simple base functions" begin

        P = PermutationMatrix(10)
        @test eltype(P) == Int64
        @test length(P) == 100
        @test size(P)   == (10, 10)

        @test sum(P)    == 10
        for k = 1:10
            @test sum(P[:,k]) == 1
            @test sum(P[k,:]) == 1
        end

        @test inv(P) == P'

    end

    @testset "Multiplication" begin

        # generate matrices
        P = PermutationMatrix(50)
        P_dense = collect(P)
        X = randn(50, 50)
        x = randn(50)

        @test P*x             == P_dense*x
        @test P'*x            == P_dense'*x
        @test transpose(P)*x  == transpose(P_dense)*x

        @test P*X             == P_dense*X
        @test P'*X            == P_dense'*X
        @test transpose(P)*X  == transpose(P_dense)*X

        @test P*X'            == P_dense*X'
        @test P'*X'           == P_dense'*X'
        @test transpose(P)*X' == transpose(P_dense)*X'

        @test P*transpose(X)             == P_dense*transpose(X)
        @test P'*transpose(X)            == P_dense'*transpose(X)
        @test transpose(P)*transpose(X)  == transpose(P_dense)*transpose(X)

        @test x'*P             == x'*P_dense
        @test x'*P'            == x'*P_dense'
        @test x'*transpose(P)  == x'*transpose(P_dense)

        @test transpose(x)*P             == transpose(x)*P_dense
        @test transpose(x)*P'            == transpose(x)*P_dense'
        @test transpose(x)*transpose(P)  == transpose(x)*transpose(P_dense)

        @test P*X             == P_dense*X
        @test P'*X            == P_dense'*X
        @test transpose(P)*X  == transpose(P_dense)*X

        @test P*X'            == P_dense*X'
        @test P'*X'           == P_dense'*X'
        @test transpose(P)*X' == transpose(P_dense)*X'

        @test P*transpose(X)             == P_dense*transpose(X)
        @test P'*transpose(X)            == P_dense'*transpose(X)
        @test transpose(P)*transpose(X)  == transpose(P_dense)*transpose(X)

        @test X*P             == X*P_dense
        @test X*P'            == X*P_dense'
        @test X*transpose(P)  == X*transpose(P_dense)

        @test X'*P            == X'*P_dense
        @test X'*P'           == X'*P_dense'
        @test X'*transpose(P) == X'*transpose(P_dense)

        @test transpose(X)*P             == transpose(X)*P_dense
        @test transpose(X)*P'            == transpose(X)*P_dense'
        @test transpose(X)*transpose(P)  == transpose(X)*transpose(P_dense)

        @test PT_X_P(X, P)    == P'*X*P
        @test PT_X_P(X, P)    == P_dense'*X*P_dense

    end

end
end


text:
using LabView0mqJl

try
    include(get_script_path(raw"C:\_LabView_projects\ZMQ\LabView0mqJl.jl\src\LabVIEW\jl-scripts\simple_examples.jl"))
    LabView0mqJl.setglobals(isOK=true)
catch excep
    LabView0mqJl.setglobals(isOK=false, excpn=(;excep, stack_trace=backtrace()))
end

server_0mq4lv(fns)


text:
using ..Helper: init

function simulate!(x, r, epsilon)
    dim, steps = size(x)
    for t = 2:steps

        # common = epsilon * mean(x[:, t-1])
        common = 0.0
        @simd for i = 1:dim
            common += x[i, t-1]
        end
        common *= epsilon / dim

        for i = 1:dim
            y = (1 - epsilon) * x[i, t-1] + common
            x[i, t] = r * y * (1 - y)
        end
    end
end

function prepare(name)
    res = init(name)
    function run()
        simulate!(res[:x], res[:r], res[:epsilon])
        res
    end
end

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
using BenchmarkTools

H1 = rand(ComplexF64,4,4); H1 = Hermitian(H1+H1');
H2 = rand(Float64,4,4);    H2 = Hermitian(H2+H2');

H = H2

w, Z = eig(H)

w2, Z2 = hermfact!(similar(w), similar(Z), copy(H))

@assert w==w2 && Z==Z2

a = @benchmark eig(H) setup=(A=copy($H))
println("BASE eigendecomposition")
show(STDOUT,MIME"text/plain"(),a)
println("\n")
b = @benchmark hermfact!(w2, Z2, A) setup=(A=copy($H))
println("In-place eigendecomposition")
show(STDOUT,MIME"text/plain"(),b)
println()


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ fb086364-68fc-11eb-32f6-ab7218186246
begin
	# using Pkg; Pkg.add("Distributions")
	using JuliaDB, Distributions
	using OnlineStats
	using PlutoUI
end

# ╔═╡ aeb3b894-68fa-11eb-18c7-516e32b88d55
md"""
## JuliaDB - tutorial

- ref. [Using table](https://www.youtube.com/watch?v=pv5zfIs2lyU)
- ref. [Selections in JuliaDB ](https://www.youtube.com/watch?v=eQORf4m_8Hk)
- ref. [Reductions in JuliaDB](https://www.youtube.com/watch?v=tCVgn9m2ajM)
- ref. [Grouping in JuliaDB to calculate descriptive statistics](https://www.youtube.com/watch?v=hU7N-EWnC2A)
- ref. [Joining 2 tables](https://www.youtube.com/watch?v=CIBxY7PUjc0)
- ref. [Importing a csv file in JuliaDB] (https://www.youtube.com/watch?v=slL_oj5Dmzc)
- ref. [JuliaDB](https://juliacomputing.github.io/JuliaDB)

"""

# ╔═╡ e0cc6636-690d-11eb-0885-9790a38923c8
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 21d810d6-68fb-11eb-3542-5906971c0c1a
md"""
### The `table()` function
"""

# ╔═╡ 53c79c1a-68fb-11eb-26a2-bf3828728dbb
begin
	ids = collect(1:10);
	ages = round.(Int8, rand(Normal(50, 20), 10)); # Normal distribution with μ=50 and σ^2 = 20
	length(ids), length(ages)
end

# ╔═╡ 715cd434-68fb-11eb-177f-eb8052e0787c
typeof(ids), typeof(ages)

# ╔═╡ 762ef460-68fb-11eb-1439-e7c2b6528c28
table_1 = table(ids, ages, names=[:ID, :AGE])

# ╔═╡ 953e99dc-68fb-11eb-30b6-5f9c7cc86257
subject_3 = table_1[3]

# ╔═╡ 9c2f284c-68fb-11eb-256c-ab5ee6c17d67
table_1[3:5]

# ╔═╡ a2e4491a-68fb-11eb-384a-652b7eb5ac71
# Another table (in_memory)

table_2 = table(ids,
    [:I, :I, :II, :III, :III, :II, :I, :III, :II, :I],
    round.(rand(Normal(16, 4), 10), digits=1), ## Normal dist μ=16 and σ^2 = 4
    round.(rand(Normal(12, 3), 10), digits=1), ## Normal dist. with μ=12 and σ^2 = 3
    names = [:ID, :GROUP, :HB, :WCC], pkey=:ID)

# ╔═╡ a3bd8a22-68fb-11eb-1333-7199c93b565d
table_3 = table(Columns(ID=ids, CRP=round.(Int16, rand(Normal(100, 20), 10))), pkey=:ID)

# ╔═╡ 09d1097e-68fc-11eb-1259-5da24a55ac7c
md"""
### Selections
"""

# ╔═╡ 2a824488-68fc-11eb-0de7-fbd2904440fb
## create new (sub-)table, using filter function passing a lambda and a table
## select all rowqs fro group :I
grouped_by_I = filter(t -> t.GROUP == :I, table_2)

# ╔═╡ 47bfa588-68fc-11eb-0365-a95c30fd3d39
grouped_by_III = filter(t -> t.GROUP == :III, table_2)

# ╔═╡ 533225c6-68fc-11eb-0d69-45914f9418db
high_WCC = filter(t -> t.WCC > 12, table_2)

# ╔═╡ 5eb97ade-68fc-11eb-18b4-d1a014439da9
md"""
### Reductions
"""

# ╔═╡ 758a03b4-68fc-11eb-3834-1710458d0191
table_2

# ╔═╡ 7f48e442-68fc-11eb-3526-45f76696ff14
## sum on column :WCC
reduce(+, table_2, select=:WCC; init=0)

# ╔═╡ 8d170fb8-68fc-11eb-197d-5fb88f9b12b2
## min and max values on column :WCC
## NOTE: we pass a tuple to reduce

reduce((min, max), table_2, select=:WCC)

# ╔═╡ a169a92e-68fc-11eb-23cd-a9c983423808
reduce((Mean(), Variance()), table_2,
	select=:WCC) # using function from OnlineStats package

# ╔═╡ 443cf75c-68fd-11eb-0ad0-25884c49c9ae
md"""
### Grouping
"""

# ╔═╡ 630e56b2-68fd-11eb-1bf5-2f69e0dd269a
## Calculate the mean of HB per group
groupreduce(Mean(), table_2, :GROUP, select = :HB)

# ╔═╡ 6b59c392-68fd-11eb-18dd-29050ae8cb39
## Calculate the variance of HB per group
groupreduce(Variance(), table_2, :GROUP, select = :HB)

# ╔═╡ 884358c4-68fd-11eb-3e8f-e11e24f1dc84
## grouping by mean, median, ...

groupby((mean, median, std, var, quantile), table_2, :GROUP, select=:HB)

# ╔═╡ b5e4fb40-68fd-11eb-3568-6754c8dcf55b
md"""
### Joining 2 tables
"""

# ╔═╡ d9c87b40-68fd-11eb-32ac-9108bff7f829
left_table = table(ids, ages, names=[:ID, :AGE], pkey=:ID)

# ╔═╡ f4fc8dd2-68fd-11eb-2309-ffe7574dc250
right_table = table(collect(1:15),
	round.(Int16, rand(Normal(100, 20), 15)),
	names=[:ID, :GROUP],
	pkey = :ID)

# ╔═╡ 14de3128-68fe-11eb-0c08-f5409401513e
join(left_table, right_table)  ## inner-join the 2 tables on ID

# ╔═╡ 28eedfd2-68fe-11eb-1b84-6dcf0ce6857b
join(right_table, left_table)  ## inner-join the 2 tables on ID, the other way...

# ╔═╡ 4d6808b6-68fe-11eb-0790-f1619d25e329
join(left_table, right_table, how=:left) ## inner-join the 2 tables on ID, from left table

# ╔═╡ 58c352f4-68fe-11eb-384d-d576e23d1d46
join(left_table, right_table, how=:outer) ## outer-join left and right table on ID

# ╔═╡ 648943d4-68fe-11eb-0b23-c71f63be5957
join(right_table, left_table, how=:anti)  ## only tuple that have no id in right table

# ╔═╡ b61dc290-68fe-11eb-1d40-e72940ae6be0
md"""
### Loading a CSV file
"""

# ╔═╡ c83fda3c-68fe-11eb-3112-3d01bf2d0693
## csv file is in the current directory
db = loadtable("data/Data.csv", header_exists=true)

# ╔═╡ 9e8aa658-68ff-11eb-1a80-f50efa89f446
length(db)

# ╔═╡ da1b2fee-690e-11eb-2ca1-1785bbccb4e5
with_terminal() do
	println(rows(db, (:CRIM, :ZN)))
end

# ╔═╡ 882f3840-6901-11eb-2329-f3ab3b4c22a1
begin
	n_t = filter(t -> t.CRIM > 2., db)

	with_terminal() do
		println(select(n_t, (:CRIM, :ZN, :AGE, :B)))
	end
end

# ╔═╡ 72b6b062-690e-11eb-1abc-d9ddf9f363e1
summarize((m = mean, s = std), db)

# ╔═╡ 89332e45-974b-4814-9bd5-b34b98eeb9d7
html"""
<style>
  main {
        max-width: calc(800px + 25px + 6px);
  }
  .plutoui-toc.aside {
    background: linen;
  }
  h3, h4 {
        background: wheat;
  }
</style>
"""

# ╔═╡ Cell order:
# ╟─aeb3b894-68fa-11eb-18c7-516e32b88d55
# ╠═fb086364-68fc-11eb-32f6-ab7218186246
# ╟─e0cc6636-690d-11eb-0885-9790a38923c8
# ╟─21d810d6-68fb-11eb-3542-5906971c0c1a
# ╠═53c79c1a-68fb-11eb-26a2-bf3828728dbb
# ╠═715cd434-68fb-11eb-177f-eb8052e0787c
# ╠═762ef460-68fb-11eb-1439-e7c2b6528c28
# ╠═953e99dc-68fb-11eb-30b6-5f9c7cc86257
# ╠═9c2f284c-68fb-11eb-256c-ab5ee6c17d67
# ╠═a2e4491a-68fb-11eb-384a-652b7eb5ac71
# ╠═a3bd8a22-68fb-11eb-1333-7199c93b565d
# ╟─09d1097e-68fc-11eb-1259-5da24a55ac7c
# ╠═2a824488-68fc-11eb-0de7-fbd2904440fb
# ╠═47bfa588-68fc-11eb-0365-a95c30fd3d39
# ╠═533225c6-68fc-11eb-0d69-45914f9418db
# ╟─5eb97ade-68fc-11eb-18b4-d1a014439da9
# ╠═758a03b4-68fc-11eb-3834-1710458d0191
# ╠═7f48e442-68fc-11eb-3526-45f76696ff14
# ╠═8d170fb8-68fc-11eb-197d-5fb88f9b12b2
# ╠═a169a92e-68fc-11eb-23cd-a9c983423808
# ╟─443cf75c-68fd-11eb-0ad0-25884c49c9ae
# ╠═630e56b2-68fd-11eb-1bf5-2f69e0dd269a
# ╠═6b59c392-68fd-11eb-18dd-29050ae8cb39
# ╠═884358c4-68fd-11eb-3e8f-e11e24f1dc84
# ╟─b5e4fb40-68fd-11eb-3568-6754c8dcf55b
# ╠═d9c87b40-68fd-11eb-32ac-9108bff7f829
# ╠═f4fc8dd2-68fd-11eb-2309-ffe7574dc250
# ╠═14de3128-68fe-11eb-0c08-f5409401513e
# ╠═28eedfd2-68fe-11eb-1b84-6dcf0ce6857b
# ╠═4d6808b6-68fe-11eb-0790-f1619d25e329
# ╠═58c352f4-68fe-11eb-384d-d576e23d1d46
# ╠═648943d4-68fe-11eb-0b23-c71f63be5957
# ╟─b61dc290-68fe-11eb-1d40-e72940ae6be0
# ╠═c83fda3c-68fe-11eb-3112-3d01bf2d0693
# ╠═9e8aa658-68ff-11eb-1a80-f50efa89f446
# ╠═da1b2fee-690e-11eb-2ca1-1785bbccb4e5
# ╠═882f3840-6901-11eb-2329-f3ab3b4c22a1
# ╠═72b6b062-690e-11eb-1abc-d9ddf9f363e1
# ╟─89332e45-974b-4814-9bd5-b34b98eeb9d7

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
module Tests

include("../src/Tokenizer.jl")
include("../src/Markovify.jl")
using .Tokenizer, .Markovify, Test

#=

Tests for the module Tokenizer

=#

text = "ABC.DAB.\n ACDA"

@test tokenize(text; on=letters) == [
    ["A", "B", "C", "."],
    ["D", "A", "B", "."],
    ["A", "C", "D", "A"],
]

@test tokenize(text; on=lines) == [
    ["A", "B", "C", ".", "D", "A", "B", "."],
    [" ", "A", "C", "D", "A"]
]

#=

Tests for the module Markovify

=#

tokens = tokenize(text; on=letters)
model = Model(tokens; order=1)

@test model.nodes ==
    Dict(
        ["A"] =>    Dict("B" => 2, "C" => 1, :end => 1),
        ["B"] =>    Dict("C" => 1, "." => 1),
        ["C"] =>    Dict("D" => 1,"." => 1),
        ["D"] =>    Dict("A" => 2),
        [:begin] => Dict("A" => 2,"D" => 1),
        ["."] =>    Dict(:end => 2)
    )

@test walk(model, ["A"])[1] == "A"

@test_throws KeyError walk(model, ["E"])

model2 = Model(tokens; order=2)

@test Markovify.states_with_suffix(model2, ["."]) == [
    ["B", "."],
    ["C", "."]
]

@test Markovify.state_with_prefix(model2, ["B"]) in [
    ["B", "."],
    ["B", "C"]
]

@test Markovify.indexin(collect(-1:12), 10) == 12
end


text:
module Comlineoption_module
    struct comlineoption_param
        inpname::String
        usethread::Bool
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants)

Multicomponent version of the compressible Euler equations
```math
\partial t
\begin{pmatrix}
\rho v_1 \\ \rho v_2 \\ E \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n}
\end{pmatrix}
+
\partial x
\begin{pmatrix}
\rho v_1^2 + p \\ \rho v_1 v_2 \\ (E+p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1
\end{pmatrix}
+
\partial y
\begin{pmatrix}
\rho v_1 v_2 \\ \rho v_2^2 + p \\ (E+p) v_2 \\ \rho_1 v_2 \\ \rho_2 v_2 \\ \vdots \\ \rho_{n} v_2
\end{pmatrix}
=
\begin{pmatrix}
0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0
\end{pmatrix}
```
for calorically perfect gas in two space dimensions.

In case of more than one component, the specific heat ratios `gammas` and the gas constants
`gas_constants` in [kJ/(kg*K)] should be passed as tuples, e.g., `gammas=(1.4, 1.667)`.

The remaining variables like the specific heats at constant volume 'cv' or the specific heats at
constant pressure 'cp' are then calculated considering a calorically perfect gas.
"""
struct CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT<:Real} <: AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP}
  gammas            ::SVector{NCOMP, RealT}
  gas_constants     ::SVector{NCOMP, RealT}
  cv                ::SVector{NCOMP, RealT}
  cp                ::SVector{NCOMP, RealT}

  function CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}(gammas       ::SVector{NCOMP, RealT},
                                                                           gas_constants::SVector{NCOMP, RealT}) where {NVARS, NCOMP, RealT<:Real}

    NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `gas_constants` have to be filled with at least one value"))

    cv = gas_constants ./ (gammas .- 1)
    cp = gas_constants + gas_constants ./ (gammas .- 1)

    new(gammas, gas_constants,cv, cp)
  end
end

function CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants)

  _gammas        = promote(gammas...)
  _gas_constants = promote(gas_constants...)
  RealT          = promote_type(eltype(_gammas), eltype(_gas_constants))

  NVARS = length(_gammas) + 3
  NCOMP = length(_gammas)

  __gammas        = SVector(map(RealT, _gammas))
  __gas_constants = SVector(map(RealT, _gas_constants))

  return CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}(__gammas, __gas_constants)
end

@inline Base.real(::CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}) where {NVARS, NCOMP, RealT} = RealT

function varnames(::typeof(cons2cons), equations::CompressibleEulerMulticomponentEquations2D)

  cons  = ("rho_v1", "rho_v2", "rho_e")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (cons..., rhos...)
end

function varnames(::typeof(cons2prim), equations::CompressibleEulerMulticomponentEquations2D)

  prim  = ("v1", "v2", "p")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (prim..., rhos...)
end

# Set initial conditions at physical location `x` for time `t`

"""
    initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f
  ini     = c + A * sin(omega * (x[1] + x[2] - t))

  v1      = 1.0
  v2      = 1.0

  rho     = ini

  # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1)
  prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))

  prim1 = rho * v1
  prim2 = rho * v2
  prim3 = rho^2

  prim_other = SVector{3, real(equations)}(prim1, prim2, prim3)

  return vcat(prim_other, prim_rho)
end

"""
    source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
@inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # Same settings as in `initial_condition`
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f

  gamma  = totalgamma(u, equations)

  x1, x2  = x
  si, co  = sincos((x1 + x2 - t)*omega)
  tmp1    = co * A * omega
  tmp2    = si * A
  tmp3    = gamma - 1
  tmp4    = (2*c - 1)*tmp3
  tmp5    = (2*tmp2*gamma - 2*tmp2 + tmp4 + 1)*tmp1
  tmp6    = tmp2 + c

  # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1
  du_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * tmp1 for i in eachcomponent(equations))

  du1 = tmp5
  du2 = tmp5
  du3 = 2*((tmp6 - 1.0)*tmp3 + tmp6*gamma)*tmp1

  du_other  = SVector{3, real(equations)}(du1, du2, du3)

  return vcat(du_other, du_rho)
end

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})

A shock-bubble testcase for multicomponent Euler equations
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  @unpack gas_constants = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.03

  # Region I
  rho1_1  = delta
  rho2_1  = 1.225 * gas_constants[1]/gas_constants[2] - delta
  v1_1    = zero(delta)
  v2_1    = zero(delta)
  p_1     = 101325

  # Region II
  rho1_2  = 1.225-delta
  rho2_2  = delta
  v1_2    = zero(delta)
  v2_2    = zero(delta)
  p_2     = 101325

  # Region III
  rho1_3  = 1.6861 - delta
  rho2_3  = delta
  v1_3    = -113.5243
  v2_3    = zero(delta)
  p_3     = 159060

  # Set up Region I & II:
  inicenter = SVector(zero(delta), zero(delta))
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  if (x[1] > 0.50)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r < 0.25)
    # Set up Region I
    rho1    = rho1_1
    rho2    = rho2_1
    v1      = v1_1
    v2      = v2_1
    p       = p_1
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
end

"""
    initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A for multicomponent adapted weak blast wave taken from
- Sebastian Hennemann, Gregor J. Gassner (2020)
  A provably entropy stable subcell shock capturing approach for high order split form DG
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter         = SVector(0.0, 0.0)
  x_norm            = x[1] - inicenter[1]
  y_norm            = x[2] - inicenter[2]
  r                 = sqrt(x_norm^2 + y_norm^2)
  phi               = atan(y_norm, x_norm)
  sin_phi, cos_phi  = sincos(phi)

  prim_rho          = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? 2^(i-1) * (1-2)/(1-2^ncomponents(equations))*1.0 : 2^(i-1) * (1-2)/(1-2^ncomponents(equations))*1.1691 for i in eachcomponent(equations))

  v1                = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2                = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  p                 = r > 0.5 ? 1.0 : 1.245

  prim_other         = SVector{3, real(equations)}(v1, v2, p)

  return prim2cons(vcat(prim_other, prim_rho),equations)
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e  = u

  rho = density(u, equations)

  v1    = rho_v1/rho
  v2    = rho_v2/rho
  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))

  if orientation == 1
    f_rho = densities(u, v1, equations)
    f1  = rho_v1 * v1 + p
    f2  = rho_v2 * v1
    f3  = (rho_e + p) * v1
  else
    f_rho = densities(u, v2, equations)
    f1  = rho_v1 * v2
    f2  = rho_v2 * v2 + p
    f3  = (rho_e + p) * v2
  end

  f_other  = SVector{3, real(equations)}(f1, f2, f3)

  return vcat(f_other, f_rho)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)

Adaption of the entropy conserving two-point flux by
- Ayoub Gouasmi, Karthik Duraisamy (2020)
  "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations""
  arXiv:1904.00972v3 [math.NA] 4 Feb 2020
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  # Unpack left and right state
  @unpack gammas, gas_constants, cv = equations
  rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
  rhok_mean   = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i+3], u_rr[i+3]) for i in eachcomponent(equations))
  rhok_avg    = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i+3] + u_rr[i+3]) for i in eachcomponent(equations))

  # Iterating over all partial densities
  rho_ll      = density(u_ll, equations)
  rho_rr      = density(u_rr, equations)

  # extract velocities
  v1_ll       = rho_v1_ll/rho_ll
  v2_ll       = rho_v2_ll/rho_ll
  v1_rr       = rho_v1_rr/rho_rr
  v2_rr       = rho_v2_rr/rho_rr
  v1_avg      = 0.5 * (v1_ll + v1_rr)
  v2_avg      = 0.5 * (v2_ll + v2_rr)
  v1_square   = 0.5 * (v1_ll^2 + v1_rr^2)
  v2_square   = 0.5 * (v2_ll^2 + v2_rr^2)
  v_sum       = v1_avg + v2_avg

  enth      = zero(v_sum)
  help1_ll  = zero(v1_ll)
  help1_rr  = zero(v1_rr)

  for i in eachcomponent(equations)
    enth      += rhok_avg[i] * gas_constants[i]
    help1_ll  += u_ll[i+3] * cv[i]
    help1_rr  += u_rr[i+3] * cv[i]
  end

  T_ll        = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll
  T_rr        = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr
  T           = 0.5 * (1.0/T_ll + 1.0/T_rr)
  T_log       = ln_mean(1.0/T_ll, 1.0/T_rr)

  # Calculate fluxes depending on orientation
  help1       = zero(T_ll)
  help2       = zero(T_rr)
  if orientation == 1
    f_rho       = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v1_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      help1     += f_rho[i] * cv[i]
      help2     += f_rho[i]
    end
    f1 = (help2) * v1_avg + enth/T
    f2 = (help2) * v2_avg
    f3 = (help1)/T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2
  else
    f_rho       = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v2_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      help1     += f_rho[i] * cv[i]
      help2     += f_rho[i]
    end
    f1 = (help2) * v1_avg
    f2 = (help2) * v2_avg + enth/T
    f3 = (help1)/T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2
  end
  f_other  = SVector{3, real(equations)}(f1, f2, f3)

  return vcat(f_other, f_rho)
end

"""
    flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
                 equations::CompressibleEulerMulticomponentEquations2D)

Adaption of the entropy conserving and kinetic energy preserving two-point flux by
- Hendrik Ranocha (2018)
  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
  for Hyperbolic Balance Laws
  [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
See also
- Hendrik Ranocha (2020)
  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
  the Euler Equations Using Summation-by-Parts Operators
  [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
"""
@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  # Unpack left and right state
  @unpack gammas, gas_constants, cv = equations
  rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
  rhok_mean   = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i+3], u_rr[i+3]) for i in eachcomponent(equations))
  rhok_avg    = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i+3] + u_rr[i+3]) for i in eachcomponent(equations))

  # Iterating over all partial densities
  rho_ll      = density(u_ll, equations)
  rho_rr      = density(u_rr, equations)

  # Calculating gamma
  gamma               = totalgamma(0.5*(u_ll+u_rr), equations)
  inv_gamma_minus_one = 1/(gamma-1)

  # extract velocities
  v1_ll               = rho_v1_ll / rho_ll
  v1_rr               = rho_v1_rr / rho_rr
  v1_avg              = 0.5 * (v1_ll + v1_rr)
  v2_ll               = rho_v2_ll / rho_ll
  v2_rr               = rho_v2_rr / rho_rr
  v2_avg              = 0.5 * (v2_ll + v2_rr)
  velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)

  # helpful variables
  help1_ll  = zero(v1_ll)
  help1_rr  = zero(v1_rr)
  enth_ll   = zero(v1_ll)
  enth_rr   = zero(v1_rr)
  for i in eachcomponent(equations)
    enth_ll   += u_ll[i+3] * gas_constants[i]
    enth_rr   += u_rr[i+3] * gas_constants[i]
    help1_ll  += u_ll[i+3] * cv[i]
    help1_rr  += u_rr[i+3] * cv[i]
  end

  # temperature and pressure
  T_ll            = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll
  T_rr            = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr
  p_ll            = T_ll * enth_ll
  p_rr            = T_rr * enth_rr
  p_avg           = 0.5 * (p_ll + p_rr)
  inv_rho_p_mean  = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)

  f_rho_sum = zero(T_rr)
  if orientation == 1
    f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v1_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      f_rho_sum += f_rho[i]
    end
    f1 = f_rho_sum * v1_avg + p_avg
    f2 = f_rho_sum * v2_avg
    f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)
  else
    f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v2_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      f_rho_sum += f_rho[i]
    end
    f1 = f_rho_sum * v1_avg
    f2 = f_rho_sum * v2_avg + p_avg
    f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + 0.5 * (p_ll*v2_rr + p_rr*v2_ll)
  end

  # momentum and energy flux
  f_other  = SVector{3, real(equations)}(f1, f2, f3)

  return vcat(f_other, f_rho)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # Get the density and gas gamma
  rho_ll   = density(u_ll, equations)
  rho_rr   = density(u_rr, equations)
  gamma_ll = totalgamma(u_ll, equations)
  gamma_rr = totalgamma(u_rr, equations)

  # Get the velocities based on direction
  if orientation == 1
    v_ll = rho_v1_ll / rho_ll
    v_rr = rho_v1_rr / rho_rr
  else # orientation == 2
    v_ll = rho_v2_ll / rho_ll
    v_rr = rho_v2_rr / rho_rr
  end

  # Compute the sound speeds on the left and right
  p_ll = (gamma_ll - 1) * (rho_e_ll - 1/2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll)
  c_ll = sqrt(gamma_ll * p_ll / rho_ll)
  p_rr = (gamma_rr - 1) * (rho_e_rr - 1/2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr)
  c_rr = sqrt(gamma_rr * p_rr / rho_rr)

  λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
end

@inline function max_abs_speeds(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  rho   = density(u, equations)
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho

  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  c     = sqrt(gamma * p / rho)

  return (abs(v1) + c, abs(v2) + c, )
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  prim_rho = SVector{ncomponents(equations), real(equations)}(u[i+3] for i in eachcomponent(equations))

  rho   = density(u, equations)
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho
  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))
  prim_other =  SVector{3, real(equations)}(v1, v2, p)

  return vcat(prim_other, prim_rho)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas, gas_constants = equations
  rho_v1, rho_v2, rho_e = u

  rho       = density(u, equations)

  # Multicomponent stuff
  help1         = zero(rho)
  gas_constant  = zero(rho)
  for i in eachcomponent(equations)
    help1         += u[i+3] * cv[i]
    gas_constant  += gas_constants[i] * (u[i+3]/rho)
  end

  v1        = rho_v1 / rho
  v2        = rho_v2 / rho
  v_square  = v1^2 + v2^2
  gamma     = totalgamma(u, equations)

  p         = (gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s         = log(p) - gamma * log(rho) - log(gas_constant)
  rho_p     = rho / p
  T         = (rho_e - 0.5 * rho * v_square) / (help1)
  entrop_rho  = SVector{ncomponents(equations), real(equations)}( gas_constant * ((gamma - s)/(gamma - 1.0) - (0.5 * v_square * rho_p)) for i in eachcomponent(equations))

  w1        = gas_constant * v1 * rho_p
  w2        = gas_constant * v2 * rho_p
  w3        = gas_constant * rho_p * (-1)

  entrop_other = SVector{3, real(equations)}(w1, w2, w3)

  return vcat(entrop_other, entrop_rho)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas = equations
  v1, v2, p = prim

  cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i+3] for i in eachcomponent(equations))
  rho     = density(prim, equations)
  gamma   = totalgamma(prim, equations)

  rho_v1  = rho * v1
  rho_v2  = rho * v2
  rho_e   = p/(gamma-1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2)

  cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e)

  return vcat(cons_other, cons_rho)
end

"""
    totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D)

Function that calculates the total gamma out of all partial gammas using the
partial density fractions as well as the partial specific heats at constant volume.
"""
@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas = equations

  help1 = zero(u[1])
  help2 = zero(u[1])

  for i in eachcomponent(equations)
    help1 += u[i+3] * cv[i] * gammas[i]
    help2 += u[i+3] * cv[i]
  end

  return help1/help2
end

@inline function density_pressure(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  rho          = density(u, equations)
  gamma        = totalgamma(u, equations)
  rho_times_p  = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))

  return rho_times_p
end

@inline function density(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho = zero(u[1])

  for i in eachcomponent(equations)
    rho += u[i+3]
  end

  return rho
 end

 @inline function densities(u, v, equations::CompressibleEulerMulticomponentEquations2D)

  return SVector{ncomponents(equations), real(equations)}(u[i+3]*v for i in eachcomponent(equations))
 end

end # @muladd


text:
# Unit tests for analysis/visualization.jl

@testset "analysis/visualization.jl" begin

end


text:
using Mimi

@defcomp malnutrition begin
	WDF = Variable() # Water depravation factor
    EF = Variable() # Effect factor
	HDF = Variable() # Human development factor related to vulnerability to malnutrition
    HHF = Variable() # human health factor

	WSI = Parameter() # Water stress index
	WU = Parameter() # Percentage of water consumption by agriculture
    WR = Parameter() # Minimum per capita water requirement of water for the agricultural sector to avoid malnutrition
	HDI = Parameter() # Human Development Index
    DF_malnutrition = Parameter() # Damage factor for malnutrition

	function run_timestep(p, v, d, t)
		v.WDF = p.WSI * p.WU
        if p.HDI < 0.3
			v.HDF = 1
		elseif p.HDI > 0.88
			v.HDF = 0
		else
			v.HDF = 2.03 * p.HDI^2 - 4.09 * p.HDI + 2.04
		end
        v.EF = p.WR^(-1) * v.HDF
        v.HHF = v.WDF * v.EF * p.DF_malnutrition
	end
end

@defcomp social_cost begin
    DALY_value = Variable() # DALY monetary value in the production site's country
    SC = Variable() # Social cost of water consumption

    DALY_ref = Parameter() # DALY reference value
    GNI = Parameter() # Gross national income per capita of the production site's country, adjusted for purchasing power parity
    GNI_ref = Parameter() # Gross national income per capita of the reference country
    e = Parameter() # Income elasticity of willingness to pay for health or life
    HHF = Parameter() # Human Health Factor

    function run_timestep(p, v, d, t)
        v.DALY_value = p.DALY_ref * (p.GNI / p.GNI_ref)^p.e
        v.SC = p.HHF * v.DALY_value
    end

end

function construct_model()
	# initiate the model
	m = Model()
    set_dimension!(m, :time, collect(1:1))

	# add components to the model
	add_comp!(m, malnutrition)
	add_comp!(m, social_cost)

	# give value to parameters
	update_param!(m, :malnutrition, :WSI, 0.4322) # for China, source Worldbank, https://data.worldbank.org/indicator/ER.H2O.FWST.ZS
	update_param!(m, :malnutrition, :WU, 0.64) # https://data.worldbank.org/indicator/ER.H2O.FWAG.ZS
	update_param!(m, :malnutrition, :WR, 1350.) # m3, per capita per year, from Pfister 2009
	update_param!(m, :malnutrition, :HDI, 0.761) #https://hdr.undp.org/en/countries/profiles/CHN
	update_param!(m, :malnutrition, :DF_malnutrition, 0.0184) #Pfister 2009
	update_param!(m, :social_cost, :DALY_ref, 185000) #VBA 2021
	update_param!(m, :social_cost, :GNI, 16201.4) #https://data.worldbank.org/indicator/NY.GNP.PCAP.PP.KD?locations=CN
	update_param!(m, :social_cost, :GNI_ref, 42135) # https://data.worldbank.org/indicator/NY.GNP.PCAP.PP.KD?locations=OE
	update_param!(m, :social_cost, :e, 0.6) # VBA 2021

    connect_param!(m, :social_cost, :HHF, :malnutrition, :HHF)

	return m
end

m = construct_model()

run(m)

m[:malnutrition, :HHF]

m[:social_cost, :DALY_value]

m[:social_cost, :SC]

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
function readInvestments(filename::String)

    df = readtable(filename, nastrings=[".", "", "NA"])

    # find columns that have been parsed as Strings by readtable
    col_to_test = Array(Symbol, 0)

    nCols = size(df, 2)
    for ii=1:nCols
        typeof(df[1, ii]) == UTF8String?
        push!(col_to_test, names(df)[ii]):
        nothing
    end

    # test each column's data to see if Datetime will parse it
    col_that_pass = Array(Symbol, 0)

    for colname in col_to_test
        d = match(r"[-|\s|\/|.]", df[1, colname])
        d !== nothing? (bar = split(df[1, colname], d.match)): (bar = [])
        if length(bar) == 3
            push!(col_that_pass, colname)
        end
    end

    # parse first column that passes the Datetime regex test
    idx = Date{ISOCalendar}[date(d) for d in
                                 df[col_that_pass[1]]] # without
                                        # Date it would fail chkIdx
                                        # in constructor

    delete!(df, [col_that_pass[1]])

    ## try whether DataFrame fits subtypes
    invs = AssetMgmt.Investments(df, idx)

    return invs
end

function writeInvestments(filename::String, invs::Investments)
    ## create large dataframe
    idxDf = DataFrame(idx = idx(invs));
    df = [idxDf invs.vals];
    writetable(filename, df)
end


text:
# using typeof built-in
# using command line
# julia> is the prompt

# simple integer division
julia> 7/7
1.0

# ans is builtin!
# similar to _ in Python
julia> typeof(ans)
Float64


text:
module Marble
# Markdown becoming Latex

export
    # Types
    MarbleEnv,
    # Document prep functions
    prepair, parse, process, render, template, build, changed, cache, save,
    # General use functions
    build_stream, build_file, build_dir, init_dir, print_settings

# External Packages
using Formatting
using FileState
using JinjaTemplates
using JSON
using SettingsBundles
using SHA
using YAML

# Internal Modules
# include("flavor/Marble_MD.jl")
# using MarbleFlavor

include("MarbleDoc.jl")
include("MarbleEnv.jl")
include("flavor/Marble_MD.jl") # Custom MD flavor
include("render.jl") # Wrapper around Jinja2
include("util.jl") # Utility functions
include("buildchain.jl") # Main application logic

end # module


text:
using Circuitscape, Documenter

makedocs(
    sitename = "Circuitscape.jl Documentation",
    pages = ["Home" => "index.md",
             "User Guide" => "usage.md",
             "Inputs, Outputs and Options" => "options.md",
             "On Solvers and Computation Time" => "compute.md",
             "Calling Circuitscape from other Programs" => "otherprog.md",
             "Logging Options" => "logging.md"
    ]
)

deploydocs(
    repo = "github.com/Circuitscape/Circuitscape.jl.git",
    devbranch = "master",
    devurl = "latest"
)


text:
using SpecialFunctions
"""
    errorf()
``errorf'' is an NLP instance taken from
https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nbsir85-3206.pdf
The best known primal feasible objective value for this problem approximately 0.
This example contains 5 nonlinear terms.
Variables: 0 <= x[1:2] <= 10;
Objective: Minimize errorf(x[1] + x[2]) + sin(x[1]) * exp(-x[2]/10)
Constraints: x[1]^2 + x[2]^2 >= 10

The reformulation using multiple additional variables is as follows:
Variables:
    0 <= x[1:2] <= 10
    x[3:4] - transformation variables
    -0.632121 <= x[3] <= 1.90484
    -1.90484 <= x[4] <= 0.632121
    y[1:8]
    0 <= y[1] <= 20
Constraints:
    y[1] = x[1] + x[2]
    y[2] = errorf(y[1])
    x[3] = y[3] + y[4]
    x[4] = y[3] - y[4]
    y[3] = sin(x[1])
    y[4] = exp(-x[2]/10)
    y[5] = x[1]^2
    y[6] = x[2]^2
    y[7] = x[3]^2
    y[8] = x[4]^2
Objective: Minimize y[2] + 0.25 * (y[7] - y[8])
"""

function errorf()
    best_known_objective = round(erf(3*π/2) + sin(3*π/2) * 1.0; digits=4)
    cplex_optimizer =
        JuMP.optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0)
    error_tolerances = [NaN64, 1e-3]
    # base partition holds the inflection points and the end-points
    milp_statistics = zeros(length(error_tolerances), 5)
    lp_statistics = zeros(length(error_tolerances), 5)

    base_partition = Dict{Int,Vector{Float64}}()
    base_partition[1] = [0.0, 1.0] # erf()
    base_partition[2] = [0.0, π, 2*π, 3*π, 10] # sin()
    base_partition[3] = [0, 10.0] # exp
    base_partition[4] = [0, 10.0] # x[1]^2
    base_partition[5] = [0, 10.0] # x[2]^2
    base_partition[6] = [-0.632121, 1.90484] # x[3]^2
    base_partition[7] = [-1.90484, 0.632121] # x[4]^2
    # functions
    functions = Dict{Int,Function}()
    functions[1] = x -> erf(x)
    functions[2] = x -> sin(x)
    functions[3] = x -> exp(-x/10)
    functions[4] = x -> x^2
    functions[5] = x -> x^2
    functions[6] = x -> x^2
    functions[7] = x -> x^2

        # solve MILP for various error tolerances
        for i = 1:length(error_tolerances)
            error_tolerance = error_tolerances[i]
            milp_var = Dict{Int,Any}()
            x_indexes = Dict{Int,Int}()
            y_indexes = Dict{Int,Int}()
            # construct MILP relaxations
            milp = Model(cplex_optimizer)
            x_lb = [0, 0, -0.632121, -1.90484]
            x_ub = [10, 10, 1.90484, 0.632121]
            @variable(milp, x_lb[i] <= x[i=1:4] <= x_ub[i])
            @variable(milp, y[1:8])
            set_lower_bound(y[1], 0)
            set_upper_bound(y[1], 20)
            @constraint(milp, y[1] == x[1] + x[2])
            @constraint(milp, x[3] == y[3] + y[4])
            @constraint(milp, x[4] == y[3] - y[4])
            @objective(milp, Min, y[2] + 0.25 * (y[7] - y[8]))
            var_map = Dict{Int,Any}()
            var_map[1] = (y[1], y[2]) # (indep, dep)
            var_map[2] = (x[1], y[3])
            var_map[3] = (x[2], y[4])
            var_map[4] = (x[1], y[5])
            var_map[5] = (x[2], y[6])
            var_map[6] = (x[3], y[7])
            var_map[7] = (x[4], y[8])
            for j = 1:7
                var, x_index, y_index = add_milp_relaxation(
                    functions[j],
                    base_partition[j],
                    error_tolerance,
                    milp,
                    name = string(j),
                )
                milp_var[j] = var
                x_indexes[j] = x_index
                y_indexes[j] = y_index
                @constraint(milp, var_map[j][1] == var[x_index])
                @constraint(milp, var_map[j][2] == var[y_index])
            end
            _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(milp)
            relaxation_objective = objective_value(milp)
            @assert relaxation_objective ≤ best_known_objective
            relative_gap =
                abs(best_known_objective - relaxation_objective) / (abs(relaxation_objective) + 1e-4)
            milp_statistics[i, :] = [
                error_tolerance,
                best_known_objective,
                relaxation_objective,
                relative_gap,
                round(solve_time, digits = 2),
            ]
        end

        for i = 1:length(error_tolerances)
            error_tolerance = error_tolerances[i]
            lp_var = Dict{Int,Any}()
            x_indexes = Dict{Int,Int}()
            y_indexes = Dict{Int,Int}()
            # construct LP relaxations
            lp = Model(cplex_optimizer)
            x_lb = [0, 0, -0.632121, -1.90484]
            x_ub = [10, 10, 1.90484, 0.632121]
            @variable(lp, x_lb[i] <= x[i=1:4] <= x_ub[i])
            @variable(lp, y[1:8])
            set_lower_bound(y[1], 0)
            set_upper_bound(y[1], 20)
            @constraint(lp, y[1] == x[1] + x[2])
            @constraint(lp, x[3] == y[3] + y[4])
            @constraint(lp, x[4] == y[3] - y[4])
            @objective(lp, Min, y[2] + 0.25 * (y[7] - y[8]))
            var_map = Dict{Int,Any}()
            var_map[1] = (y[1], y[2]) # (indep, dep)
            var_map[2] = (x[1], y[3])
            var_map[3] = (x[2], y[4])
            var_map[4] = (x[1], y[5])
            var_map[5] = (x[2], y[6])
            var_map[6] = (x[3], y[7])
            var_map[7] = (x[4], y[8])
            for j = 1:7
                var, x_index, y_index = add_lp_relaxation(
                    functions[j],
                    base_partition[j],
                    error_tolerance,
                    lp,
                    name = string(j),
                )
                lp_var[j] = var
                x_indexes[j] = x_index
                y_indexes[j] = y_index
                @constraint(lp, var_map[j][1] == var[x_index])
                @constraint(lp, var_map[j][2] == var[y_index])
            end
            _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(lp)
            relaxation_objective = objective_value(lp)
            @assert relaxation_objective ≤ best_known_objective
            relative_gap =
                abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
            lp_statistics[i, :] = [
                error_tolerance,
                best_known_objective,
                relaxation_objective,
                relative_gap,
                round(solve_time, digits = 2),
            ]
        end

        open("csv/erf_milp.csv", "w") do io
            writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
            writedlm(io, milp_statistics, ",")
        end

        open("csv/erf_lp.csv", "w") do io
            writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
            writedlm(io, lp_statistics, ",")
        end

end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo

# ## Mechanism
mech = get_mechanism(:tippetop,
    timestep=0.01,
    gravity=-9.81,
    friction_coefficient=0.4,
    contact=true,
    contact_type=:nonlinear)

# ## Simulate
initialize!(mech, :tippetop,
    x=[0.0, 0.0, 1.0],
    θ=[0,0.01,0],
    ω=[0.0, 0.01, 50.0])

storage = simulate!(mech, 10.0,
    record=true,
    verbose=false,
    opts=SolverOptions(verbose=false, btol=1e-5))

# ## Open visualizer
vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis)


text:
function CROSSOVER_binary_uniform(p1::Sol, p2::Sol)

    debug && DEBUG_feasible_solution(p1)
    debug && DEBUG_feasible_solution(p2)

    momkp = p1.prob
    n = momkp.n

    uniform_sample = rand(1:2, n)

    offspring = SOLUTION_empty_solution(momkp)

    for idx_var in 1:n
        # get the value of the gene
        value = p1.x[idx_var]
        if uniform_sample[idx_var] == 2
            value = p2.x[idx_var]
        end

        # if the item has been added, modify the solution
        if value == 1
            SOLUTION_add_item!(offspring, idx_var)
        end
    end

    debug && DEBUG_correct_solution(offspring)

    return offspring

end

function CROSSOVER_one_point(p1::Sol, p2::Sol)

    debug && DEBUG_feasible_solution(p1)
    debug && DEBUG_feasible_solution(p2)

    momkp = p1.prob
    n = momkp.n

    point = rand(1:n)

    offspring = SOLUTION_empty_solution(momkp)

    for idx_var in 1:n
        # get the value of the gene
        if idx_var <= point
            value = p1.x[idx_var]
        else
            value = p2.x[idx_var]
        end

        # if the item should be been added, modify the solution
        if value == 1
            SOLUTION_add_item!(offspring, idx_var)
        end
    end

    debug && DEBUG_correct_solution(offspring)

    return offspring

end

function CROSSOVER_two_points(p1::Sol, p2::Sol)

    debug && DEBUG_feasible_solution(p1)
    debug && DEBUG_feasible_solution(p2)

    momkp = p1.prob
    n = momkp.n

    point_1 = rand(1:n-1)
    point_2 = rand(point_1+1:n)

    offspring = SOLUTION_empty_solution(momkp)

    for idx_var in 1:n
        # get the value of the gene
        if idx_var <= point_1
            value = p1.x[idx_var]
        elseif idx_var <= point_2
            value = p2.x[idx_var]
        else
            value = p1.x[idx_var]
        end

        # if the item should be been added, modify the solution
        if value == 1
            SOLUTION_add_item!(offspring, idx_var)
        end
    end

    debug && DEBUG_correct_solution(offspring)

    return offspring

end


text:
struct BilinearForms
    cellmatrices::Any
    celltomatrix::Any
    ncells::Any
    function BilinearForms(cellmatrices, celltomatrix)
        nphase, ncells = size(celltomatrix)
        @assert nphase == 2
        @assert all(celltomatrix .>= 0)
        @assert all(celltomatrix .<= length(cellmatrices))
        new(cellmatrices, celltomatrix, ncells)
    end
end

function BilinearForms(basis, cellquads, stiffnesses, cellsign, cellmap)
    ncells = length(cellsign)

    uniformquad = uniform_cell_quadrature(cellquads)
    uniformbf1 = bilinear_form(basis, uniformquad, stiffnesses[+1], cellmap)
    uniformbf2 = bilinear_form(basis, uniformquad, stiffnesses[-1], cellmap)

    cellmatrices = [uniformbf1, uniformbf2]
    celltomatrix = zeros(Int, 2, ncells)

    for cellid = 1:ncells
        if cellsign[cellid] == +1
            celltomatrix[1, cellid] = 1
        elseif cellsign[cellid] == -1
            celltomatrix[2, cellid] = 2
        else
            pquad = cellquads[+1, cellid]
            pbf = bilinear_form(basis, pquad, stiffnesses[+1], cellmap)
            push!(cellmatrices, pbf)
            celltomatrix[1, cellid] = length(cellmatrices)

            nquad = cellquads[-1, cellid]
            nbf = bilinear_form(basis, nquad, stiffnesses[-1], cellmap)
            push!(cellmatrices, nbf)
            celltomatrix[2, cellid] = length(cellmatrices)
        end
    end
    return BilinearForms(cellmatrices, celltomatrix)
end

function BilinearForms(basis, cellquads, stiffnesses, cutmesh)
    cellsign = cell_sign(cutmesh)
    cellmap = cell_map(cutmesh, 1)
    return BilinearForms(basis, cellquads, stiffnesses, cellsign, cellmap)
end

function Base.getindex(cbf::BilinearForms, s, cellid)
    row = cell_sign_to_row(s)
    return cbf.cellmatrices[cbf.celltomatrix[row, cellid]]
end

function Base.getindex(cbf::BilinearForms, s)
    idx = cell_sign_to_row(s)
    return cbf.cellmatrices[idx]
end

function Base.show(io::IO, bf::BilinearForms)
    ncells = bf.ncells
    nuniquematrices = length(bf.cellmatrices)
    str = "BilinearForms\n\tNum. Cells: $ncells\n\tNum. Unique Cell Matrices: $nuniquematrices"
    print(io, str)
end


text:
# Run tests on the fiber problem in NLEVP (bessel function nonlinearity)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

push!(LOAD_PATH, string(@__DIR__, "/../src/gallery_extra"))
using GalleryNLEVP

@testset "NLEVP cd_player" begin
    nep_org = nep_gallery(NLEVP_NEP,"cd_player")
    nep_nat = nep_gallery("nlevp_native_cd_player")
    n = size(nep_org,1)

    @bench @testset "compare to native" begin
        for λ = [0.8236475079774124 + 0.9103565379264364im, 0.28604423575660154 + 0.5273506342660073im]
            for d = [0,1,2,3,4]
                M_org = compute_Mder(nep_org, λ, d)
                M_nat = compute_Mder(nep_nat, λ, d)
                @test opnorm(M_org - M_nat)/(opnorm(M_org)+eps()) < 1e-14
            end
        end

    end

end


text:
function freqresp(sys::DelayLtiSystem, ω::AbstractVector{T}) where {T <: Real}
    ny = noutputs(sys)
    nu = ninputs(sys)

    P_fr = freqresp(sys.P.P, ω);

    G_fr = zeros(eltype(P_fr), length(ω), ny, nu)

    for ω_idx=1:length(ω)
        P11_fr = P_fr[ω_idx, 1:ny, 1:nu]
        P12_fr = P_fr[ω_idx, 1:ny, nu+1:end]
        P21_fr = P_fr[ω_idx, ny+1:end, 1:nu]
        P22_fr = P_fr[ω_idx, ny+1:end, nu+1:end]

        delay_matrix_inv_fr = Diagonal(exp.(im*ω[ω_idx]*sys.Tau)) # Frequency response of the diagonal matrix with delays
        # Inverse of the delay matrix, so there should not be any minus signs in the exponents

        G_fr[ω_idx,:,:] .= P11_fr + P12_fr/(delay_matrix_inv_fr - P22_fr)*P21_fr # The matrix is invertible (?!)
    end

    return G_fr
end

function evalfr(sys::DelayLtiSystem, s)
    (ny, nu) = size(sys)

    P_fr = evalfr(sys.P.P, s)

    P11_fr = P_fr[1:ny, 1:nu]
    P12_fr = P_fr[1:ny, nu+1:end]
    P21_fr = P_fr[ny+1:end, 1:nu]
    P22_fr = P_fr[ny+1:end, nu+1:end]

    delay_matrix_inv_fr = Diagonal(exp.(s*sys.Tau))

    return P11_fr + P12_fr/(delay_matrix_inv_fr - P22_fr)*P21_fr
end

"""
    `y, t, x = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...)`

    Simulate system `sys`, over time `t`, using input signal `u`, with initial state `x0`, using method `alg` .

    Arguments:

    `t`: Has to be an `AbstractVector` with equidistant time samples (`t[i] - t[i-1]` constant)
    `u`: Function to determine control signal `ut` at a time `t`, on any of the following forms:
        Can be a constant `Number` or `Vector`, interpreted as `ut .= u` , or
        Function `ut .= u(t)`, or
        In-place function `u(ut, t)`. (Slightly more effienct)
    `kwargs...`: These are sent to `solve` from DelayDiffEq.

    Returns: times `t`, and `y` and `x` at those times.
"""
function lsim(sys::DelayLtiSystem{T,S}, u, t::AbstractArray{<:Real}; x0=fill(zero(T), nstates(sys)), alg=MethodOfSteps(Tsit5()), kwargs...) where {T,S}
    # Make u! in-place function of u
    u! = if isa(u, Number) || isa(u,AbstractVector) # Allow for u to be a constant number or vector
        (uout, t) -> uout .= u
    elseif DiffEqBase.isinplace(u, 2)               # If u is an inplace (more than 1 argument function)
        u
    else                                            # If u is a regular u(t) function
        (out, t) -> (out .= u(t))
    end
    _lsim(sys, u!, t, x0, alg; kwargs...)
end

function dde_param(dx, x, h!, p, t)
    A, B1, B2, C2, D21, Tau, u!, uout, hout, tmp = p

    u!(uout, t)     # uout = u(t)

    #dx .= A*x + B1*ut
    mul!(dx, A, x)
    mul!(tmp, B1, uout)
    dx .+= tmp

    @views for k=1:length(Tau)     # Add each of the delayed signals
        u!(uout, t-Tau[k])      # uout = u(t-tau[k])
        h!(hout, p, t-Tau[k])
        dk_delayed = dot(C2[k,:], hout) + dot(D21[k,:], uout)
        dx .+= B2[:, k] .* dk_delayed
    end
    return
end

# TODO Discontinuities in u are not handled well yet.
function _lsim(sys::DelayLtiSystem{T,S}, Base.@nospecialize(u!), t::AbstractArray{<:Real}, x0::Vector{T}, alg; kwargs...) where {T,S}
    P = sys.P

    if !iszero(P.D22)
        throw(ArgumentError("non-zero D22-matrix block is not supported")) # Due to limitations in differential equations
    end

    t0 = first(t)
    dt = t[2] - t[1]
    if !all(diff(t) .≈ dt) # QUESTION Does this work or are there precision problems?
        error("The t-vector should be uniformly spaced, t[2] - t[1] = $dt.") # Perhaps dedicated function for checking this?
    end

    # Get all matrices to save on allocations
    A, B1, B2, C1, C2, D11, D12, D21, D22 = P.A, P.B1, P.B2, P.C1, P.C2, P.D11, P.D12, P.D21, P.D22
    Tau = sys.Tau

    hout = fill(zero(T), nstates(sys))  # in place storage for h
    uout = fill(zero(T), ninputs(sys)) # in place storage for u
    tmp = similar(x0)

    h!(out, p, t) = (out .= 0)      # History function

    p = (A, B1, B2, C2, D21, Tau, u!, uout, hout, tmp)
    prob = DDEProblem{true}(dde_param, x0, h!, (T(t[1]), T(t[end])), p, constant_lags=sys.Tau)

    sol = DelayDiffEq.solve(prob, alg; saveat=t, kwargs...)

    x = sol.u::Vector{Vector{T}} # the states are labeled u in DelayDiffEq

    y = Matrix{T}(undef, noutputs(sys), length(t))
    d = Matrix{T}(undef, size(C2,1), length(t))
    # Build y signal (without d term)
    for k = 1:length(t)
        u!(uout, t[k])
        y[:,k] = C1*x[k] + D11*uout
        #d[:,k] = C2*x[k] + D21*uout
    end

    dtmp = Vector{T}(undef, size(C2,1))
    # Function to evaluate d(t)_i at an arbitrary time
    # X is continuous, so interpoate, u is not
    function dfunc!(tmp::Array{T1}, t, i) where T1
        tmp .= if t < t0
            T1(0)
        else
            sol(t) # solution object has built in interpolator of same order as solver.
        end
        u!(uout, t)
        return dot(view(C2,i,:),tmp) + dot(view(D21,i,:),uout)
    end

    # Account for the effect of the delayed d-signal on y
    for k=1:length(Tau)
        for j=1:length(t)
            di = dfunc!(tmp, t[j] - Tau[k], k)
            y[:, j] .+= view(D12,:, k) .* di
        end
    end

    return y', t, reduce(hcat ,x)'
end

# We have to default to something, look at the sys.P.P and delays
function _bounds_and_features(sys::DelayLtiSystem, plot::Val)
    ws, pz =  _bounds_and_features(sys.P.P, plot)
    logtau = log10.(abs.(sys.Tau))
    logtau = filter(x->x>4, logtau) # Ignore low frequency
    if isempty(logtau)
        return ws, pz
    end
    extreme = extrema(logtau)
    return [min(ws[1], floor(extreme[1]-0.2)), max(ws[2], ceil(extreme[2]+0.2))], pz
end

# Againm we have to do something for default vectors, more or less a copy from timeresp.jl
function _default_Ts(sys::DelayLtiSystem)
    if !isstable(sys.P.P)
        return 0.05   # Something small
    else
        ps = pole(sys.P.P)
        r = minimum([abs.(real.(ps));0]) # Find the fastest pole of sys.P.P
        r = min(r, minimum([sys.Tau;0])) # Find the fastest delay
        if r == 0.0
            r = 1.0
        end
        return 0.07/r
    end
end

iscontinuous(sys::DelayLtiSystem) = true

function Base.step(sys::DelayLtiSystem{T}, t::AbstractVector; kwargs...) where T
    nu = ninputs(sys)
    if t[1] != 0
        throw(ArgumentError("First time point must be 0 in step"))
    end
    u = (out, t) -> (t < 0 ? out .= 0 : out .= 1)
    x0=fill(zero(T), nstates(sys))
    if nu == 1
        y, tout, x = lsim(sys, u, t; x0=x0, kwargs...)
    else
        x = Array{T}(undef, length(t), nstates(sys), nu)
        y = Array{T}(undef, length(t), noutputs(sys), nu)
        for i=1:nu
            y[:,:,i], tout, x[:,:,i] = lsim(sys[:,i], u, t; x0=x0, kwargs...)
        end
    end
    return y, tout, x
end

function impulse(sys::DelayLtiSystem{T}, t::AbstractVector; alg=MethodOfSteps(BS3()), kwargs...) where T
    nu = ninputs(sys)
    iszero(sys.P.D12) || @warn("Impulse with a direct term from input to delay vector leads to poor accuracy.")
    if t[1] != 0
        throw(ArgumentError("First time point must be 0 in impulse"))
    end
    u = (out, t) -> (out .= 0)
    if nu == 1
        y, tout, x = lsim(sys, u, t; alg=alg, x0=sys.P.B[:,1], kwargs...)
    else
        x = Array{T}(undef, length(t), nstates(sys), nu)
        y = Array{T}(undef, length(t), noutputs(sys), nu)
        for i=1:nu
            y[:,:,i], tout, x[:,:,i] = lsim(sys[:,i], u, t; alg=alg, x0=sys.P.B[:,i], kwargs...)
        end
    end
    return y, tout, x
end

# Used for pade approximation
"""
`p2 = _linscale(p::Polynomial, a)`

Given a polynomial `p` and a number `a, returns the polynomials `p2` such that
`p2(s) == p(a*s)`.
"""
function _linscale(p::Polynomial, a)
    # This function should perhaps be implemented in Polynomials.jl
    coeffs_scaled = similar(p.coeffs, promote_type(eltype(p), typeof(a)))
    a_pow = 1
    coeffs_scaled[1] = p.coeffs[1]
    for k=2:length(p.coeffs)
        a_pow *= a
        coeffs_scaled[k] = p.coeffs[k]*a_pow
    end
    return Polynomial(coeffs_scaled)
end

# Coefficeints for Padé approximations
# PADE_Q_COEFFS = [Polynomial([binomial(N,i)*prod(N+1:2*N-i) for i=0:N]) for N=1:10]
const PADE_Q_COEFFS =  [[2, 1],
 [12, 6, 1],
 [120, 60, 12, 1],
 [1680, 840, 180, 20, 1],
 [30240, 15120, 3360, 420, 30, 1],
 [665280, 332640, 75600, 10080, 840, 42, 1],
 [17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1],
 [518918400, 259459200, 60540480, 8648640, 831600, 55440, 2520, 72, 1],
 [17643225600, 8821612800, 2075673600, 302702400, 30270240, 2162160, 110880, 3960, 90, 1],
 [670442572800, 335221286400, 79394515200, 11762150400, 1210809600, 90810720, 5045040, 205920, 5940, 110, 1]]

"""
    pade(τ::Real, N::Int)

Compute the `N`th order Padé approximation of a time-delay of length `τ`.
"""
function pade(τ::Real, N::Int)
    if !(1 <= N <= 10); error("Order of Padé approximation must be between 1 and 10. Got $N."); end

    Q = Polynomials.Polynomial(PADE_Q_COEFFS[N])

    return tf(_linscale(Q, -τ), _linscale(Q, τ)) # return Q(-τs)/Q(τs)
end

"""
    pade(G::DelayLtiSystem, N)

Approximate all time-delays in `G` by Padé approximations of degree `N`.
"""
function pade(G::DelayLtiSystem, N)
    ny, nu = size(G)
    nTau = length(G.Tau)
    X = append([ss(pade(τ,N)) for τ in G.Tau]...) # Perhaps append should be renamed blockdiag
    return lft(G.P.P, X)
end


text:
abstract type Parameterization{T} <: Map{T} end

gradientnorm(m::AbstractMap, t) = norm(gradient(m, t))

gradientnorm(m::DomainSets.UnitCircleMap{T}, t) where {T} = one(T)

export normal
"Evaluate the normal of the parameterized domain at the given point."
function normal(p::Parameterization, x)
    a = gradient(p, x)
    SVector(-a[2], a[1])./norm(a)
end


text:
#==============================================================================#
# http.jl
#
# HTTP Requests with retry/back-off and HTTPException.
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

import Base: show, UVError

http_status(e::HTTP.StatusError) = e.status
header(e::HTTP.StatusError, k, d="") = HTTP.header(e.response, k, d)
http_message(e::HTTP.StatusError) = String(e.response.body)
content_type(e::HTTP.StatusError) = HTTP.header(e.response, "Content-Type")

function http_request(request::AWSRequest)

    @repeat 4 try

        options = []
        if get(request, :return_stream, false)
            io = BufferStream()
            request[:response_stream] = io
            push!(options, (:response_stream, io))
        end

        verbose = debug_level - 1

        http_stack = HTTP.stack(redirect=false,
                                retry=false,
                                aws_authorization=false,
                                verbose=verbose)

        return HTTP.request(http_stack,
                            request[:verb],
                            HTTP.URI(request[:url]),
                            HTTP.mkheaders(request[:headers]),
                            request[:content];
                            #aws_service = request[:service],
                            #aws_region = request[:region],
                            #aws_access_key_id = request[:creds].access_key_id,
                            #aws_secret_access_key = request[:creds].secret_key,
                            #aws_session_token = request[:creds].token,
                            verbose = verbose,
                            require_ssl_verification=false,
                            options...)

    catch e

        @delay_retry if isa(e, Base.DNSError) ||
                        isa(e, HTTP.ParseError) ||
                        isa(e, HTTP.IOError) ||
                        #isa(e, EOFError) || FIXME needed ?
                       (isa(e, HTTP.StatusError) && http_status(e) >= 500) end
    end

    assert(false) # Unreachable.
end

function http_get(url::String)

    host = HTTP.URI(url).host

    http_request(@SymDict(verb = "GET",
                          url = url,
                          headers = ["Host" => host],
                          content = UInt8[]))
end

#==============================================================================#
# End of file.
#==============================================================================#


text:
struct Connect <: Packet
    header::UInt8
    protocol_name::String
    protocol_level::UInt8
    clean_session::Bool
    keep_alive::UInt16
    client_id::String
    will::Nullable{Message}
    username::Nullable{String}
    password::Nullable{Array{UInt8}}
    id::UInt16
end
Connect(clean_session::Bool, keep_alive::UInt16, client_id::String, will::Nullable{Message}, username::Nullable{String}, password::Nullable{Array{UInt8}}) = Connect(CONNECT, "MQTT", 0x04, clean_session, keep_alive, client_id, will, username, password, 0)

function write(s::IO, packet::Connect)
    mqtt_write(s, packet.protocol_name)
    mqtt_write(s, packet.protocol_level)
    clean_session_flag = convert(UInt8, packet.clean_session) << 1
    will_flag = convert(UInt8, !isnull(packet.will)) << 2
    will_qos = (isnull(packet.will) ? 0x00 : convert(UInt8, get(packet.will).qos)) << 3
    will_retain = (isnull(packet.will) ? 0x00 : convert(UInt8, get(packet.will).retain)) << 5
    username_flag = convert(UInt8, !isnull(packet.username)) << 6
    password_flag = convert(UInt8, !isnull(packet.password)) << 7
    connect_flags = clean_session_flag | will_flag | will_qos | will_retain | username_flag | password_flag
    mqtt_write(s, connect_flags)
    mqtt_write(s, packet.keep_alive)
    mqtt_write(s, packet.client_id)
    if !isnull(packet.will)
        mqtt_write(s, get(packet.will).topic)
        mqtt_write(s, get(packet.will).payload)
    end
    if !isnull(packet.username)
        mqtt_write(s, get(packet.username))
    end
    if !isnull(packet.password)
        mqtt_write(s, get(packet.password))
    end
end

has_id(packet::Connect) = true

Base.show(io::IO, x::Connect) = print(io, "connect[protocol_name: '", x.protocol_name, "'",
", protocol_level: ", x.protocol_level,
", clean_session: ", x.clean_session,
", keep_alive: ", x.keep_alive,
", client_id: '", x.client_id, "'",
", will: ", get(x.will, "none"),
", username: ", get(x.username, "none"),
", password: ", get(x.password, "none"), "]")

struct Connack <: Packet
    header::UInt8
    session_present::Bool
    return_code::UInt8
end
Connack(session_present::Bool, return_code::UInt8) = Connack(CONNACK, session_present, return_code)

function read(s::IO, flags::UInt8, ::Type{Connack})
    session_present = read(s, UInt8)
    return_code = read(s, UInt8)
    return Connack(convert(Bool, session_present), return_code)
end

Base.show(io::IO, x::Connack) = print(io, "connack[session_present: ", x.session_present, ", return_code: ", x.return_code ,"]")

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLMesh.jl" begin

   # function GLMesh()
   @testset "GLMesh" begin
      @test
      @test
      @test
      @test
   end

   # function GLMesh(primitive)
   @testset "GLMesh" begin
      @test
      @test
      @test
      @test
   end

   # function releaseGpuResources(mesh::GLMesh)
   @testset "releaseGpuResources" begin
      @test
      @test
      @test
      @test
   end

   # function computeNormal(p1::Point2d, p2::Point2d)
   @testset "computeNormal" begin
      @test
      @test
      @test
      @test
   end

   # function computeNormal(p0::Point3d,p1::Point3d,p2::Point3d)
   @testset "computeNormal" begin
      @test
      @test
      @test
      @test
   end

   # function getBoundingBox(mesh::GLMesh)
   @testset "getBoundingBox" begin
      @test
      @test
      @test
      @test
   end

   # function GLCuboid(box::Box3d)
   @testset "GLCuboid" begin
      @test
      @test
      @test
      @test
   end

   # function GLAxis(p0::Point3d,p1::Point3d)
   @testset "GLAxis" begin
      @test
      @test
      @test
      @test
   end

end


text:
import ..internal
function init()
    allocator = internal.rcl_get_default_allocator()
    internal.checkcall(internal.rcl_init(Int32(0), C_NULL, allocator))
end

shutdown() = internal.checkcall(internal.rcl_shutdown())

ok() :: Bool = internal.rcl_ok()

instance_id() :: UInt64 = internal.rcl_get_instance_id()


text:
"""
    ExternalDocstrings

ExternalDocstrings.jl is a helper for writing docstrings in markdown files.

# Usage

Create markdown files in `src/docs/` (where `src/MyPackage.jl` is the main package source
file).  Put the following line inside of `MyPackage` namespace:

```julia
ExternalDocstrings.@define_docstrings
```

It defines a docstring for name `MyPackage.\$name` using the markdown content in
`src/docs/\$name.md`.

# Extended help

## Markdown transformations

To use standard markdown (and CommonMark) while supporting special syntaxes for Julia
docstring, ExternalDocstrings.jl performs a couple of transformations:

(1) Code fence notation

    ```julia
    # ...
    ```

is transformed into

    ```jldoctest LABEL
    # ...
    ```

where `LABEL` is a label unique to the markdown file (i.e., all code blocks in one markdown
file are executed in the same session).

(2) To help doctests for non-REPL code block,

    nothing  # hide
    # output

is inserted at the end of code block if it does not look like a REPL session and does not
already have `# output`.

(3) `<kbd>KEY</kbd>` is replaced by `_KEY_`.

## Tips

### Disable doctest

To enable syntax highlighting without doctest, use use slightly different code fence
notations such as

    ```JULIA
    this_is_not_doctested() = nothing
    ```

### Vendoring

ExternalDocstrings.jl written as a single-source package to help vendoring.  For example, it
can be installed simply by:

```bash
wget https://raw.githubusercontent.com/tkf/ExternalDocstrings.jl/main/src/ExternalDocstrings.jl -O src/ExternalDocstrings.jl
```
"""
module ExternalDocstrings
#=
MIT License

Copyright (c) 2022 Takafumi Arakaki <aka.tkf@gmail.com> and contributors

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

function transform_docstring(doc::AbstractString, label)
    output = IOBuffer()
    input = IOBuffer(doc)
    while !eof(input)
        ln = readline(input)
        if startswith(ln, "```julia")
            print(output, "```jldoctest ", label, "\n")
            isrepl = false
            hasoutput = false
            while !eof(input)
                ln = readline(input)
                if startswith(ln, "```")
                    if !isrepl && !hasoutput
                        print(
                            output,
                            """
                            nothing  # hide
                            # output
                            """,
                        )
                    end
                    print(output, ln, "\n")
                    break
                end
                print(output, ln, "\n")
                if startswith(ln, "julia> ")
                    isrepl = true
                elseif startswith(ln, "# output")
                    hasoutput = true
                end
            end
        else
            ln = replace(ln, r"<kbd>(.*?)</kbd>" => s"_\1_")
            print(output, ln, "\n")
        end
    end
    return String(take!(output))
end

function define_docstrings(pkg::Module)
    if pathof(pkg) === nothing
        @warn """
        `define_docstrings` called on a non-package module `$pkg`.
        Not defining docstrings.
        """
        return
    end
    srcdir = dirname(pathof(pkg))
    docstrings = [nameof(pkg) => joinpath(dirname(srcdir), "README.md")]
    docsdir = joinpath(srcdir, "docs")
    if isdir(docsdir)
        for filename in readdir(docsdir)
            stem, ext = splitext(filename)
            ext == ".md" || continue
            name = Symbol(stem)
            name in names(pkg, all = true) || continue
            push!(docstrings, name => joinpath(docsdir, filename))
        end
    end
    n_auto_labels = 0
    for (name, path) in docstrings
        label = string(name)
        if match(r"^[a-z0-9_]+$"i, label) === nothing
            label = "$(nameof(pkg))$n_auto_labels"
            n_auto_labels += 1
        end

        include_dependency(path)
        doc = read(path, String)
        doc = transform_docstring(doc, label)

        ex = :($Base.@doc $doc $name)
        ex.args[2]::LineNumberNode
        ex.args[2] = LineNumberNode(1, Symbol(path))

        Base.eval(pkg, ex)
    end
end

macro define_docstrings()
    :(define_docstrings($(QuoteNode(__module__))))
end

end


text:
module TestGroupBy
include("preamble.jl")
using Transducers: DefaultInit

@testset begin
    @test foldl(right, GroupBy(string, Map(last), push!!), [1, 2, 1, 2, 3]) ==
          Dict("1" => [1, 1], "2" => [2, 2], "3" => [3])

    @test foldl(
        right,
        GroupBy(
            identity,
            opcompose(Map(last), Scan(+)),
            (_, x) -> x > 3 ? reduced(x) : x,
            nothing,
        ),
        [1, 2, 1, 2, 3],
    ) == Dict(2 => 4, 1 => 2)
end

@testset "post-groupby filtering" begin
    d = foldl(right, GroupBy(isodd, opcompose(Map(last), Filter(isodd)), +), 1:10)
    @test d == Dict(true => 25)
    @test d.state == Dict(true => 25, false => DefaultInit(+))
    @test valtype(d) <: Int
    @test valtype(d.state) <: Union{Int,typeof(DefaultInit(+))}
end

@testset "automatic asmonoid" begin
    @test foldl(right, GroupBy(identity, Map(first), (a, b) -> a + b), [1, 2, 1, 2, 3]) ==
          Dict(1 => 2, 2 => 4, 3 => 3)

    @test foldl(right, GroupBy(identity, (_, x) -> x), [1, 2, 1, 2, 3]) ==
          Dict(1 => 1 => 1, 2 => 2 => 2, 3 => 3 => 3)
end

@testset "GroupByViewDict" begin
    gd1 = foldl(right, GroupBy(string, Map(last), push!!), [1, 2, 1, 2, 3])
    gd2 = foldl(
        right,
        GroupBy(x -> gcd(x, 6), opcompose(Map(last), Filter(isodd)), push!!),
        1:10,
    )
    gd3 = foldl(
        right,
        opcompose(Filter(isodd), GroupBy(x -> gcd(x, 6), Map(last), push!!)),
        1:10,
    )
    gd4 = foldl(
        right,
        GroupBy(x -> gcd(x, 6), opcompose(Map(last), Unique(x -> gcd(x, 4))), push!!),
        1:10,
    )
    @test gd1 == Dict("1" => [1, 1], "2" => [2, 2], "3" => [3])
    @test gd2 == gd3
    @test gd4 == Dict(2 => [2, 4], 3 => [3], 6 => [6], 1 => [1])
    @test Dict(gd1) == Dict("1" => [1, 1], "2" => [2, 2], "3" => [3])
    @test Dict(gd2) == Dict(gd3)
    @test Dict(gd4) == Dict(2 => [2, 4], 3 => [3], 6 => [6], 1 => [1])
end

end  # module


text:
# Simple test on the weights code

using SubmodularMaximization
using SubmodularMaximization: channel_capacity_by_target,
                              channel_capacity_by_target_time

grid_size = 9
horizon = SubmodularMaximization.default_horizon

# place the target in the center of the grid
target_state = (5,5)

# Some representative and interesting weights sets
state_sets = [[(5,5), (5,5)],
              [(5,5), (5,10)],
              [(5,3), (5,7)],
              [(5,1), (5,9)]]

grid = Grid(grid_size, grid_size)
histogram_filter = Filter(grid, target_state)

# Create some uncertainty. Two steps is probably a reasonable approximation of
# the kind of target uncertainty that we will encounter
process_update!(histogram_filter, transition_matrix(grid))

for states in state_sets
  println()

  @show states

  problem = MultiRobotTargetTrackingProblem(states, [histogram_filter],
                                            grid=grid,
                                            horizon=horizon)

  # Investigate weights

  for ii = 1:length(states)
    @show target_weight = channel_capacity_by_target(problem, histogram_filter,
                                                     robot_index=ii)

    step_weights = map(1:horizon) do x
      @show channel_capacity_by_target_time(problem, histogram_filter,
                                            robot_index=ii, step=x)
    end
  end

  # Compute weight matrices
  methods = [channel_capacities_mcts,
             channel_capacities_by_target,
             channel_capacities_by_target_time]

  for method in methods
    @time weights = compute_weight_matrix(problem, channel_capacity_method=method)
    total_weight = sum(weights) / 2

    println("Weight: ", total_weight, " (", method, ")")
  end

end

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
push!(LOAD_PATH,"../src/")
using Documenter, ReinforcementLearningDatasets

makedocs(sitename="ReinforcementLearningDatasets")


text:
@testset "78.subsets.jl" begin
    @test subsets([1,2,3]) == [[],[1],[2],[1,2],[3],[1,3],[2,3],[1,2,3]]
    @test subsets([0]) == [[],[0]]
end


text:
using MagneticResonanceSignals: FixedString

@testset "FixedString" begin
    # Constructor
    @test FixedString{5}("abcd").data == (0x61, 0x62, 0x63, 0x64, 0x00)
    @test_throws ArgumentError("String \"α\" is not ascii") FixedString{5}("α")
    @test_throws ArgumentError("String \"abcd\" does not fit into 4-1 ascii characters") FixedString{4}("abcd")
    @test convert(FixedString{4}, FixedString{4}("abc")) === FixedString{4}("abc")
    @test convert(FixedString{4}, "abc") === FixedString{4}("abc")
    @test convert(String, FixedString{4}("abc"))::String == "abc"

    # show
    @test sprint(show, FixedString{10}("abcd")) == "FixedString{10}(\"abcd\")"
    @test sprint(show, FixedString{10}("ab\"cd")) == "FixedString{10}(\"ab\\\"cd\")"

    # IO
    writebuf = IOBuffer()
    @test write(writebuf, FixedString{6}("abcd")) == 6
    @test take!(writebuf) == UInt8['a', 'b', 'c', 'd', '\0', '\0']

    readbuf = IOBuffer(UInt8['a', 'b', 'c', 'd', '\0', '\0'])
    @test read(readbuf, FixedString{6}) == FixedString{6}("abcd")
    @test position(readbuf) == 6

    # AbstractString interface
    @test length(FixedString{10}("abc")) == 3
    @test codeunit(FixedString{10}("abc")) == UInt8
    @test ncodeunits(FixedString{10}("abc")) == 3
    @test ncodeunits(FixedString{2}((UInt8('a'),UInt8('b')))) == 2

    # Iteration
    @test collect(FixedString{10}("abc")) == ['a', 'b', 'c']
    @test collect(FixedString{10}("abc")) == ['a', 'b', 'c']
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# ---
# title: Regression plots
# cover: assets/regression_plots.png
# description: Linear and nonlinear regressions.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

x = rand(100)
y = @. randn() + x
df = (; x, y)
xy = data(df) * mapping(:x, :y)
layers = linear() + visual(Scatter)
draw(layers * xy)

#

x = rand(100)
y = @. randn() + 5 * x ^ 2
df = (; x, y)
xy = data(df) * mapping(:x, :y)
layers = smooth() + visual(Scatter)
fg = draw(layers * xy)

# save cover image #src
mkpath("assets") #src
save("assets/regression_plots.png", fg) #src