text:
# module TestBasics

using Test, Dates, Temporal

@testset "Initialization" begin
    @test Temporal.autocol(1:27)[end] == :AA
    @test Temporal.autoidx(2) == [today()-Day(1), today()]
    @test isa(X, TS)
    @test isa(Y, TS)
    @test size(X) == size(Y) == (N, K)
    @test X.values == Y.values == data
    @test X.index == collect(dates)
    @test Y.index == collect(times)
    @test size(TS(rand(N,K))) == (N,K)
    @test TS(X.values, X.index, ['A','B','C','D']) == X
    @test TS(X.values, X.index, ["A","B","C","D"]) == X
    @test TS(X.values[:,1], X.index, 'A') == X[:,1]
    @test TS(X.values[:,1], X.index, "A") == X[:,1]
    @test TS(X.values[1,:]', X.index[1], X.fields) == X[1]
    @test size(TS()) == (0,0)
    # underscores in colnames (issue 13)
    u = TS(rand(N), X.index, "under_score")
    @test u.fields[1] == :under_score
    show(stdout, TS(rand(252,4)))
    print("\n")
    # iterator protocol
    for (t, x) in X
        @test t in X.index
        @test size(x,1) == size(X,2)
    end
end

# end


text:
"""
Verlet integration is an integration method used to integrate newtons - law of motion. It is frequently used to find trajectories in molecular dynamics simulation.
The function takes `four` inputs viz,
- `f` : the differential equation
- `x0` : the initial condition. This is a Vector with the first element as initial value for position (x_0) and the second initial condition for velocity (v_0)
- `tspan`: is the time span for integration. It is a tuple (initial time, final time)
This functionr returns a tuple (x,t):
- `x` is the solution
- `t` is the array containing the time points
Reference:
- https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html
Contributed by: [Ved Mahajan](https://github.com/Ved-Mahajan)
"""
function verlet_integration(f, x0, tspan, Δt = 1.0e-3)
    s, e = tspan
    N = Int(floor((e - s) / Δt))
    x = Vector{Float64}(undef, N)
    t = collect(s:Δt:(e-Δt))
    x[1] = x0[1]
    x[2] = x0[1] + x0[2] * Δt + 0.5 * f(x0[1]) * (Δt)^2

    for i in 2:(N-1)
        x[i+1] = 2 * x[i] - x[i-1] + f(x[i]) * (Δt)^2
    end
    return x, t
end

ext: jl

max_issues_repo_name: Nikola-Mircic/Julia

max_forks_repo_name: Nikola-Mircic/Julia


text:
type class_homogeneous_newton <: abstract_newton_direction
	delta::Float64
	delta_mod::Array{Float64,1}

	direction::class_variables
	residuals::class_homogeneous_residuals
  nlp_vals::class_nlp_cache

	K::SparseMatrixCSC{Float64,Int64}
  W::woodbury_identity
  W_updated::Bool
  F_updated::Bool

	linear_system_solver::abstract_linear_system_solver

  function class_homogeneous_newton()
      return new();
  end
end

###################
# methods
###################

function initialize_newton!(newt::class_homogeneous_newton, nlp_eval::internal_AbstractNLPEvaluator, vars::class_variables, settings::class_settings)
      try
				dim = n(vars) + m(vars) + 1;
				newt.K = spzeros(dim, dim);

				newt.direction = class_variables(n(vars),m(vars));

				newt.linear_system_solver = settings.linear_system_solver;
				initialize!(newt.linear_system_solver, newt.K, :symmetric);

        newt.nlp_vals = class_nlp_cache();

        newt.delta = 0.0;
				update_delta_mod!(newt, vars, settings)

				newt.residuals = class_homogeneous_residuals();
        update_residuals!(newt.residuals, nlp_eval, vars, newt);
        newt.W_updated = false;
			catch e
				println("ERROR in class_homogeneous_newton.intialize")
				throw(e)
			end
end

function update_newton!(newt::class_homogeneous_newton, vars::class_variables, settings::class_settings)
    try
      nlp_vals = newt.nlp_vals;

      start_advanced_timer("Matrix creation");

      val_x_scaled = x_scaled(vars);
      H = nlp_vals.val_hesslag_prod;
      h = sparse(H * val_x_scaled);
      c = nlp_vals.val_gradc;
      A = nlp_vals.val_jac_a;
      b = sparse(nlp_vals.val_b);

      #res.H += newt.delta * speye(length(val_x_scaled))
      D_x = H + spdiagm( s(vars) ./ x(vars) ) #+ ( newt.delta + settings.diagonal_modification ) * speye( n(vars) )
      D_g = sparse(val_x_scaled' * H * val_x_scaled + kappa(vars) / tau(vars)) #+ settings.diagonal_modification;
      #D_g = val_x_scaled' * nlp_vals.val_hesslag_prod * val_x_scaled + vars.kappa() / vars.tau() + this.delta;
      D_z = settings.diagonal_modification * speye(m(vars));
      #v_1 = sparse(-c - h);
      #v_2 = sparse(c - h);
      #v_3 = a - A * val_x_scaled;

     # newt.K_true = [
     #   [ D_x  	c - h	-A' 	];
     #   [ -c -h' 	D_g 	b' 	];
     #   [ A 	 -b	    D_z	]
     #   ];

      newt.K[:,:] = [
        [ D_x  		-h 	   A' 	];
        [ -h'	    D_g 	 -b' 	];
        [ A 		-b	   -D_z	]
        ];

      newt.F_updated = false

      pause_advanced_timer("Matrix creation");
  catch e
      println("ERROR in class_homogeneous_newton.update_newton!")
      throw(e)
  end
end

function update_delta_mod!(newt::class_homogeneous_newton, vars::class_variables, settings::class_settings)

		# choose modification to diagonally dominant
		mod = zeros(n(vars)+1)
		for i = 1:(n(vars)+1)
			#mod[i] = max(sum(abs(newt.K[i,:])) - abs(newt.K[i,i]) - newt.K[i,i] + 1.0, 0.0)
			mod[i] = max(sum(abs(newt.K[i,(1:n(vars)+1)])) - abs(newt.K[i,i]) - newt.K[i,i] + 1e-8, 0.0)
		end

		dd_mod = [newt.delta * mod; -settings.diagonal_modification*ones(m(vars))]
		simple_mod = [newt.delta * ones(n(vars)); newt.delta * norm(x_scaled(vars),2)^2; -settings.diagonal_modification*ones(m(vars))]

		#@show simple_mod[n(vars)+1], dd_mod[n(vars)+1]

		if true
				newt.delta_mod = dd_mod
		else
				# simple delta choice that only affects
				newt.delta_mod = simple_mod
		end
end

function update_newton_diag!(newt::class_homogeneous_newton, vars::class_variables, settings::class_settings)
		update_delta_mod!(newt, vars,settings)

    H = newt.nlp_vals.val_hesslag_prod;
    D_x_diag = diag(H) + s(vars) ./ x(vars) +	newt.delta_mod[1:n(vars)]

    val_x_scaled = x_scaled(vars);
    D_g = val_x_scaled' * H * val_x_scaled + kappa(vars) / tau(vars) + newt.delta_mod[n(vars)+1]
    diag_mod = [ D_x_diag; D_g; newt.delta_mod[(n(vars)+2):(n(vars)+1+m(vars))]];

    for i = 1:size(newt.K,1)
        newt.K[i,i] = diag_mod[i];
    end

    return factorize_newton!(newt, vars)
end

function update_newton_diag_affine!(newt::class_homogeneous_newton, vars::class_variables, settings::class_settings)
    H = newt.nlp_vals.val_hesslag_prod;
    D_x_diag = diag(H) + s(vars) ./ x(vars) + newt.delta * x(vars).^(-2)

    val_x_scaled = x_scaled(vars);
    D_g = val_x_scaled' * H * val_x_scaled + kappa(vars) / tau(vars) + newt.delta / tau(vars)^2 #+ newt.delta * norm(val_x_scaled,2)^2;
    diag_mod = [ D_x_diag; D_g; -settings.diagonal_modification * ones(m(vars))];

    for i = 1:size(newt.K,1)
        newt.K[i,i] = diag_mod[i];
    end

    return factorize_newton!(newt, vars)
end

function factorize_newton!(newt::class_homogeneous_newton, vars::class_variables)
      start_advanced_timer("Factor");
      inertia = ls_factor!(newt.linear_system_solver, n(vars) + 1, m(vars));
      pause_advanced_timer("Factor");

      newt.W_updated = false;
      newt.F_updated = true;

      return inertia
end

function form_woodbury!(newt::class_homogeneous_newton, vars::class_variables)
    try
        @assert(newt.F_updated)

        start_advanced_timer("woodbury_factor");
        c = newt.nlp_vals.val_gradc;

        vec = [c; zeros( 1 + m(vars) )];
        len = n(vars) + m(vars) + 1;
        U = [-e_( 1 + n(vars), len) vec];
        V = [vec e_( 1 + n(vars), len)]';

        newt.W = woodbury_identity(newt.linear_system_solver, U, V);

        newt.W_updated = true;
        pause_advanced_timer("woodbury_factor");
  catch e
      println("ERROR form_woodbury")
      throw(e)
  end
end

function compute_newton_direction!(newt::class_homogeneous_newton, vars::class_variables, theta::class_theta)
			try
        num_vars = n(vars);
        num_const = m(vars);
        res = newt.residuals;

				mu = res.mu;

				xs = -x(vars) .* s(vars) + theta.mu * mu * ones(num_vars);
				tk = -tau(vars) * kappa(vars) + theta.mu * mu;

				rhs =
        [	(1.0 - theta.dual) * res.r_D + xs ./ x(vars);
				(1.0 - theta.dual) * res.r_G + tk / tau(vars);
				(1.0 - theta.primal) * res.r_P ];

        @assert(newt.W_updated == true)
        start_advanced_timer("Solve");

				sol = ls_solve(newt.W, rhs);

        pause_advanced_timer("Solve");

        # update direction from linear system solution
				dir = newt.direction;
				x(dir, sol[1:num_vars] );
				tau(dir, sol[ num_vars + 1] );
				y(dir, -sol[( num_vars + 2):( num_const +  num_vars + 1)] );

				s(dir,(xs - s(vars) .* x(dir)) ./ x(vars) );
				kappa(dir, (tk - kappa(vars) * tau(dir)) / tau(vars) );

				tmp = [x(dir); tau(dir); zeros(m(dir))]
				#@show tmp' * newt.K * tmp #+ norm(y(dir),2)*1e-8

				# is this direction valid ?
				check_for_wrong_vals(dir);

        return true
			catch e
				println("ERROR in class_newton_solver.compute_newton_direction!")
        #@show full(newt.K + newt.W.U * newt.W.V)
				throw(e)
			end
		end

function mu(newton_direction::class_homogeneous_newton, vars::class_variables)
    return (dot(s(vars),x(vars)) + dot(tau(vars),kappa(vars)))/(n(vars) + 1)
end

max_forks_repo_name: ohinder/Everest


text:
module SpectraJu

include("diffusion.jl")
include("integrale.jl")
include("functions.jl")
include("baseline.jl")

export peak_diffusion, model, IRdataprep, trapz

#From functions.jl
export poly, gaussiennes, normal_dist

export baseline

end # module


text:
@testset "lz4" begin

    test_in = """
        Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor
        incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud
        exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat
        """
    test_size = convert(Int32, length(test_in))
    bound = CodecLz4.LZ4_compressBound(test_size)
    bufsize = bound
    half_bufsize = floor(Cint, bufsize/2)

    @testset "bounds" begin
        bound = CodecLz4.LZ4_compressBound(test_size)
        @test bound > 0

        bound = CodecLz4.LZ4_compressBound(CodecLz4.LZ4_MAX_INPUT_SIZE)
        @test bound > 0

        bound = CodecLz4.LZ4_compressBound(CodecLz4.LZ4_MAX_INPUT_SIZE + 1)
        @test bound == 0
    end

    @testset "fast" begin
        buffer = Vector{UInt8}(undef, bufsize)
        dec_buffer = Vector{UInt8}(undef, test_size)
        result = CodecLz4.LZ4_compress_fast(test_in, buffer, test_size, bufsize, 12)

        result = CodecLz4.LZ4_decompress_safe(buffer, dec_buffer, result, test_size)
        @test unsafe_string(pointer(dec_buffer), result) == test_in

        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_compress_fast(test_in, buffer, test_size, half_bufsize, 12)
    end

    @testset "destSize" begin
        buffer = Vector{UInt8}(undef, bufsize)
        dec_buffer = Vector{UInt8}(undef, test_size)
        sz = Ref{Int32}(test_size)

        compressed = CodecLz4.LZ4_compress_destSize(test_in, buffer, sz, half_bufsize)
        @test compressed == half_bufsize
        @test sz[] > 0

        result = CodecLz4.LZ4_decompress_safe(buffer, dec_buffer, compressed, sz[])
        @test result == sz[]
        @test unsafe_string(pointer(dec_buffer), result) == test_in[1:sz[]]
    end

    @testset "streams" begin
        buffer = Vector{UInt8}(undef, bufsize)
        dec_buffer = Vector{UInt8}(undef, test_size)
        comp_str = CodecLz4.LZ4_createStream()
        dec_str = CodecLz4.LZ4_createStreamDecode()
        try
            @test CodecLz4.LZ4_resetStream(comp_str) === nothing
            @test CodecLz4.LZ4_setStreamDecode(dec_str, Ptr{UInt8}(C_NULL), Cint(0)) == 1

            result = CodecLz4.LZ4_compress_fast_continue(comp_str, test_in, buffer, test_size, bufsize, 12)
            @test result > 0
            result = CodecLz4.LZ4_decompress_safe_continue(dec_str, buffer, dec_buffer, result, test_size)

            @test result == test_size
            @test unsafe_string(pointer(dec_buffer), result) == test_in
        finally
            @test CodecLz4.LZ4_freeStream(comp_str) == 0
            @test CodecLz4.LZ4_freeStreamDecode(dec_str) == 0
        end
    end

    @testset "uninitialized" begin
        buffer = Vector{UInt8}(undef, bufsize)
        dec_buffer = Vector{UInt8}(undef, test_size)
        stream = Ref{Ptr{CodecLz4.LZ4_stream_t}}(C_NULL)

        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_resetStream(stream[])
        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_compress_fast_continue(stream[], test_in, buffer, test_size, bufsize)
        @test CodecLz4.LZ4_freeStream(stream[]) == 0

        stream = Ref{Ptr{CodecLz4.LZ4_streamDecode_t}}(C_NULL)

        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_setStreamDecode(stream[])
        @test_throws CodecLz4.LZ4Exception CodecLz4.LZ4_decompress_safe_continue(stream[], test_in, buffer, test_size, bufsize)
        @test CodecLz4.LZ4_freeStreamDecode(stream[]) == 0
    end

end


text:
using Test
using DSP
using LinearAlgebra
using Statistics
using Unitful
using Unitful: ns, ms, µs, s, Hz, kHz, MHz, GHz, THz, m
using SampleArrays
using Loudspeakers

test_dir = "../data/test/vituix"

@testset "Loudspeakers.jl" begin
    include("utils.jl")

    # check interpolation on orbitals, e.g., 10° and 20° to estimate 15°
    include("interpolation.jl")

    # test spinorama curves (listening window, power response, ...)
    include("vituix_spinorama.jl")

    # test single driver transforms (translation, rotation, tilt, listening distance)
    # TODO: different listening distances
    include("vituix_transforms.jl")

    # test FIR filters applied to single drivers
    include("vituix_fir.jl")

    # test multiple drivers
    # TODO: check different listening distances - it seems there some differences (try even < 1m and < 2m)
    include("vituix_multi_driver.jl")

    # TODO: projection point transforms? How to get data from Vituix? Use reference angle, or Room?
end


text:
using FFTW

@testset "Chebyshev Routines" begin
    import PotentialFlow.Plates.Chebyshev

    @testset "Quadrature" begin
        N = 127
        w = Chebyshev.clenshaw_curtis_weights(N)
        @test sum(w) ≈ 2

        x = Chebyshev.nodes(N)
        @test w ⋅ x ≤ eps()
        @test w ⋅ x.^2 ≈ 2/3
        @test w ⋅ exp.(x) ≈ (exp(1) - exp(-1))
    end

    @testset "Series" begin
        # Check evaluation off offset Chebyshev series
        s = rand(2)
        @test Chebyshev.firstkind([rand(), 2.0, 1.0], s, 1)     ≈ @. 2s^2 + 2s - 1
        @test Chebyshev.firstkind([rand(), rand(), 1.0], s, 2)  ≈ @. 2s^2 - 1
        @test Chebyshev.secondkind([rand(), 2.0, 1.0], s, 1)    ≈ @. 4s^2 + 4s - 1
        @test Chebyshev.secondkind([rand(), rand(), 1.0], s, 2) ≈ @. 4s^2 - 1
    end

    @testset "Transforms" begin

        N = rand(128:512)
        R = zeros(Float64, N)
        y = zeros(Float64, N)
        @test all(rand(0:(N÷2), 10)) do n
            x = [cos(n*θ) for θ in range(π, 0, length=N)]
            Chebyshev.transform!(R, x)
            Chebyshev.inv_transform!(y, R)
            (y ≈ x) && (R[n+1] ≈ 1.0) && (sum(R) ≈ 1.0)
        end

        n₁ = rand(0:(N÷2))
        n₂ = rand(0:(N÷2))
        x₁ = [cos(n₁*θ) for θ in range(π, 0, length=N)]
        x₂ = [cos(n₂*θ) for θ in range(π, 0, length=N)]

        C = x₁ .+ im.*x₂
        I = zeros(Float64, N)

        Chebyshev.transform!(R, x₁)
        Chebyshev.transform!(I, x₂)
        Chebyshev.transform!(C)

        @test real.(C) ≈ R
        @test imag.(C) ≈ I

        plan! = FFTW.plan_r2r!(C, FFTW.REDFT00)
        @allocated Chebyshev.transform!(C, plan!)
        @test 0 == (@allocated Chebyshev.transform!(C, plan!))

        @testset "Planned Transforms" begin
            x = Chebyshev.nodes(N)

            K  = Chebyshev.plan_transform(ones(N))
            K! = Chebyshev.plan_transform!(ones(N))
            @test_throws BoundsError K*ones(N+1)
            @test_throws BoundsError K\ones(N+1)

            C = K*ones(N)
            @test C[1] ≈ 1
            @test norm(C[2:end]) ≤ 128eps()

            C = K*x
            @test norm(C[1]) ≤ 128eps()
            @test C[2] ≈ 1
            @test norm(C[3:end]) ≤ 128eps()

            y = @. 16x^5 - 20x^3 + 5x
            C = K*y
            @test norm(C[1:5]) ≤ 128eps()
            @test C[6] ≈ 1
            @test norm(C[7:end]) ≤ 128eps()

            ỹ = K \ C
            @test y ≈ ỹ

            # Inverse transforms should be the same as simply
            # evaluating the series
            @test Chebyshev.firstkind(C, x) ≈ ỹ

            K! \ C
            @test y ≈ C

            # Make sure transforms work with even number of nodes
            N = 128
            x = Chebyshev.nodes(N)
            K! = Chebyshev.plan_transform!(x)
            y = @. 16x^5 - 20x^3 + 5x
            C = copy(y)
            K! * C
            ỹ = Chebyshev.inv_transform(C)
            @test y ≈ ỹ

            K! \ C
            @test y ≈ C

            # Make sure transforms work with odd number of nodes
            N = 127
            x = Chebyshev.nodes(N)
            K! = Chebyshev.plan_transform!(x)
            y = @. 16x^5 - 20x^3 + 5x
            C = copy(y)
            K! * C
            ỹ = Chebyshev.inv_transform(C)
            @test y ≈ ỹ

            K! \ C
            @test y ≈ C
        end
    end

end


text:
module AutodiffTests

using Test

using Gridap.Arrays

function user_cell_energy(cell_u)
  function f(u)
    e = zero(eltype(u))
    for ui in u
      e += ui^2
    end
    e
  end
  lazy_map(f,cell_u)
end

function user_cell_residual(cell_u)
  lazy_map(Broadcasting(i->2*i),cell_u)
end

function user_cell_jacobian(cell_u)
  function f(u)
    n = length(u)
    j = zeros(n,n)
    for i in 1:n
      j[i,i] = 2
    end
    j
  end
  lazy_map(f,cell_u)
end

L = 10
l = 4 # Do not use a number <13 (too easy for ForwardDiff)

cell_u = [rand(l) for i in 1:L]
cell_e = user_cell_energy(cell_u)
cell_r = user_cell_residual(cell_u)
cell_j = user_cell_jacobian(cell_u)
cell_h = cell_j

cell_j_auto = autodiff_array_jacobian(user_cell_residual,cell_u)
cell_r_auto = autodiff_array_gradient(user_cell_energy,cell_u)
cell_h_auto = autodiff_array_hessian(user_cell_energy,cell_u)
test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)
test_array(cell_h_auto,cell_h)

ids = [3,4,1,2]

function user_cell_energy_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_energy(cell_u_Γ)
end

function user_cell_residual_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_residual(cell_u_Γ)
end

function user_cell_jacobian_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_jacobian(cell_u_Γ)
end

cell_e_Γ = user_cell_energy_Γ(cell_u)
cell_r_Γ = user_cell_residual_Γ(cell_u)
cell_j_Γ = user_cell_jacobian_Γ(cell_u)
cell_h_Γ = cell_j_Γ

cell_r_Γ_auto = autodiff_array_gradient(user_cell_energy_Γ,cell_u,ids)
cell_j_Γ_auto = autodiff_array_jacobian(user_cell_residual_Γ,cell_u,ids)
cell_h_Γ_auto = autodiff_array_hessian(user_cell_energy_Γ,cell_u,ids)

test_array(cell_r_Γ_auto,cell_r_Γ)
test_array(cell_j_Γ_auto,cell_j_Γ)
test_array(cell_h_Γ_auto,cell_h_Γ)

end # module


text:
import Base: show

function show(io::IO, mime::MIME"text/plain", cur_plot::AbstractPlot)
  display(_show(cur_plot))
end

function show(io::IO, mime::MIME, cur_plot::AbstractPlot)
  show(io, mime, _show(cur_plot))
end

function _show(cur_plot::AbstractPlot)
  cur_size, cur_data, cur_layout, cur_config = _show_helper(cur_plot)
  render(cur_size, cur_data, cur_layout, cur_config)
end

function _show_helper(simple_plot::SimplePlot)
  global _plot = simple_plot

  simple_layout = parse_layout(simple_plot)

  simple_data = simple_plot.data
  simple_config = simple_plot.config

  simple_size = simple_plot.size

  return (
    simple_size, simple_data, simple_layout, simple_config
  )
end

function _show_helper(composite_plot::CompositePlot)
  composite_config = nothing
  composite_data = []

  composite_layout_base = nothing
  composite_layout_extra = Dict()

  for (cur_index, simple_plot) in enumerate(composite_plot.plots)
    if isnothing(composite_config)
      composite_config = simple_plot.config
    else
      @assert composite_config == simple_plot.config
    end

    cur_layout_base = deepcopy(simple_plot.layout)

    if haskey(cur_layout_base, "xaxis")
      cur_xaxis = pop!(cur_layout_base, "xaxis")
    else
      cur_xaxis = Dict()
    end

    if haskey(cur_layout_base, "yaxis")
      cur_yaxis = pop!(cur_layout_base, "yaxis")
    else
      cur_yaxis = Dict()
    end

    if ( cur_index == 1 )
      @assert isnothing(composite_layout_base)
      composite_layout_base = cur_layout_base

      isempty(cur_xaxis) || ( composite_layout_extra["xaxis"] = cur_xaxis )
      isempty(cur_yaxis) || ( composite_layout_extra["yaxis"] = cur_yaxis )

      append!(composite_data, simple_plot.data)
      continue
    end

    for (cur_key, cur_value) in cur_layout_base
      if haskey(composite_layout_base, cur_key)
        @assert composite_layout_base[cur_key] == cur_value
      else
        composite_layout_base[cur_key] = cur_value
      end
    end

    isempty(cur_xaxis) || ( composite_layout_extra["xaxis$( cur_index )"] = cur_xaxis )
    isempty(cur_yaxis) || ( composite_layout_extra["yaxis$( cur_index )"] = cur_yaxis )

    tmp_data = deepcopy(simple_plot.data)
    for cur_datum in tmp_data
      cur_datum["xaxis"] = "x" * string(cur_index)
      cur_datum["yaxis"] = "y" * string(cur_index)
    end
    append!(composite_data, tmp_data)
  end

  composite_layout = merge(composite_layout_base, composite_layout_extra)

  width_per_height = composite_plot.size[1]
  width_per_height /= composite_plot.size[2]

  plots_count = length(composite_plot.plots)

  cur_columns = nothing
  best_blanks = Inf
  for cur_method in [ceil, floor]
    tmp_columns = Int(cur_method(sqrt( plots_count * width_per_height )))
    tmp_rows = Int(ceil(plots_count / tmp_columns))

    num_blanks = tmp_columns * tmp_rows - plots_count
    ( num_blanks < best_blanks ) || continue

    best_blanks = num_blanks
    cur_columns = tmp_columns
  end
  cur_rows = Int(ceil(plots_count / cur_columns))

  composite_layout["grid"] = Dict(
    "rows" => string(cur_rows),
    "columns" => string(cur_columns),
    "pattern" => "independent"
  )

  composite_size = composite_plot.size

  return (
    composite_size, composite_data, composite_layout, composite_config
  )
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m8_opt_rho13_5(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m8_opt_rho13_5!(A_copy)
end

@inline function exp_ps_m8_opt_rho13_5!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 Ba9 Bb9 B9 T2k11
    # Computing Bb2 = x*I+x*A
    coeff1=0.42031771622120523
    coeff2=0.154626052690562
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.42031771622120523
    coeff2=0.154626052690562
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.09654809779701808
    coeff2=0.4657361097626277
    coeff3=0.19731022777550086
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.9062397883450928
    coeff2=0.05669626301307128
    coeff3=0.12302744324936216
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.06178357797354583
    coeff2=0.06738540542304026
    coeff3=0.7329752440067754
    coeff4=0.5453304515439815
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.7648819021317048
    coeff2=0.15551270544310056
    coeff3=0.13872927228389606
    coeff4=0.031072070853003377
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0037226968956481066
    coeff2=0.1513971118834432
    coeff3=0.16951280499357504
    coeff4=0.7173684223758208
    coeff5=0.6151562256858113
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.4984122280039027
    coeff2=0.3084467039511569
    coeff3=0.32464232547450916
    coeff4=0.11022586483948126
    coeff5=0.00562561032811143
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-2.47688825760804e-9
    coeff2=0.0002423691784933533
    coeff3=-0.008759615364757254
    coeff4=-0.024971360914667298
    coeff5=0.20141416233600964
    coeff6=0.9838034227912541
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-7.332920860083383e-8
    coeff2=3.6275534419311994e-10
    coeff3=1.599720660809858e-11
    coeff4=5.386947932212966e-11
    coeff5=3.696432230206371e-12
    coeff6=1.7604767541951467e-14
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=5.241383560289876e-5
    coeff2=0.007243536894625013
    coeff3=0.01054251273835961
    coeff4=-0.12590635459023522
    coeff5=-0.35076815913299725
    coeff6=0.9784008568588043
    coeff7=4.756683002701541e-7
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.004732216406789163
    coeff2=1.2338220978385536e-5
    coeff3=1.805471807800756e-5
    coeff4=6.965400723989983e-6
    coeff5=1.4212591399083807e-6
    coeff6=1.2164967481538e-7
    coeff7=1.0000000002008125
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9900217594371604
    coeff2=-0.02132723543710295
    coeff3=0.0895038182070465
    coeff4=0.4332397850086141
    coeff5=0.6230790774477177
    coeff6=0.836798518267107
    coeff7=6.715240076479665e-7
    coeff8=0.03410244422083362
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.09171588796170577
    coeff2=0.00978234693719379
    coeff3=0.01590736865892108
    coeff4=0.010821193755227018
    coeff5=0.002975380215829149
    coeff6=0.005577510585906947
    coeff7=-7.0435154500906185e-6
    coeff8=1.0199299719677057
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing Ba9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.5991807949797006
    coeff2=0.15278813100979383
    coeff3=0.3877431664280993
    coeff4=0.15662767864214414
    coeff5=0.3206562914411985
    coeff6=0.06358374798917776
    coeff7=-0.4636979253428603
    coeff8=-0.005883392744500455
    coeff9=1.0870379046019574
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.8838850843889047
    coeff2=0.13220125487009166
    coeff3=0.3359270098679286
    coeff4=0.13565211594408094
    coeff5=0.27797106138779576
    coeff6=0.05519918645255441
    coeff7=0.39661780513350103
    coeff8=0.007845590149896397
    coeff9=0.9425352289352358
    # Smart lincomb recycle B5
    memslots7 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots7,true,I*coeff1,true,true)
    # Computing B9 with operation: mult
    mul!(memslots3,memslots2,memslots7)
    # Deallocating Ba9 in slot 2
    # Deallocating Bb9 in slot 7
    # Computing T2k11 = x*I+x*A+x*B2+x*B3+x*B4+x*B6+x*B7+x*B8+x*B9
    coeff1=0.03877713825381048
    coeff2=-9.60542079781634e-5
    coeff3=-0.0001512696450918944
    coeff4=-7.050561377429653e-5
    coeff5=-9.005499714502881e-5
    coeff6=-0.335908979618945
    coeff7=0.0029363288440026106
    coeff8=-0.00026219521556040433
    coeff9=1.0181552414422717
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots8 .+ coeff7.*memslots9 .+ coeff8.*memslots10 .+ coeff9.*memslots3
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    # Deallocating B9 in slot 3
    return memslots1 # Returning T2k11
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using OrdinaryDiffEq, ModelingToolkit, Test
using GalacticOptim, Optim, RecursiveArrayTools

N = 32
const xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = ModelingToolkit.ifelse(a == N+1, 1, ModelingToolkit.ifelse(a == 0, N, a))
function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end

# Test with tuple parameters
p = (3.4, 1., 10., step(xyd_brusselator))

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
u0 = init_brusselator_2d(xyd_brusselator)

# Test with 3-tensor inputs
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0,(0.,11.5),p)

modelingtoolkitize(prob_ode_brusselator_2d)

## Optimization

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

prob = OptimizationProblem(rosenbrock,x0,p)
sys = modelingtoolkitize(prob) # symbolicitize me captain!

prob = OptimizationProblem(sys,x0,p,grad=true,hess=true)
sol = solve(prob,NelderMead())
@test sol.minimum < 1e-8

sol = solve(prob,BFGS())
@test sol.minimum < 1e-8

sol = solve(prob,Newton())
@test sol.minimum < 1e-8

## SIR System Regression Test

β = 0.01# infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, λ_R, λ_D]...)

# regional contact matrix and regional population

## regional contact matrix
regional_all_contact_matrix = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age
N = [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...)

# Time
final_time = 20
𝒯 = (1.0,final_time);

function SIRD_ac!(du,u,p,t)
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initialize this parameter (death probability stratified by age, taken from literature)

    δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
    δ = vcat(repeat([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))

    C = regional_all_contact_matrix

    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1:4*5]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
    dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1:4*5]

    # Force of infection
    Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]]

    # System of equations
    @. dS = -Λ*S
    @. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
    @. dR = λ_R*(1-δ)*I
    @. dD = λ_D*δ*I
    @. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]

end;

# create problem and check it works
problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
@time solution = solve(problem, Tsit5(), saveat = 1:final_time);

problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem)
fast_problem = ODEProblem(sys,ℬ, 𝒯, 𝒫 )
@time solution = solve(fast_problem, Tsit5(), saveat = 1:final_time)

## Issue #778

r0 = [1131.340, -2282.343, 6672.423]
v0 = [-5.64305, 4.30333, 2.42879]
Δt = 86400.0*365
μ = 398600.4418
rv0 = ArrayPartition(r0,v0)

function f(dy, y, μ, t)
    r = sqrt(sum(y[1,:].^2))
    dy[1,:] = y[2,:]
    dy[2,:] = -μ .* y[1,:] / r^3
end

prob = ODEProblem(f, rv0, (0.0, Δt), μ)
sol = solve(prob, Vern8())

modelingtoolkitize(prob)


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

using SlepcWrap
using Documenter
using Literate

# Generate examples
example_src = joinpath(@__DIR__,"..","example")
example_dir = joinpath(@__DIR__,"src","example")
Sys.rm(example_dir; recursive=true, force=true)
Literate.markdown(joinpath(example_src, "helmholtz.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
Literate.markdown(joinpath(example_src, "helmholtz_fancy.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
Literate.markdown(joinpath(example_src, "complex.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells

makedocs(;
    modules=[SlepcWrap],
    authors="bmxam",
    sitename="SlepcWrap.jl",
    clean=true,doctest=false,
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/bmxam/SlepcWrap.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "example/helmholtz.md",
        ],
        "Fancy examples" => Any[
            "example/helmholtz_fancy.md",
            "example/complex.md",
        ],
        "API Reference" => Any[
            "api/init.md",
            "api/eps.md",
        ],
        "API fancy" => "api/fancy/fancy.md",
    ],
)

deploydocs(;
    repo="github.com/bmxam/SlepcWrap.jl.git",
    push_preview = true
)


text:
function nominal_errors(network, X, Y; coeffs = [-0.74, -0.44])
    """ Computes absolute errors with no adversarial perturbations
        Args:
            - network: neural network model (neural verification format)
            - X: training examples (each column is an example)
            - Y: labels (each column is a label)
            ---------------
            - coeffs: coefficients to get output (linear combination of network outputs)
    """

    N_samples = size(Y, 2)

    truevals = vec(coeffs' * Y)
    predvals = [coeffs' * compute_output(network, X[:, i]) for i = 1:N_samples]

    return abs.(truevals - predvals)
end


text:
# export export_SDP_Instance

function export_SDP_Instance(sdp_pb::SDP_Problem, path::String)
    # Export blocks of constraints
    blocks_file = joinpath(path, "matrix.sdp")
    !isfile(blocks_file) || rm(blocks_file)

    open(blocks_file, "a") do fblocks
        print_matrices(fblocks, sdp_pb.matrices)
    end

    # Export linear part of constraints
    lin_file = joinpath(path, "lin.sdp")
    !isfile(lin_file) || rm(lin_file)

    open(lin_file, "a") do flin
        print_scalarlinear(flin, sdp_pb.linear)
    end

    # Export constants of constraints
    cst_file = joinpath(path, "const.sdp")
    !isfile(cst_file) || rm(cst_file)

    open(cst_file, "a") do fcst
        print_constant(fcst, sdp_pb.cst_ctr)
    end

    # Export bloc types
    types_file = joinpath(path, "types.sdp")
    !isfile(types_file) || rm(types_file)

    open(types_file, "a") do ftypes
        print_sdpblocks(ftypes, sdp_pb.name_to_sdpblock)
    end
end

function print_matrices(io::IO, matrices::SortedDict{Tuple{SDP_CtrObjName, String, String, String}, Float64}, indentedprint=true)
    MPC.print_blocksfile_header(io)

    ctrkey1len =  max( maximum(x->length(x[1][1]), keys(matrices)), length("#j_conj"))
    ctrkey2len =  max( maximum(x->length(x[1][2]), keys(matrices)), length("j_expl"))
    ctrkey3len = max( maximum(x->length(x[1][3]), keys(matrices)), length("clique"))
    blocklen =  max( maximum(x->length(x[2]), keys(matrices)), length("Zi"))
    rowlen =    max( maximum(x->length(x[3]), keys(matrices)), length("k"))
    collen =    max( maximum(x->length(x[4]), keys(matrices)), length("l"))
    MPC.print_string(io, "#j_conj", ctrkey1len, indentedprint=indentedprint)
    MPC.print_string(io, "j_expl", ctrkey2len, indentedprint=indentedprint)
    MPC.print_string(io, "clique", ctrkey3len, indentedprint=indentedprint)
    MPC.print_string(io, "Zi", blocklen, indentedprint=indentedprint)
    MPC.print_string(io, "k", rowlen, indentedprint=indentedprint)
    MPC.print_string(io, "l", collen, indentedprint=indentedprint)
    MPC.print_string(io, "Real(A_ij[k, l])", 23, indentedprint=indentedprint)
    MPC.print_string(io, "Imag(A_ij[k, l])", 23, indentedprint=indentedprint)
    println(io)

    for (((ctrname1, ctrname2, ctrname3), blockname, var1, var2), λ) in matrices

        MPC.print_string(io, ctrname1, ctrkey1len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname2, ctrkey2len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname3, ctrkey3len, indentedprint=indentedprint)
        MPC.print_string(io, blockname, blocklen, indentedprint=indentedprint)
        MPC.print_string(io, var1, rowlen, indentedprint=indentedprint)
        MPC.print_string(io, var2, collen, indentedprint=indentedprint)
        @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
    end
end

function print_scalarlinear(io::IO, linear::SortedDict{Tuple{SDP_CtrObjName, String}, Float64}, indentedprint=true)
    MPC.print_linfile_header(io)

    (length(linear) == 0) && return

    ctrkey1len =  max( maximum(x->length(x[1][1]), keys(linear)), length("#j_conj"))
    ctrkey2len =  max( maximum(x->length(x[1][2]), keys(linear)), length("j_expl"))
    ctrkey3len = max( maximum(x->length(x[1][3]), keys(linear)), length("clique"))
    varlen =    max( maximum(x->length(x[2]), keys(linear)), length("x[k]"))
    MPC.print_string(io, "#j_conj", ctrkey1len, indentedprint=indentedprint)
    MPC.print_string(io, "j_expl", ctrkey2len, indentedprint=indentedprint)
    MPC.print_string(io, "clique", ctrkey3len, indentedprint=indentedprint)
    MPC.print_string(io, "k", varlen, indentedprint=indentedprint)
    MPC.print_string(io, "Real(b_j[k])", 23, indentedprint=indentedprint)
    MPC.print_string(io, "Imag(b_j[k])", 23, indentedprint=indentedprint)
    println(io)

    for (((ctrname1, ctrname2, ctrname3), varname), λ) in linear

        MPC.print_string(io, ctrname1, ctrkey1len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname2, ctrkey2len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname3, ctrkey3len, indentedprint=indentedprint)
        MPC.print_string(io, varname, varlen, indentedprint=indentedprint)
        @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
    end
end

function print_constant(io::IO, cst::SortedDict{SDP_CtrObjName, Float64}, indentedprint=true)
    MPC.print_cstfile_header(io)

    (length(cst) == 0) && return

    ctrkey1len =  max( maximum(x->length(x[1]), keys(cst)), length("#j_conj"))
    ctrkey2len =  max( maximum(x->length(x[2]), keys(cst)), length("j_expl"))
    ctrkey3len = max( maximum(x->length(x[3]), keys(cst)), length("clique"))

    MPC.print_string(io, "#j_conj", ctrkey1len, indentedprint=indentedprint)
    MPC.print_string(io, "j_expl", ctrkey2len, indentedprint=indentedprint)
    MPC.print_string(io, "clique", ctrkey3len, indentedprint=indentedprint)
    MPC.print_string(io, "Real(c_j)", 23, indentedprint=indentedprint)
    MPC.print_string(io, "Imag(c_j)", 23, indentedprint=indentedprint)
    println(io)

    for ((ctrname1, ctrname2, ctrname3), λ) in cst

        MPC.print_string(io, ctrname1, ctrkey1len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname2, ctrkey2len, indentedprint=indentedprint)
        MPC.print_string(io, ctrname3, ctrkey3len, indentedprint=indentedprint)
        @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
    end
end

function print_sdpblocks(io::IO, name_to_sdpblock::SortedDict{String, SDP_Block}, indentedprint=true)
    MPC.print_typesfile_header(io)

    cstrlen = max( maximum(x->length(x), keys(name_to_sdpblock)), length("#Zi"))

    MPC.print_string(io, "#Zi", cstrlen, alignright=false); println(io, " type")

    for (blockname, block) in name_to_sdpblock
        MPC.print_string(io, blockname, cstrlen)
        println(io, " SDP")
    end
end


text:
function primer74BFS(G, ime, prezime)
	s = -1;
	suggestions = [];

	for i = 1:length(G.V)
		if G.V[i].ime==ime && G.V[i].prezime==prezime
			s = i
			break
		end
	end

	if s == -1
		println("Ne postoji korisnik za zadatim imenom.")
		return []
	else
		println("Id zadatog korisnika je: $s")
	end

	for u = 1:length(G.V)
	   if u!=s
		G.V[u].color = 'W'
		G.V[u].d = typemax(Int64)
		G.V[u].pred = -1
	   end
	end

        G.V[s].color = 'G'
        G.V[s].d = 0
        G.V[s].pred = -1

        Q = Int[]
        push!(Q, s)
	while ~isempty(Q)
		u = Q[1]
		deleteat!(Q, 1)
		if G.V[u].d < 2
		     for v in findall(G.AdjMatrix[u,:] .==1)
			   if G.V[v].color=='W'
			      G.V[v].color = 'G'
			      G.V[v].d = G.V[u].d + 1
			      G.V[v].pred = u
			      push!(Q, v)
			   end
		     end
		end
		if G.V[u].d == 2 #prijatelji od prijatelja se nalaze na distanci 2
			push!(suggestions, G.V[u]);
		end
		G.V[u].color = 'B';
	end
	return suggestions
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains package wide type aliases, constants, macros, utility methods, etc.

################################################# IMPORT/EXPORT ############################################################

export
# Type Aliases
VertexID, EdgeID, VertexList, EdgeList

################################################# TYPE ALIASES #############################################################

""" Datatype used to store vertex id numbers """
typealias VertexID Int

""" Datatype used to store edges """
typealias EdgeID Pair{VertexID,VertexID}

""" A list of Vertex IDs """
typealias VertexList AbstractVector{VertexID}

""" A list of Edge IDs """
typealias EdgeList AbstractVector{EdgeID}


text:
# This file was generated, do not modify it. # hide
pprint(first(df,10))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"The empty domain with elements of type `T`."
struct EmptySpace{T} <: Domain{T}
end

const AnyEmptySpace = EmptySpace{Any}

EmptySpace() = EmptySpace{Float64}()
EmptySpace(::Type{T}) where {T} = EmptySpace{T}()

similardomain(::EmptySpace, ::Type{T}) where {T} = EmptySpace{T}()

emptyspace(d) = emptyspace(eltype(d))
emptyspace(::Type{T}) where {T} = EmptySpace{T}()

indomain(x::T, d::EmptySpace{T}) where {T} = false
approx_indomain(x, d::EmptySpace, tolerance) = in(x, d)

show(io::IO, d::EmptySpace) = print(io, "{} (empty domain)")

isempty(d::EmptySpace) = true

isopenset(d::EmptySpace) = true
isclosedset(d::EmptySpace) = true

boundary(d::EmptySpace) = d
interior(d::EmptySpace) = d
closure(d::EmptySpace) = d
boundingbox(d::EmptySpace) = d

# Arithmetic operations

issubset1(d1::EmptySpace, d2) = true
issubset2(d1, d2::EmptySpace) = isempty(d1)

# setdiffdomain2(d1, d2::EmptySpace) = d2

map_domain(map::Map{T}, d::EmptySpace{T}) where {T} = d
mapped_domain(map::Map, d::EmptySpace) = EmptySpace{codomaintype(map)}()

==(d1::EmptySpace, d2::EmptySpace) = true
isequal1(d1::EmptySpace, d2) = isempty(d2)
isequal2(d1, d2::EmptySpace) = isempty(d1)

"The full space of elements of type `T`."
struct FullSpace{T} <: Domain{T} end

const AnyFullSpace = FullSpace{Any}

FullSpace() = FullSpace{Float64}()
FullSpace(d) = FullSpace{eltype(d)}()

fullspace(d) = fullspace(eltype(d))
fullspace(::Type{T}) where {T} = FullSpace{T}()

isfullspace(d::FullSpace) = true
isfullspace(d::Domain) = false

similardomain(::FullSpace, ::Type{T}) where {T} = FullSpace{T}()

euclideanspace(n::Val{N}) where {N} = euclideanspace(n, Float64)
euclideanspace(::Val{N}, ::Type{T}) where {N,T} = FullSpace{SVector{N,T}}()

indomain(x::T, d::FullSpace{T}) where {T} = true
approx_indomain(x, d::FullSpace, tolerance) = in(x, d)

show(io::IO, d::FullSpace) = print(io, "{x} (full space)")

# We choose the origin as a point in the full space
point_in_domain(d::FullSpace) = zero(eltype(d))

isempty(::FullSpace) = false

isopenset(d::FullSpace) = true
isclosedset(d::FullSpace) = true

boundary(d::FullSpace{T}) where {T} = EmptySpace{T}()
interior(d::FullSpace) = d
closure(d::FullSpace) = d

# Arithmetic operations

issubset1(d1::FullSpace, d2) = isfullspace(d2)
issubset2(d1, d2::FullSpace) = true

map_domain(m::AbstractAffineMap{T}, d::FullSpace{T}) where {T} = d

==(d1::FullSpace, d2::FullSpace) = true
isequal1(d1::FullSpace, d2) = isfullspace(d2)
isequal2(d1, d2::FullSpace) = isfullspace(d1)

convert(::Type{Domain}, ::Type{T}) where T = FullSpace{T}()
convert(::Type{Domain{S}}, ::Type{T}) where {T,S} = convert(Domain{S}, convert(Domain, T))

infimum(d::FullSpace{T}) where {T} = typemin(T)
supremum(d::FullSpace{T}) where {T} = typemax(T)

# Some convenient complete domains

"The set of integers of type Int (ℤ = \\BbbZ)."
const ℤ = FullSpace{Int}()
"The set of rational numbers of type Rational{Int} (ℚ = \\BbbQ)."
const ℚ = FullSpace{Rational{Int}}()
"The set of reals of type Float64 (ℝ = \\BbbR)."
const ℝ = FullSpace{Float64}()
"The complex plane with Float64 real and imaginar parts (ℂ = \\BbbC)."
const ℂ = FullSpace{Complex{Float64}}()

"The space ℝ^1"
const ℝ1 = FullSpace{SVector{1,Float64}}()
"The space ℝ^2"
const ℝ2 = FullSpace{SVector{2,Float64}}()
"The space ℝ^3"
const ℝ3 = FullSpace{SVector{3,Float64}}()
"The space ℝ^4"
const ℝ4 = FullSpace{SVector{4,Float64}}()