text:
function applycustomics!(system::CVSystem)
    # custom peripheral volumes for better initial blood distribution
    system.branches.term[12].P[1,2] = 73*mmHgToPa; # splenic artery
    system.branches.term[12].V[1,2] = (system.branches.term[12].P[1,2]*
        system.branches.term[12].C[2] + system.branches.term[12].V0[2]);
    system.branches.term[12].P[1,3] = 37.2*mmHgToPa;
    system.branches.term[12].V[1,3] = (system.branches.term[12].P[1,3]*
        system.branches.term[12].C[3] + system.branches.term[12].V0[3]);
    system.branches.term[12].P[1,4] = 5.66*mmHgToPa;
    system.branches.term[12].V[1,4] = (system.branches.term[12].P[1,4]*
        system.branches.term[12].C[4] + system.branches.term[12].V0[4]);
    system.branches.term[12].P[1,5] = 5.2*mmHgToPa;
    system.branches.term[12].V[1,5] = (system.branches.term[12].P[1,5]*
        system.branches.term[12].C[5] + system.branches.term[12].V0[5]);
    system.branches.term[13].P[1,2] = 62*mmHgToPa; # gastric artery
    system.branches.term[13].V[1,2] = (system.branches.term[13].P[1,2]*
        system.branches.term[13].C[2] + system.branches.term[13].V0[2]);
    system.branches.term[13].P[1,3] = 30.5*mmHgToPa;
    system.branches.term[13].V[1,3] = (system.branches.term[13].P[1,3]*
        system.branches.term[13].C[3] + system.branches.term[13].V0[3]);
    system.branches.term[13].P[1,4] = 5.46*mmHgToPa;
    system.branches.term[13].V[1,4] = (system.branches.term[13].P[1,4]*
        system.branches.term[13].C[4] + system.branches.term[13].V0[4]);
    system.branches.term[13].P[1,5] = 5.09*mmHgToPa;
    system.branches.term[13].V[1,5] = (system.branches.term[13].P[1,5]*
        system.branches.term[13].C[5] + system.branches.term[13].V0[5]);
    system.branches.term[14].P[1,2] = 69*mmHgToPa; # hepatic artery
    system.branches.term[14].V[1,2] = (system.branches.term[14].P[1,2]*
        system.branches.term[14].C[2] + system.branches.term[14].V0[2]);
    system.branches.term[14].P[1,3] = 34.2*mmHgToPa;
    system.branches.term[14].V[1,3] = (system.branches.term[14].P[1,3]*
        system.branches.term[14].C[3] + system.branches.term[14].V0[3]);
    system.branches.term[14].P[1,4] = 5.21*mmHgToPa;
    system.branches.term[14].V[1,4] = (system.branches.term[14].P[1,4]*
        system.branches.term[14].C[4] + system.branches.term[14].V0[4]);
    system.branches.term[14].P[1,5] = 4.78*mmHgToPa;
    system.branches.term[14].V[1,5] = (system.branches.term[14].P[1,5]*
        system.branches.term[14].C[5] + system.branches.term[14].V0[5]);
    system.branches.term[15].P[1,2] = 74*mmHgToPa; # superior mesenteric artery
    system.branches.term[15].V[1,2] = (system.branches.term[15].P[1,2]*
        system.branches.term[15].C[2] + system.branches.term[15].V0[2]);
    system.branches.term[15].P[1,3] = 39*mmHgToPa;
    system.branches.term[15].V[1,3] = (system.branches.term[15].P[1,3]*
        system.branches.term[15].C[3] + system.branches.term[15].V0[3]);
    system.branches.term[15].P[1,4] = 5.62*mmHgToPa;
    system.branches.term[15].V[1,4] = (system.branches.term[15].P[1,4]*
        system.branches.term[15].C[4] + system.branches.term[15].V0[4]);
    system.branches.term[15].P[1,5] = 5.13*mmHgToPa;
    system.branches.term[15].V[1,5] = (system.branches.term[15].P[1,5]*
        system.branches.term[15].C[5] + system.branches.term[15].V0[5]);
    system.branches.term[16].P[1,2] = 78*mmHgToPa;
    system.branches.term[16].V[1,2] = (system.branches.term[16].P[1,2]*
        system.branches.term[16].C[2] + system.branches.term[16].V0[2]);
    system.branches.term[16].P[1,3] = 38.5*mmHgToPa;
    system.branches.term[16].V[1,3] = (system.branches.term[16].P[1,3]*
        system.branches.term[16].C[3] + system.branches.term[16].V0[3]);
    system.branches.term[18].P[1,2] = 78*mmHgToPa;
    system.branches.term[18].V[1,2] = (system.branches.term[18].P[1,2]*
        system.branches.term[18].C[2] + system.branches.term[18].V0[2]);
    system.branches.term[18].P[1,3] = 38.5*mmHgToPa;
    system.branches.term[18].V[1,3] = (system.branches.term[18].P[1,3]*
        system.branches.term[18].C[3] + system.branches.term[18].V0[3]);
    system.branches.term[21].P[1,2] = 57.5*mmHgToPa; # inferior mesenteric
    system.branches.term[21].V[1,2] = (system.branches.term[21].P[1,2]*
        system.branches.term[21].C[2] + system.branches.term[21].V0[2]);
    system.branches.term[21].P[1,3] = 28.4*mmHgToPa;
    system.branches.term[21].V[1,3] = (system.branches.term[21].P[1,3]*
        system.branches.term[21].C[3] + system.branches.term[21].V0[3]);
    system.branches.term[21].P[1,4] = 5.62*mmHgToPa; # inferior mesenteric
    system.branches.term[21].V[1,4] = (system.branches.term[21].P[1,4]*
        system.branches.term[21].C[4] + system.branches.term[21].V0[4]);
    system.branches.term[21].P[1,5] = 5.3*mmHgToPa;
    system.branches.term[21].V[1,5] = (system.branches.term[21].P[1,5]*
        system.branches.term[21].C[5] + system.branches.term[21].V0[5]);
    system.branches.term[23].P[1,2] = 60*mmHgToPa;
    system.branches.term[23].V[1,2] = (system.branches.term[23].P[1,2]*
        system.branches.term[23].C[2] + system.branches.term[23].V0[2]);
    system.branches.term[23].P[1,3] = 30*mmHgToPa;
    system.branches.term[23].V[1,3] = (system.branches.term[23].P[1,3]*
        system.branches.term[23].C[3] + system.branches.term[23].V0[3]);
    system.branches.term[27].P[1,2] = 61*mmHgToPa;
    system.branches.term[27].V[1,2] = (system.branches.term[27].P[1,2]*
        system.branches.term[27].C[2] + system.branches.term[27].V0[2]);
    system.branches.term[27].P[1,3] = 30.5*mmHgToPa;
    system.branches.term[27].V[1,3] = (system.branches.term[27].P[1,3]*
        system.branches.term[27].C[3] + system.branches.term[27].V0[3]);
    system.branches.term[33].P[1,2] = 42*mmHgToPa;
    system.branches.term[33].V[1,2] = (system.branches.term[33].P[1,2]*
        system.branches.term[33].C[2] + system.branches.term[33].V0[2]);
    system.branches.term[33].P[1,3] = 21*mmHgToPa;
    system.branches.term[33].V[1,3] = (system.branches.term[33].P[1,3]*
        system.branches.term[33].C[3] + system.branches.term[33].V0[3]);
    system.branches.term[35].P[1,2] = 16.4*mmHgToPa;
    system.branches.term[35].V[1,2] = (system.branches.term[35].P[1,2]*
        system.branches.term[35].C[2] + system.branches.term[35].V0[2]);
    system.branches.term[35].P[1,3] = 9.7*mmHgToPa;
    system.branches.term[35].V[1,3] = (system.branches.term[35].P[1,3]*
        system.branches.term[35].C[3] + system.branches.term[35].V0[3]);
    system.branches.term[37].P[1,2] = 60*mmHgToPa;
    system.branches.term[37].V[1,2] = (system.branches.term[37].P[1,2]*
        system.branches.term[37].C[2] + system.branches.term[37].V0[2]);
    system.branches.term[37].P[1,3] = 29*mmHgToPa;
    system.branches.term[37].V[1,3] = (system.branches.term[37].P[1,3]*
        system.branches.term[37].C[3] + system.branches.term[37].V0[3]);
    system.branches.term[39].P[1,2] = 62*mmHgToPa;
    system.branches.term[39].V[1,2] = (system.branches.term[39].P[1,2]*
        system.branches.term[39].C[2] + system.branches.term[39].V0[2]);
    system.branches.term[39].P[1,3] = 30.5*mmHgToPa;
    system.branches.term[39].V[1,3] = (system.branches.term[39].P[1,3]*
        system.branches.term[39].C[3] + system.branches.term[39].V0[3]);
    system.branches.term[47].P[1,2] = 16.3*mmHgToPa;
    system.branches.term[47].V[1,2] = (system.branches.term[47].P[1,2]*
        system.branches.term[47].C[2] + system.branches.term[47].V0[2]);
    system.branches.term[47].P[1,3] = 9.7*mmHgToPa;
    system.branches.term[47].V[1,3] = (system.branches.term[47].P[1,3]*
        system.branches.term[47].C[3] + system.branches.term[47].V0[3]);
    system.branches.term[48].P[1,2] = 42*mmHgToPa;
    system.branches.term[48].V[1,2] = (system.branches.term[48].P[1,2]*
        system.branches.term[48].C[2] + system.branches.term[48].V0[2]);
    system.branches.term[48].P[1,3] = 20.8*mmHgToPa;
    system.branches.term[48].V[1,3] = (system.branches.term[48].P[1,3]*
        system.branches.term[48].C[3] + system.branches.term[48].V0[3]);
    system.branches.term[90].P[1,2] = 16.6*mmHgToPa;
    system.branches.term[90].V[1,2] = (system.branches.term[90].P[1,2]*
        system.branches.term[90].C[2] + system.branches.term[90].V0[2]);
    system.branches.term[90].P[1,3] = 9.9*mmHgToPa;
    system.branches.term[90].V[1,3] = (system.branches.term[90].P[1,3]*
        system.branches.term[90].C[3] + system.branches.term[90].V0[3]);
    system.branches.term[54].P[1,2] = 16.4*mmHgToPa;
    system.branches.term[54].V[1,2] = (system.branches.term[54].P[1,2]*
        system.branches.term[54].C[2] + system.branches.term[54].V0[2]);
    system.branches.term[54].P[1,3] = 9.75*mmHgToPa;
    system.branches.term[54].V[1,3] = (system.branches.term[54].P[1,3]*
        system.branches.term[54].C[3] + system.branches.term[54].V0[3]);
    system.branches.term[55].P[1,2] = 16.4*mmHgToPa;
    system.branches.term[55].V[1,2] = (system.branches.term[55].P[1,2]*
        system.branches.term[55].C[2] + system.branches.term[55].V0[2]);
    system.branches.term[55].P[1,3] = 9.75*mmHgToPa;
    system.branches.term[55].V[1,3] = (system.branches.term[55].P[1,3]*
        system.branches.term[55].C[3] + system.branches.term[55].V0[3]);
    system.branches.term[57].P[1,2] = 60*mmHgToPa;
    system.branches.term[57].V[1,2] = (system.branches.term[57].P[1,2]*
        system.branches.term[57].C[2] + system.branches.term[57].V0[2]);
    system.branches.term[57].P[1,3] = 29.2*mmHgToPa;
    system.branches.term[57].V[1,3] = (system.branches.term[57].P[1,3]*
        system.branches.term[57].C[3] + system.branches.term[57].V0[3]);
    system.branches.term[58].P[1,2] = 57*mmHgToPa;
    system.branches.term[58].V[1,2] = (system.branches.term[58].P[1,2]*
        system.branches.term[58].C[2] + system.branches.term[58].V0[2]);
    system.branches.term[58].P[1,3] = 27.7*mmHgToPa;
    system.branches.term[58].V[1,3] = (system.branches.term[58].P[1,3]*
        system.branches.term[58].C[3] + system.branches.term[58].V0[3]);
    system.branches.term[67].P[1,2] = 44*mmHgToPa;
    system.branches.term[67].V[1,2] = (system.branches.term[67].P[1,2]*
        system.branches.term[67].C[2] + system.branches.term[67].V0[2]);
    system.branches.term[67].P[1,3] = 21.7*mmHgToPa;
    system.branches.term[67].V[1,3] = (system.branches.term[67].P[1,3]*
        system.branches.term[67].C[3] + system.branches.term[67].V0[3]);
    system.branches.term[68].P[1,2] = 18.2*mmHgToPa;
    system.branches.term[68].V[1,2] = (system.branches.term[68].P[1,2]*
        system.branches.term[68].C[2] + system.branches.term[68].V0[2]);
    system.branches.term[68].P[1,3] = 10.5*mmHgToPa;
    system.branches.term[68].V[1,3] = (system.branches.term[68].P[1,3]*
        system.branches.term[68].C[3] + system.branches.term[68].V0[3]);
    system.branches.term[70].P[1,2] = 54*mmHgToPa;
    system.branches.term[70].V[1,2] = (system.branches.term[70].P[1,2]*
        system.branches.term[70].C[2] + system.branches.term[70].V0[2]);
    system.branches.term[70].P[1,3] = 26.5*mmHgToPa;
    system.branches.term[70].V[1,3] = (system.branches.term[70].P[1,3]*
        system.branches.term[70].C[3] + system.branches.term[70].V0[3]);
    system.branches.term[73].P[1,2] = 58.8*mmHgToPa;
    system.branches.term[73].V[1,2] = (system.branches.term[73].P[1,2]*
        system.branches.term[73].C[2] + system.branches.term[73].V0[2]);
    system.branches.term[73].P[1,3] = 28.8*mmHgToPa;
    system.branches.term[73].V[1,3] = (system.branches.term[73].P[1,3]*
        system.branches.term[73].C[3] + system.branches.term[73].V0[3]);
    system.branches.term[83].P[1,2] = 17.9*mmHgToPa;
    system.branches.term[83].V[1,2] = (system.branches.term[83].P[1,2]*
        system.branches.term[83].C[2] + system.branches.term[83].V0[2]);
    system.branches.term[83].P[1,3] = 10.3*mmHgToPa;
    system.branches.term[83].V[1,3] = (system.branches.term[83].P[1,3]*
        system.branches.term[83].C[3] + system.branches.term[83].V0[3]);
    system.branches.term[84].P[1,2] = 54*mmHgToPa;
    system.branches.term[84].V[1,2] = (system.branches.term[84].P[1,2]*
        system.branches.term[84].C[2] + system.branches.term[84].V0[2]);
    system.branches.term[84].P[1,3] = 26.4*mmHgToPa;
    system.branches.term[84].V[1,3] = (system.branches.term[84].P[1,3]*
        system.branches.term[84].C[3] + system.branches.term[84].V0[3]);
    system.branches.term[85].P[1,2] = 43.3*mmHgToPa;
    system.branches.term[85].V[1,2] = (system.branches.term[85].P[1,2]*
        system.branches.term[85].C[2] + system.branches.term[85].V0[2]);
    system.branches.term[85].P[1,3] = 21.5*mmHgToPa;
    system.branches.term[85].V[1,3] = (system.branches.term[85].P[1,3]*
        system.branches.term[85].C[3] + system.branches.term[85].V0[3]);
    system.branches.term[91].P[1,2] = 16*mmHgToPa;
    system.branches.term[91].V[1,2] = (system.branches.term[91].P[1,2]*
        system.branches.term[91].C[2] + system.branches.term[91].V0[2]);
    system.branches.term[91].P[1,3] = 9.52*mmHgToPa;
    system.branches.term[91].V[1,3] = (system.branches.term[91].P[1,3]*
        system.branches.term[91].C[3] + system.branches.term[91].V0[3]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
include("leap_frog_definitions.jl")
max_hit=10000
barrier=20

condition_max_hits(u,t,integrator)= u[6]>max_hit
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS(u,t,integrator) # Event when event_f(u,t) == 0
   u[4]
end
function affect_update_iterator!(integrator)
    integrator.u[6]=integrator.u[6]+1
end

callback_max_hits=DiscreteCallback(condition_max_hits,affect_stop!)
callback_hits_PSS=ContinuousCallback(condition_hits_PSS, affect_update_iterator!,nothing)
cb=CallbackSet(callback_hits_PSS, callback_max_hits)

function PSS_function(Q2,P2, H,  t_end)
    #for a given Q,P,H with X=0
    Q1=Q1_find_dimer(Q2,P2,H)
    if ~isempty(Q1)
        # println(Q1)
        q0,p0=[zeros(3) for i in 1:2]
        q0[1]=Q1[1]
        q0[2]=Q2
        q0[3]=0
        p0[1]=0
        p0[2]=P2
        p0[3]=0
       #constructor for ODE
        prob= HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, RK4(),maxiters=1e20, reltol=1e-6,abstol=1e-6,callback=cb,save_start=true,save_end=true,save_everystep=false)
        #output Q, P and dH
        return sol[:,2:end-1][2,:],sol[:,2:end-1][5,:]
    else
    #need to return 3 values, dH=1 flags that there is no Y
        return 0, 0
    end
end

ext: jl


text:
module BibInternal

"""
Abstract entry supertype.
"""
abstract type AbstractEntry end

# Includes
include("constant.jl")
include("utilities.jl")
include("bibtex.jl")
include("entry.jl")

end # module


text:
@testset begin
    config = _default_config()

    project_contents = """
        [deps]
        Example = "7876af07-990d-54b4-ab0e-23690620f79a"
    """

    manifest_contents = project_contents_to_manifest_contents(
        config;
        project_file_contents = project_contents,
    )

    projs_and_mnfsts = [
        Manifest(FileContents(manifest_contents)),
    ]

    @test warmup_cache(config, projs_and_mnfsts) isa Nothing
end


text:
const TensorSubGrid = ProductGrid{<:NTuple{N,GRID}} where N where {GRID<:AbstractSubGrid}
function tensorproductbitarray(vectors::Union{BitVector,Vector{Bool}}...)
	N = length(vectors)
	R = falses(map(length, vectors))
	for i in CartesianIndices(size(R))
		R[i] = reduce(&,map(k->vectors[k][i.I[k]],1:length(vectors)))
	end
	R
end

mask(grid::TensorSubGrid) = tensorproductbitarray(map(mask, elements(grid))...)
subindices(grid::TensorSubGrid) = findall(mask(grid))
supergrid(grid::TensorSubGrid) = ProductGrid(map(supergrid, elements(grid))...)
issubindex(i, g::TensorSubGrid) = all(map(issubindex, i, elements(g)))
issubindex(i::CartesianIndex, g::TensorSubGrid) = issubindex(i.I, g)


text:
abstract type StopCriteria end
struct NoStoppingCriteria <: StopCriteria end
mutable struct MeanVelocityStoppingCriteria{T <: Real} <: StopCriteria
    old_mean_velocity::T
    tolerance::T
    problem::FluidFlowProblem
end
StopCriteria(problem::FluidFlowProblem) = NoStoppingCriteria()
StopCriteria(problem::PoiseuilleFlow) = MeanVelocityStoppingCriteria(0.0, 1e-12, problem)
StopCriteria(problem::CouetteFlow) = MeanVelocityStoppingCriteria(0.0, 1e-7, problem)
StopCriteria(problem::LidDrivenCavityFlow) =
    MeanVelocityStoppingCriteria(0.0, 1e-5, problem)
StopCriteria(problem::DecayingShearFlow) =
    problem.static ? MeanVelocityStoppingCriteria(0.0, 1e-8, problem) : NoStoppingCriteria()

should_stop!(::StopCriteria, q, f_in) = false
function should_stop!(stop_criteria::MeanVelocityStoppingCriteria{T}, q, f_in) where {T}
    f = Array{T}(undef, size(f_in, 3))
    u = zeros(dimension(q))
    Nx = size(f_in, 1)
    Ny = size(f_in, 2)

    divide_by = 0
    u_mean = 0.0

    @inbounds for x_idx in 1:Nx, y_idx in 1:Ny
        divide_by += 1

        # Calculated
        @inbounds for f_idx in 1:size(f_in, 3)
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end
        ρ = density(q, f)
        velocity!(q, f, ρ, u)

        u_mean += u[1]
    end

    u_mean /= divide_by

    converged = abs(u_mean / stop_criteria.old_mean_velocity - 1)

    if (converged < stop_criteria.tolerance)
        return true
    end

    if (isnan(u_mean))
        @warn "nan in velocity profile"
        return true
    end

    stop_criteria.old_mean_velocity = u_mean

    return false
end

struct VelocityConvergenceStoppingCriteria{T <: Real, VT <: AbstractVector{T}} <:
       StopCriteria
    old_velocity::Array{VT, 2}
    tolerance::T
    problem::FluidFlowProblem
end
function VelocityConvergenceStoppingCriteria(tolerance::T, problem) where {T <: Real}
    return VelocityConvergenceStoppingCriteria(
        [zeros(T, 2) for x_idx in 1:(problem.NX), y_idx in 1:(problem.NY)],
        tolerance,
        problem,
    )
end

function should_stop!(
    stop_criteria::VelocityConvergenceStoppingCriteria{T},
    q,
    f_in,
) where {T}
    f = Array{T}(undef, size(f_in, 3))
    u = zeros(dimension(q))
    Nx = size(f_in, 1)
    Ny = size(f_in, 2)

    old_velocity_norm = 0.0 # norm(stop_criteria.old_velocity)
    error = 0.0

    @inbounds for x_idx in 1:Nx, y_idx in 1:Ny
        # Calculated
        @inbounds for f_idx in 1:size(f_in, 3)
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end

        ρ = density(q, f)
        velocity!(q, f, ρ, u)

        u_old = stop_criteria.old_velocity[x_idx, y_idx]

        error += (((u[1] - u_old[1])^2 + (u[2] - u_old[2])^2))

        old_velocity_norm += u_old[1]^2 + u_old[2]^2
        stop_criteria.old_velocity[x_idx, y_idx] = copy(u)
    end

    converged = sqrt(error) / old_velocity_norm

    if (converged < stop_criteria.tolerance)
        return true
    end

    if (isnan(converged))
        @warn "nan in velocity profile"
        return true
    end

    # stop_criteria.old_mean_velocity = u_mean

    return false
end


text:
T = Complex{Float64}
fs = 100.
N = 90000
expect_output_type = Float64

# Test calculate_psd function
@testset "calculate PSD" begin
    data = rand(T, N)
    psd = calculate_psd(data, fs)
    @test size(psd, 1) == N
    @test eltype(psd) == expect_output_type
end


text:
import JSON2, LibGit2
const ROOT_PATH = (@__DIR__) * "/.."
const META_FILENAME = ROOT_PATH * "/component_sources.json"
sources = open(META_FILENAME, "r") do f
    JSON2.read(f)
end
for source in sources
    if isdir(ROOT_PATH*"/"*source.path)
        rm(ROOT_PATH*"/"*source.path, force = true, recursive = true)
    end
    repo = LibGit2.clone(source.git, ROOT_PATH*"/"*source.path)
    LibGit2.checkout!(repo, source.hash)
end


text:
# Code used to generate SDP problems with a quadratic objective save in Testproblems/SDPQuad with an accurante
# solution from MOSEK
# The problem data is generated based on the Optimality conditions:
# Ax*+s*=b, s in S+, Px*+q+A'y*=0, y in S+

# Create different objective functions
# 1) P = rand pos def
# 2) P = I
# 3) P = 0

# A) Create problems like Ax+s=b, mat(S) in S+
# B) Create problems like <Ai,X> = bi, X in S+

# workspace()
# include("./Helper.jl")
# include("../../OSSDP/Code/src/Solver.jl")

# using JuMP, Mosek, JLD, HelperFunctions, OSSDP, Base.Test

rng = MersenneTwister(12345)
dirPath = "../DataFiles/Julia/SDPQuad/"
!ispath(dirPath) && mkdir(dirPath)
nn = 45

  # @testset "Random SDP Mosek vs. QOCS" begin

 for iii =1:1:nn
  # choose size of problem
  n = rand(rng,10:100)
  r = rand(rng,2:15)
  m = r^2
  A = spzeros(m,n)
  # construct A with several symmetric Ai that are vectorized and make the cols of A
  for kkk=1:n
    Ai = sprandn(rng,r,r,0.4)
    Ai = 0.5*(Ai+Ai')
    A[:,kkk] = vec(Ai)
  end
  xtrue = randn(rng,n,1)

  # create primal feasible point -> b
  Strue = generatePosDefMatrix(r,rng)
  strue = vec(Strue)
  b = A*xtrue + strue

  # create dual feasible point -> P,q
  # consider three different types of the quadratic cost
  if iii <= 15
     P = generatePosDefMatrix(n,rng)
     Phalf = P^(1/2)
  elseif iii <=30
     P = speye(n)
     Phalf = P
  else
   P = Phalf = spzeros(n,n)
  end

  Ytrue = generatePosDefMatrix(r,rng)
  ytrue = vec(Ytrue)
  q = (-P*xtrue -  A'*ytrue)[:]
  ra = 0.
  Kf = 0
  Kl = 0
  Kq = []
  Ks = [r^2]

  iii<=30 && (PhalfInvQ=Phalf\q)

  # K = Cone(0,0,[],[r^2])
  # setOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 0,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=false)
  # setON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 10,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=false)
  # setAdaptiveOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 0,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=true)
  # setAdaptiveON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 10,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=true)
  # res1,nothing = OSSDP.solve(P,q,A,b,K,setOFF);
  # print("\n.")
  # res2,nothing = OSSDP.solve(P,q,A,b,K,setON);
  # print(".")
  # res3,nothing = OSSDP.solve(P,q,A,b,K,setAdaptiveOFF);
  # print(".")
  # res4,nothing = OSSDP.solve(P,q,A,b,K,setAdaptiveON);
  # print(".")

  # println(">>>QOCS: Problem $(iii), res1: cost: $(res1.cost), status: $(res1.status)")
  # println(">>>QOCS: Problem $(iii), res2: cost: $(res2.cost), status: $(res2.status)")
  # println(">>>QOCS: Problem $(iii), res3: cost: $(res3.cost), status: $(res3.status)")
  # println(">>>QOCS: Problem $(iii), res4: cost: $(res4.cost), status: $(res4.status)")

  # solve accurately once with mosek s
  model = Model(solver=MosekSolver())
  @variable(model, x[1:n])
  @variable(model, S[1:r,1:r],SDP)
  # if P!= 0 add soc-constraint to model
  if P == spzeros(n,n)
    @objective(model, Min, q'*x)
  else
    # in case of quadratic objective the following problem is solved:
    # min x'Px + 2q'x
    # s.t. Ax+s=b, s in S+
    # instead solve equivalent problem
    # min t
    # s.t. ||Phalf*x+inv(Phalf)*q|| <= t
    # Ax+s=b, s in S+
    @variable(model, t)
    @objective(model, Min, t)
    @constraint(model,norm(Phalf*x+PhalfInvQ) <= t)
  end
  s = vec(S)
  @constraint(model,A*x+s.==b)
  status = JuMP.solve(model)
  objTrue = getobjectivevalue(model)

  (P!=spzeros(n,n)) && (objTrue = 0.5*( objTrue^2-q'*(P\q) ))
  solTrue = getvalue(x)

  # @test abs(objTrue-res1.cost) < 1e-2
  # @test norm(solTrue-res1.x,Inf) < 1e-2

  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end

  fn = "SDPQuad"*nr*".jld"
  problemType = "Random SDP with quadratic Objective"
  problemName = "RandomSDP"*nr

  JLD.save(dirPath*fn,"n",n,"m",m,"A",A,"b",b,"P",P,"q",q,"r",ra,"objTrue",objTrue,"solTrue",solTrue,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks)
  println("$(iii)/$(nn) completed!")
end
  # end


text:
"""
Struct with the following fields
- df
- horizons
- models
- transformfunc
- grwthoffset
- formula
- regstyle
- ynames
"""
@with_kw struct LocalProjection
    df                  ::DataFrame
    horizons            ::Array{Int,1} = [1]
    models              ::Array{FixedEffectModel,1} = Array{FixedEffectModel,1}(undef, length(collect(horizons)))
    transformfunc       ::Symbol = :identity
    grwthoffset         ::Int = 0
    formula             ::FormulaTerm
    regstyle            ::Symbol = :rates
    ynames              ::Array{Symbol,1} = Array{Symbol,1}(undef, length(collect(horizons)))
end

#######################################################################
# LocalProjection functions
#######################################################################
"""
    makedepvars(lp::LocalProjection)
Takes the lhs variable from `lp.formula` and makes the cumlative growth  variables
for each horizon specified by `horizons` for the y-variable. This updates the
dataframe stored in `lp`. Also updates the `ynames` array with the name of the
created variables
"""
function makedepvars(lp::LocalProjection)
    if lp.regstyle == :rates
        FΔ!(lp.df, Symbol(lp.formula.lhs), lp.horizons, getfield(Main, lp.transformfunc), lp.grwthoffset)
        lp.ynames .= [Symbol("FΔ", i, lp.formula.lhs) for i in lp.horizons]
    elseif lp.regstyle == :levels
        F!(lp.df, Symbol(lp.formula.lhs), lp.horizons, getfield(Main, lp.transformfunc))
        lp.ynames .= [Symbol("F", i, lp.formula.lhs) for i in lp.horizons]
    end
end

"""
    runregs(lp::LocalProjection)
Estimates the regression for each horizon using `FixedEffectModels.reg(lp.df, Term(lp.ynames[index]) ~ lp.formula.rhs)`
and store the estimated models in `lp.models`. The estimations are done lazely, i.e. each regression starts
from scratch without saving intermediate computations common across each regression
"""
function runregs(lp::LocalProjection)
    for (index, value) in enumerate(lp.horizons)
        form = Term(lp.ynames[index]) ~ lp.formula.rhs
        lp.models[index] = reg(lp.df, form)
    end
end

"""
    coef(lp::LocalProjection, xname::Symbol)
Returns a vector of coefficients for the variable `xname` for each regression. `coef(lp, xname)[i] `
corresponds to the coefficient for the horizon `lp.horizons[i]`.
"""
function coef(lp::LocalProjection, xname::Symbol)
    idx = findfirst(Symbol.(coefnames(lp.models[1])) .== xname)
    β = [m.coef[idx] for m in lp.models]
    return β
end

"""
    se(lp::LocalProjection, xname::Symbol)
Returns a vector of standard errors for each coefficient for the variable `xname` for each regression.
`se(lp, xname)[i] ` corresponds to the se for the horizon `lp.horizons[i]`.
"""
function se(lp::LocalProjection, xname::Symbol)
    idx = findfirst(Symbol.(coefnames(lp.models[1])) .== xname)
    se = [sqrt(vcov(m)[idx, idx]) for m in lp.models]
    return se
end

"""
    coef_se(lp::LocalProjection, xname::Symbol)
Returns a dataframe with the columns `[horizons, β, se]` aligning the
estimated coefficient, standard error and horizon that these correspond to.
Constructed from
`DataFrame(horizons = lp.horizons, β = coef(lp, xname), se = se(lp, xname))`
"""
function coef_se(lp::LocalProjection, xname::Symbol)
    return DataFrame(horizons = lp.horizons, β = coef(lp, xname), se = se(lp, xname))
end


text:
# Interfaces
function cor(x::AbstractArray, cosine::Bool, verbose::Bool)::Void
    """
    Pearson correlation of a matrix with itself.
    This one uses matrix multiplication in BLAS to speed up calculation when there are no NA's
    and uses threading to speed up the rest of the calculation.
    """
    double * result, int *nNA
    nr, nc = size(x)
    nNA_ext = 0

    # This matrix will hold preprocessed entries that can be simply multiplied together to get the
    # numerator
    multMat = zeros(nr, nc)
    # Number of NA entries in each column
    nNAentries = zeros(nc)
    # Flag indicating whether the mean of each column is NA
    NAmean = falses(nc)

    # Decide how many threads to use
    nt = useNThreads( nc*nc, *nThreads)

    # Put the general data of the correlation calculation into a structure that can be passed on to
    # threads.
    cor1ThreadData thrdInfo[MxThreads];  # initailized

    # Column preparation (calculation of the matrix to be multiplied) in a threaded form.
    cptd = Array{colPrepThreadData}(MxThreads)
    thr = Array{pthread_t}(MxThreads)
    status = Array{Int64}(MxThreads)  # normal: 0
    pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
    progressCounter pc; pc.i = 0; pc.n = nc;

    # calculate threadPrepColCor using data cptd in each thr

    result = multMat' * multMat  # set NA means to zero

    # Here I need to recalculate results that have NA's in them.

    nNA = (int) nNA_ext;
    free(NAmean);
    free(nNAentries);
    free(multMat);
end

function bicor(x::AbstractArray, maxPOutliers::Float64, fallback::Int64, cosine::Bool,
               verbose::Bool)::Void
    """
    bicorrelation of a matrix with itself.
    This one uses matrix multiplication in BLAS to speed up calculation when there are no NA's
    and is threaded to speed up the rest of the calculation.
    """
    double * result, int *nNA
    nr, nc = size(x)
    nNA = 0
    nNA_ext = 0

    # This matrix will hold preprocessed entries that can be simply multiplied together to get the
    # numerator
    multMat = zeros(nr, nc)
    # Number of NA entries in each column
    nNAentries = zeros(nc)
    # Flag indicating whether the mean of each column is NA
    NAmed = falses(nc)

    # Decide how many threads to use
    int nt = useNThreads( nc*nc, *nThreads);

    double * aux[MxThreads];

    for (int t=0; t < nt; t++)
    {
       if ( (aux[t] = (double *) malloc(6*nr * sizeof(double)))==NULL)
       {
         *err = 1;
         Rprintf("cor1: memory allocation error. The needed block is very small... suspicious.\n");
         for (int tt = t-1; tt>=0; tt--) free(aux[tt]);
         free(NAmed); free(nNAentries); free(multMat);
         return;
       }
    }

    # Put the general data of the correlation calculation into a structure that can be passed on to
    # threads.

    cor1ThreadData thrdInfo[MxThreads];

    # Column preparation (calculation of the matrix to be multiplied) in a threaded form.
    cptd = Array{colPrepThreadData}(MxThreads)
    thr = Array{pthread_t}(MxThreads)
    status = Array{Int64}(MxThreads)  # normal: 0
    pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
    progressCounter pc; pc.i = 0; pc.n = nc;

    # Rprintf("Preparing columns...\n");
    # calculate threadPrepColBicor using data cptd in each thr

    int pearson = 0;

    if (*fallback==3)
    {
      for (int t=0; t<nt; t++) if (thrdInfo[t].zeroMAD > 0)
      {
        pearson = 1;
        if (*verbose)
          Rprintf("Warning in bicor(x): Thread %d (of %d) reported zero MAD in column %d. %s",
                  t, nt, thrdInfo[t].zeroMAD, "Switching to Pearson correlation.\n");
      }
      if (pearson==1) # Re-do all column preparations using Pearson preparation.
      {
        # Set fallback to 4 for slow calculations below.
        for (int t = 0; t < nt; t++) thrdInfo[t].fallback = 4;

        pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER;
        pc.i = 0; pc.n = nc;
        # calculate threadPrepColCor using data cptd in each thr with thrdInfo
      }
    }

    result = multMat' * multMat  # lower triangle only

    # Here I need to recalculate results that have NA's in them.

    # Symmetrize the result and set all rows and columns with NA means to zero

    for (int t=nt-1; t >= 0; t--) free(aux[t]);

    *nNA = (int) nNA_ext;

    free(NAmed);
    free(nNAentries);
    free(multMat);
end

function bicor(x::AbstractArray, y::AbstractArray, robustX::Bool, robustY::Bool,
               maxPOutliers::Float64, fallback::Int64, cosineX::Bool, cosineY::Bool,
               verbose::Bool)::Void
    """
    Two-variable bicorrelation. Basically the same as bicor1, just must calculate the whole matrix.
    If robustX,Y is zero, the corresponding variable will be treated as in pearson correlation.
    """
    double * result, int *nNA
    assert(size(x, 1) == size(y, 1))
    nr = size(x, 1)
    ncx = size(x, 2)
    ncy = size(y, 2)
    nNA_ext = 0

    # This matrix will hold preprocessed entries that can be simply multiplied together to get the
    # numerator
    multMatX = zeros(nr, ncx)
    multMatY = zeros(nr, ncy)
    # Number of NA entries in each column
    nNAentriesX = zeros(ncx)
    nNAentriesY = zeros(ncy)
    # Flag indicating whether the mean of each column is NA
    NAmedX = falses(ncx)
    NAmedY = falses(ncy)

    # Decide how many threads to use
    int nt = useNThreads( ncx* ncy, *nThreads);

    double * aux[MxThreads];
    for (int t=0; t < nt; t++)
    {
       if ( (aux[t] = (double *) malloc(6*nr * sizeof(double)))==NULL)
       {
         *err = 1;
         Rprintf("cor1: memory allocation error. The needed block is very small... suspicious.\n");
         for (int tt = t-1; tt>=0; tt--) free(aux[tt]);
         free(NAmedY); free(NAmedX); free(nNAentriesY); free(nNAentriesX); free(multMatY); free(multMatX);
         return;
       }
    }

    cor1ThreadData thrdInfoX[MxThreads];
    cor1ThreadData thrdInfoY[MxThreads];
    cor2ThreadData thrdInfo[MxThreads];

    # Prepare the multMat columns in X and Y

    # Rprintf(" ..preparing columns in x\n");
    cptd = Array{colPrepThreadData}(MxThreads)
    thr = Array{pthread_t}(MxThreads)
    status = Array{Int64}(MxThreads)  # normal: 0

    progressCounter pcX, pcY;
    int pearsonX = 0, pearsonY = 0;

    # Prepare columns in X

    pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
    pcX.i = 0; pcX.n = ncx;

    if robustX
        # calculate threadPrepColBicor using data cptd in each thr
    else
        # calculate threadPrepColCor using data cptd in each thr

    # If the fallback method is to re-do everything in Pearson, check whether any columns had zero MAD.
    if (*fallback==3)
    {
      for (int t=0; t<nt; t++) if (thrdInfoX[t].zeroMAD > 0)
      {
        pearsonX = 1;
        if (*verbose)
          Rprintf("Warning in bicor(x, y): thread %d of %d reported zero MAD in column %d of x. %s",
                  t, nt, thrdInfoX[t].zeroMAD, "Switching to Pearson calculation for x.\n");
      }
      if (pearsonX==1) # Re-do all column preparations
      {
        for (int t = 0; t < nt; t++) thrdInfoX[t].fallback = 4;

        pthread_mutex_t mutex2 = PTHREAD_MUTEX_INITIALIZER;
        pcX.i = 0; pcX.n = ncx;

        # calculate threadPrepColCor using data cptd in each thr
      }
    }

    # Prepare columns in Y

    # Rprintf(" ..preparing columns in y\n");
    pthread_mutex_t mutex1Y = PTHREAD_MUTEX_INITIALIZER;
    pcY.i = 0; pcY.n = ncy;

    if robustY
        # calculate threadPrepColBicor using data cptd in each thr
    else
        # calculate threadPrepColCor using data cptd in each thr

    # If the fallback method is to re-do everything in Pearson, check whether any columns had zero MAD.
    if (*fallback==3)
    {
      for (int t=0; t<nt; t++) if (thrdInfoY[t].zeroMAD > 0)
      {
        pearsonY = 1;
        if (*verbose)
          Rprintf("Warning in bicor(x, y): thread %d of %d reported zero MAD in column %d of y. %s",
                  t, nt, thrdInfoY[t].zeroMAD, "Switching to Pearson calculation for y.\n");
      }
      if (pearsonY==1) # Re-do all column preparations
      {
        for (int t = 0; t < nt; t++) thrdInfoY[t].fallback = 4;

        pthread_mutex_t mutex2Y = PTHREAD_MUTEX_INITIALIZER;
        pcY.i = 0;
        pcY.n = ncy;

        # calculate threadPrepColCor using data cptd in each thr
      }
    }

    result = multMatX' * multMatY

    # Remedial calculations

    # NA out all rows and columns that need it and check for values outside of [-1, 1]

    NA2ThreadData  natd[MxThreads];
    pcX.i = 0; pcY.i = 0;  # reset the progress counter
    pthread_t  thr2[MxThreads];

    for (int t=0; t<nt; t++)
    {
      natd[t].x = &thrdInfo[t];
      natd[t].pci = &pcX;
      natd[t].pcj = &pcY;
      status[t] = pthread_create_c(&thr2[t], NULL, threadNAing, (void *) &natd[t], thrdInfoX[t].threaded);
      if (status[t]!=0)
      {
         Rprintf("Error in bicor(x,y): thread %d could not be started successfully. Error code: %d.\n%s",
                 t, status[t], "*** WARNING: RETURNED RESULTS WILL BE INCORRECT. ***");
            *err = 2;
      }
    }

    for (int t=0; t<nt; t++)
       if (status[t]==0) pthread_join_c(thr2[t], NULL, thrdInfoX[t].threaded);

    nNA = nNA_ext

    # Clean up
    for (int t=nt-1; t >= 0; t--) free(aux[t]);
    free(NAmedY);
    free(NAmedX);
    free(nNAentriesY);
    free(nNAentriesX);
    free(multMatY);
    free(multMatX);
end

function cor(x::AbstractArray, y::AbstractArray, cosineX::Bool, cosineY::Bool,
             verbose::Bool)::Void
    """
    This can actually be relatively slow, since the search for calculations that need to be done is not
    parallel, so one thread may have to traverse the whole matrix. I can imagine parallelizing even that
    part, but for now leave it as is as this will at best be a minuscule improvement.
    """
    double * result, int *nNA
    assert(size(x, 1) == size(y, 1))
    nr = size(x, 1)
    ncx = size(x, 2)
    ncy = size(y, 2)
    size_t nNA_ext = 0

    # This matrix will hold preprocessed entries that can be simply multiplied together to get the
    # numerator
    multMatX = zeros(nr, ncx)
    multMatY = zeros(nr, ncy)
    # Number of NA entries in each column
    nNAentriesX = zeros(ncx)
    nNAentriesY = zeros(ncy)
    # Flag indicating whether the mean of each column is NA
    NAmeanX = falses(ncx)
    NAmeanY = falses(ncy)

    # Decide how many threads to use
    int nt = useNThreads( ncx* ncy, *nThreads);

    cor1ThreadData thrdInfoX[MxThreads];
    cor1ThreadData thrdInfoY[MxThreads];
    cor2ThreadData thrdInfo[MxThreads];

    # Prepare the multMat columns in X and Y
    cptd = Array{colPrepThreadData}(MxThreads)
    thr = Array{pthread_t}(MxThreads)
    status = Array{Int64}(MxThreads)  # normal: 0

    progressCounter pcX, pcY;

    # Prepare columns in X
    pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
    pcX.i = 0; pcX.n = ncx;
    # calculate threadPrepColCor using data cptd in each thr with thrdInfoX

    # Prepare columns in Y
    pthread_mutex_t mutex1Y = PTHREAD_MUTEX_INITIALIZER;
    pcY.i = 0; pcY.n = ncy;
    # calculate threadPrepColCor using data cptd in each thr with thrdInfoY

    result = multMatX' * multMatY

    # Remedial calculations

    # NA out all rows and columns that need it and check for values outside of [-1, 1]
    NA2ThreadData  natd[MxThreads];
    pcX.i = 0; pcY.i = 0;  # reset the progress counters
    pthread_t  thr2[MxThreads];

    for (int t=0; t<nt; t++)
    {
      natd[t].x = &thrdInfo[t];
      natd[t].pci = &pcX;
      natd[t].pcj = &pcY;
      status[t] = pthread_create_c(&thr2[t], NULL, threadNAing, (void *) &natd[t], thrdInfoX[t].threaded);
      if (status[t]!=0)
      {
         Rprintf("Error in cor(x,y): thread %d could not be started successfully. Error code: %d.\n%s",
                 t, status[t], "*** WARNING: RETURNED RESULTS WILL BE INCORRECT. ***");
            *err = 2;
      }
    }

    for (int t=0; t<nt; t++)
      if (status[t]==0)  pthread_join_c(thr2[t], NULL, thrdInfoX[t].threaded);

    nNA = nNA_ext;

    # clean up and return
    free(NAmeanY);
    free(NAmeanX);
    free(nNAentriesY);
    free(nNAentriesX);
    free(multMatY);
    free(multMatX);
end


text:
#md # !!! note
#
#     This tutorial is under construction, but the code below is already functional.
#
# Driver that computes the lid-driven cavity benchmark at low Reynolds numbers
# when using a mixed FE Q(k)/Pdisc(k-1).

# Load Gridap library
using Gridap

# Discrete model
n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)

# Define Dirichlet boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])

# Define reference FE (Q2/P1(disc) pair)
order = 2
reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffeₚ = ReferenceFE(lagrangian,Float64,order-1;space=:P)

# Define test FESpaces
V = TestFESpace(model,reffeᵤ,labels=labels,dirichlet_tags=["diri0","diri1"],conformity=:H1)
Q = TestFESpace(model,reffeₚ,conformity=:L2,constraint=:zeromean)
Y = MultiFieldFESpace([V,Q])

# Define trial FESpaces from Dirichlet values
u0 = VectorValue(0,0)
u1 = VectorValue(1,0)
U = TrialFESpace(V,[u0,u1])
P = TrialFESpace(Q)
X = MultiFieldFESpace([U,P])

# Define triangulation and integration measure
degree = order
Ωₕ = Triangulation(model)
dΩ = Measure(Ωₕ,degree)

# Define bilinear and linear form
f = VectorValue(0.0,0.0)
a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )dΩ
l((v,q)) = ∫( v⋅f )dΩ

# Build affine FE operator
op = AffineFEOperator(a,l,X,Y)

# Solve
uh, ph = solve(op)

# Export results to vtk
writevtk(Ωₕ,"results",order=2,cellfields=["uh"=>uh,"ph"=>ph])


text:
# Change the double **matrix declarations and access in the C code to use a flat
# column-major float array.

# In function emd:
#  "double **cost" into "double *cost"
#  "cost[var->row][var->col]" => "cost[(var->row) + n_x * (var->col)]"
#  "cost[i][j]" => "cost[i + n_x * j]"
#  "cost[basis[i]->row][basis[i]->col]" => "cost[(basis[i]->row) + n_x * (basis[i]->col)]"
# and similar for double **flows...

# In function initialize_flow:
#  "double **cost" into "double *cost"

function translate_c_file(srcfile, destfile)
    str = open(srcfile, "r") do fh
        read(fh, String)
    end
    str = replace(str, "#include <emd.h>",                   "#include <emd_translated.h>")

    str = replace(str, "double **cost",                      "double *cost")

    str = replace(str, "cost[var->row][var->col]",           "cost[(var->row) + n_x * (var->col)]")
    str = replace(str, "cost[i][j]",                         "cost[i + n_x * j]")
    str = replace(str, "cost[basis[i]->row][basis[i]->col]", "cost[(basis[i]->row) + n_x * (basis[i]->col)]")

    str = replace(str, "double **flows",                      "double *flows")
    str = replace(str, "memset(flows[i], 0, n_y*sizeof(double));", "for(j = 0; j < n_y; j++) {flows[i + n_x * j]=0.0;};")
    str = replace(str, "flows[basis[i]->row][basis[i]->col]", "flows[(basis[i]->row) + n_x * (basis[i]->col)]")

    open(destfile, "w") do fh
        print(fh, str)
    end

    printstyled("Translated double matrix references in file $(srcfile) and saved in $(destfile)\n", color=:green)
end

translate_c_file(ARGS[1], ARGS[2])


text:
using Test
using WaveProp.Geometry
using WaveProp.Interpolation
using WaveProp.Integration
using LinearAlgebra

@testset "Reference segment" begin
    qstd = GaussLegendre(20)
    d     = domain(qstd)
    # simple test on smooth integrand
    for shand in [IMT(), Kress(), KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        @test isapprox(integrate(cos,q),sin(1),rtol=1e-2)
    end
    # non-smooth integrand, singularity at 0
    f     = (x) -> log(abs(x))
    I     = integrate(f,d;rtol=1e-16)
    for shand in [IMT(), Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        Ia    = integrate(f,q)
        @test isapprox(Ia,I,rtol=1e-3)
        # check that the `naive` integration woudl have failed the test
        Istd    = integrate(f,qstd)
        @test !isapprox(Istd,I,rtol=1e-3)
    end
    # non-smooth integrand, singularity somewhere inside
    s     = rand()
    f     = (x) -> x[1]==s ? 0.0 : log(abs(x[1]-s))
    I     = integrate(f,d;rtol=1e-8)
    for shand in [IMT(), Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        x,w   = singular_quadrature(q,s)
        Ia    = integrate(f,x,w)
        @debug I - Ia
        @test isapprox(Ia,I,rtol=1e-3)
        # check that the `naive` integration woudl have failed the test
        Istd    = integrate(f,q)
        @test !isapprox(Istd,I,rtol=1e-3)
    end
    # quadrature for integrating a factored kernel singular at location s
    k = (x) -> x[1]==s ? 0. : log(abs(x[1]-s))
    ϕ = (x) -> cos(x[1])
    f = (x) -> k(x)*ϕ(x)
    I = integrate(f,d)
    for shand in [IMT(), Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        x,w   = singular_quadrature(k,q,s)
        Ia    = integrate(ϕ,x,w)
        @debug I - Ia
        @test isapprox(Ia,I,rtol=1e-5)
        # check that the `naive` integration woudl have failed the test
        Istd    = integrate(f,q)
        @test !isapprox(Istd,I,rtol=1e-5)
    end
    # weights for integrating a factored kernel singular at location s
    x  = qnodes(qstd)
    for shand in [IMT(), Kress(),KressP()]
        q     = SingularQuadratureRule(qstd,shand)
        ws    = singular_weights(k,qstd,q,s)
        Ia    = integrate(ϕ,x,ws)
        @test isapprox(Ia,I,rtol=1e-3)
    end
end

@testset "Duffy" begin
    q1d   = GaussLegendre(5)
    qstd  = TensorProductQuadrature(q1d,q1d)
    duffy = Duffy()
    qsin  = SingularQuadratureRule(qstd,duffy)
    x,w  = qsin()
    # regular kernel
    k    = x -> cos(x[1])*sin(x[2])
    Ie    = integrate(k,ReferenceTriangle)
    Ia   = integrate(k,qsin)
    @test isapprox(Ie,Ia,rtol=1e-5)
    # singular kernel at right vertex
    s    = SVector(1,0)
    k    = x -> 1/norm(x-s)
    Ia   = integrate(k,qsin)
    Ie    = integrate(k,ReferenceTriangle) # uses quadgk
    @test isapprox(Ie,Ia,rtol=1e-5)
    # singular kernel at wrong vertex (should not be accurate, see test below)
    s    = SVector(0,1)
    k    = x -> 1/norm(x-s)
    Ie    = integrate(k,ReferenceTriangle)
    Ia   = integrate(k,qsin)
    @test !isapprox(Ie,Ia,rtol=1e-5)
end

@testset "2d Kress" begin
    q1d   = GaussLegendre{10}()
    qstd  = TensorProductQuadrature(q1d,q1d)
    s1d   = Kress(order=2)
    sing_handler = TensorProductSingularityHandler(s1d,s1d)
    qsin = SingularQuadratureRule(qstd,sing_handler)
    x,w  = qsin()
    # regular kernel
    k(x) = cos(x[1])
    Ie    = integrate(k,ReferenceSquare)
    Isin = sum(k.(x).*w)
    @test isapprox(Ie,Isin,rtol=1e-2)
    # singular kernel at origin
    k(x) = 1/norm(x)
    Ie    = integrate(k,ReferenceSquare)
    Isin = sum(k.(x).*w)
    @test isapprox(Ie,Isin,rtol=1e-2)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
@testset "Simple Concretization Gala" begin
  # Did you pregame the gala?
  @test @concretization(Pregame) == Set{Type}([String])
end

module SimpleConcretization

using Test
using Hyperspecialize

@testset "Simple Concretization" begin

  # Can we concretize an abstract type to a subset of its concrete types?
  @test (@concretize Integer [Int64]) == Set{Type}([Int64])
  @test @concretization(Integer) == Set{Type}([Int64])
  # Can we widen our type?
  @test @widen(Integer, Int8) == Set{Type}([Int8, Int64])
  @test @concretization(Integer) == Set{Type}([Int8, Int64])
  # If we widen again with a duplicate type, it's not repeated right?
  @test (@widen Integer (Int8, Int16)) == Set{Type}([Int8, Int16, Int64])
  @test @concretization(Integer) == Set{Type}([Int8, Int16, Int64])

  # What is the concretization without any previous calls to concretize?
  @test @concretization(Unsigned) == nothing
  # We should be able to widen the concretization of an abstract type to include a non-subtype.
  @test @widen(Unsigned, Bool) == Set{Type}([UInt128, UInt16, UInt32, UInt64, UInt8, Bool])
  @test @concretization(Unsigned) == Set{Type}([UInt128, UInt16, UInt32, UInt64, UInt8, Bool])

  # Can we concretize a non-type without any previous calls to concretize?
  @test_throws ErrorException @concretize(NotAType)
  @test @concretize(NotAType, []) == Set{Type}()
  # We should be able to widen the concretization of a non-type
  @test @widen(NotAType, Bool) == Set{Type}([Bool])
  @test @concretization(NotAType) == Set{Type}([Bool])

  # What happens if we widen a non-concretized type?
  @test @widen(Signed, Bool) == Set{Type}([BigInt, Int128, Int16, Int32, Int64, Int8, Bool])
  @test @concretization(Signed) == Set{Type}([BigInt, Int128, Int16, Int32, Int64, Int8, Bool])

  # Can we concretize a concrete type to be something else?
  @test @concretize(Int8, Bool) == Set{Type}([Bool])
  @test @concretization(Int8) == Set{Type}([Bool])

  # Can we widen a non-existent non-concretized type?
  @test_throws ErrorException @widen(StillNotAType, Float32)
  @test @concretize(StillNotAType, []) == Set{Type}([])
  @test @widen(StillNotAType, Float32) == Set{Type}([Float32])
  @test @concretization(StillNotAType) == Set{Type}([Float32])

  # Can we concretize a type to be nothing?
  @test @concretize(Int16, []) == Set{Type}()
  @test @concretization(Int16) == Set{Type}()

  # We can concretize using a variable.
  x = Int64
  @test @concretize(Int32, [x]) == Set{Type}([Int64])
  @test @concretization(Int32) == Set{Type}([Int64])

  # We cannot concretize the contents of a variable.
  x = Int64
  @test @concretize(x, (Int32,)) == Set{Type}([Int32])
  @test @concretization(x) == Set{Type}([Int32])
  @test @concretization(Int64) == nothing

  # What if we set the concretization before we define the type?
  @test @concretize(Wobble, []) == Set{Type}([])
  @test @concretization(Wobble) == Set{Type}([])

  # Cannot reconcretize
  @test_throws ErrorException (@concretize Signed Int64)

end

struct Wobble
  w::Float32
end

@testset "Simple Concretization Afterparty" begin

  # Concretization is still messed up.
  @test @concretization(Wobble) == Set{Type}([])

  # We can add the type after...
  @test @widen(Wobble, Wobble) == Set{Type}([Wobble])

end

end #module


text:
# Paper_Experiment_Functions.jl
#= ---------------- #
Set of datatype and functions for efficient testing.
# ---------------- =#

module TestFunctions

using BayesianSVM
using Distributions
using KernelFunctions

export TestingModel
export DefaultParameters, BSVMParameters
export CreateModel, TrainModel, RunTests, ProcessResults, PrintResults, WriteResults
export ComputePrediction, ComputePredictionAccuracy

#Datatype for containing the model, its results and its parameters
type TestingModel
  MethodName::String #Name of the method
  DatasetName::String #Name of the dataset
  ExperimentType::String #Type of experiment
  MethodType::String #Type of method used ("SVM","BSVM","ECM","GPC")
  Param::Dict{String,Any} #Some paramters to run the method
  Results::Dict{String,Any} #Saved results
  Model::Any
  TestingModel(methname,dataset,exp,methtype) = new(methname,dataset,exp,methtype,Dict{String,Any}(),Dict{String,Any}())
  TestingModel(methname,dataset,exp,methtype,params) = new(methname,dataset,exp,methtype,params,Dict{String,Any}())
end

#Create a default dictionary
function DefaultParameters()
  param = Dict{String,Any}()
  param["ϵ"]= 1e-8 #Convergence criteria
  param["BatchSize"] = 10 #Number of points used for stochasticity
  param["Kernel"] = "rbf" # Kernel function
  param["Θ"] = 1.0 # Hyperparameter for the kernel function
  param["γ"] = 1.0 #Variance of introduced noise
  param["M"] = 32 #Number of inducing points
  param["Window"] = 5 #Number of points used to check convergence (smoothing effect)
  param["Verbose"] = 0 #Verbose
  return param
end

#Create a default parameters dictionary for BSVM
function BSVMParameters(;Stochastic=true,NonLinear=true,Sparse=true,ALR=true,Autotuning=false,main_param=DefaultParameters())
  param = Dict{String,Any}()
  param["Stochastic"] = Stochastic #Is the method stochastic
  param["Sparse"] = Sparse #Is the method using inducing points
  param["NonLinear"] = NonLinear #Is the method using kernels
  param["Intercept"] = NonLinear ? false : true;
  param["ALR"] = ALR #Is the method using adpative learning rate (in case of the stochastic case)
  param["AutoTuning"] = Autotuning #Is hyperoptimization performed
  param["ATFrequency"] = 1 #How even autotuning is performed
  param["κ"] = 1.0;  param["τ"] = 100; #Parameters for learning rate of autotuning
  param["κ_s"] = 1.0;  param["τ_s"] = 100; #Parameters for learning rate of Stochastic gradient descent when ALR is not used
  param["ϵ"] = main_param["ϵ"]; param["Window"] = main_param["Window"]; #Convergence criteria (checking parameters norm variation on a window)
  param["Kernels"] = [Kernel(main_param["Kernel"],1.0,params=main_param["Θ"])] #Kernel creation (standardized for now)
  param["Verbose"] = main_param["Verbose"] ? 3 : 0 #Verbose
  param["BatchSize"] = main_param["BatchSize"] #Number of points used for stochasticity
  param["M"] = main_param["M"] #Number of inducing points
  param["γ"] = main_param["γ"] #Variance of introduced noise
  return param
end

#Create a model given the parameters passed in p
function CreateModel(tm::TestingModel,X_training,y_training) #tm testing_model, p parameters
  tm.Model = BSVM(X_training,y_training;Stochastic=tm.Param["Stochastic"],batchSize=tm.Param["BatchSize"],Sparse=tm.Param["Sparse"],m=tm.Param["M"],NonLinear=tm.Param["NonLinear"],
  kernels=tm.Param["Kernels"],Autotuning=tm.Param["AutoTuning"],autotuningfrequency=tm.Param["ATFrequency"],AdaptativeLearningRate=tm.Param["ALR"],Intercept=tm.Param["Intercept"],κ_s=tm.Param["κ_s"],τ_s = tm.Param["τ_s"],ϵ=tm.Param["ϵ"],γ=tm.Param["γ"],
  κ_Θ=tm.Param["κ"],τ_Θ=tm.Param["τ"],smoothingWindow=tm.Param["Window"],VerboseLevel=tm.Param["Verbose"])
end

#Train the model on trainin set (X,y) for #iterations
function TrainModel(tm::TestingModel,iterations)
  time_training = 0;
  tm.Model.nEpochs = iterations
  time_training = @elapsed tm.Model.Train()
  return time_training;
end

#Run tests accordingly to the arguments and save them
function RunTests(tm::TestingModel,X_test,y_test;accuracy::Bool=false,brierscore::Bool=false,logscore::Bool=false)
  if accuracy
    push!(tm.Results["accuracy"],TestAccuracy(y_test,ComputePrediction(tm,X_test)))
  end
  y_predic_acc = 0
  if brierscore
    y_predic_acc = ComputePredictionAccuracy(tm::TestingModel, X_test)
    push!(tm.Results["brierscore"],TestBrierScore(y_test,y_predic_acc))
  end
  if logscore
    if y_predic_acc == 0
      y_predic_acc = ComputePredictionAccuracy(tm::TestingModel,X_test)
    end
    push!(tm.Results["logscore"],TestLogScore(y_test,y_predic_acc))
  end
end

#Compute the mean and the standard deviation and assemble in one result
function ProcessResults(tm::TestingModel,writing_order)
  all_results = Array{Float64,1}()
  names = Array{String,1}()
  for name in writing_order
    result = [mean(tm.Results[name]), std(tm.Results[name])]
    all_results = vcat(all_results,result)
    names = vcat(names,name)
  end
  if haskey(tm.Results,"allresults")
    tm.Results["allresults"] = vcat(tm.Results["allresults"],all_results')
  else
    tm.Results["allresults"] = all_results'
  end
  if !haskey(tm.Results,"names")
    tm.Results["names"] = names
  end
end

function PrintResults(results,method_name,writing_order)
  i = 1
  for category in writing_order
    println("$category : $(results[i*2-1]) ± $(results[i*2])")
    i+=1
  end
end

function WriteResults(tm::TestingModel,location)
  fold = String(location*"/"*tm.ExperimentType*"Experiment_"*tm.DatasetName*"Dataset")
  if !isdir(fold); mkdir(fold); end;
  writedlm(String(fold*"/Results_"*tm.MethodName*".txt"),tm.Results["allresults"])
end

#Return predicted labels (-1,1) for test set X_test
function ComputePrediction(tm::TestingModel, X_test)
  return sign.(tm.Model.Predict(X_test))
end

#Return prediction certainty for class 1 on test set X_test
function ComputePredictionAccuracy(tm::TestingModel, X_test)
  return tm.Model.PredictProba(X_test)
end

#Return Accuracy on test set
function TestAccuracy(y_test, y_predic)
  return 1-sum(1-y_test.*y_predic)/(2*length(y_test))
end
#Return Brier Score
function TestBrierScore(y_test, y_predic)
  return sum(((y_test+1)./2 - y_predic).^2)/length(y_test)
end
#Return Log Score
function TestLogScore(y_test, y_predic)
  return sum((y_test+1)./2.*log(y_predic)+(1-(y_test+1)./2).*log(1-y_predic))/length(y_test)
end
#Return ROC
function TestROC(y_test,y_predic)
    nt = length(y_test)
    truepositive = zeros(npoints); falsepositive = zeros(npoints)
    truenegative = zeros(npoints); falsenegative = zeros(npoints)
    thresh = collect(linspace(0,1,npoints))
    for i in 1:npoints
      for j in 1:nt
        truepositive[i] += (yp[j]>=thresh[i] && y_test[j]>=0.9) ? 1 : 0;
        truenegative[i] += (yp[j]<=thresh[i] && y_test[j]<=-0.9) ? 1 : 0;
        falsepositive[i] += (yp[j]>=thresh[i] && y_test[j]<=-0.9) ? 1 : 0;
        falsenegative[i] += (yp[j]<=thresh[i] && y_test[j]>=0.9) ? 1 : 0;
      end
    end
    return (truepositive./(truepositive+falsenegative),falsepositive./(truenegative+falsepositive))
end

end #end of module


text:
include("core.jl")
include("scaling.jl")
include("algorithms.jl")
include("edge.jl")
include("io.jl")
include("readnrrd.jl")
@linux_only include("readremote.jl")
@windows_only include("readremote.jl")
@osx_only include("readOSX.jl")


text:
module MochaRegression

using Compat
using Mocha
using MochaRegressionWrapper
using MLUtilities
using Utilities
using DataModel

export MochaRegressionModel

"""
Neural Network Regression using Mocha

fields:
    - backend: Mocha Backend (CPU or GPU)
    - mochaNet: Mocha wrapper
    - labels
    - ntrain: number of data items
    - batchsize
"""
type MochaRegressionModel <: AbstractDataModel
  backend::Backend
  mochaNet::MochaRegressionWrapper.MWrap
  ground_truth::Array{Float64,1}
  ntrain::Float64 # effective training set size
  batchsize::Int64
  temperature::Float64
end

function MochaRegressionModel(
  datax,datac,
  modelfactory::Function,
  backend::Backend;
  ntrain::Int = length(datac),
  batchsize::Int = 100,
  temperature::Float64 = 1.0,
  do_shuffle::Bool = true,
  do_storage::Bool = false
  )

  data_layer = MemoryDataLayer(name = "data",
                               data = Array[datax,datac],
                               tops = [:data,:ground_truth],
                               batch_size = batchsize,
                               shuffle = do_shuffle)

  mochaNet = MochaRegressionWrapper.MWrap(data_layer,
                                 modelfactory,
                                 "MochaRegressionNet",
                                 backend,
                                 do_storage
                                 )
  MochaRegressionModel(backend,mochaNet,datac[:],ntrain,batchsize,temperature)
end

"""
`DataModel.getllik(dms::MochaRegressionModel)` returns a function
that returns a subsampled likelihood estimate.
"""
function DataModel.getllik(dms::MochaRegressionModel)
    function llik(x)
        #@assert dms.ntrain == dms.batch_size
        loss = MochaRegressionWrapper.evaluateTestNN(dms.mochaNet, x, dms.batchsize)
        return loss
    end
end

"""
`DataModel.getgrad(dms::MochaRegressionModel)` returns function returning the stochastic gradient function.
"""
function DataModel.getgrad(dms::MochaRegressionModel)
    function grad(x)
        (llik, grad) = MochaRegressionWrapper.forwardbackward(dms.mochaNet, x, regu_coef = 0.)
        return grad * dms.ntrain
    end
end

"""
`getparams(dms::MochaRegressionModel)` reads out the current weights.
"""
function getparams(dms::MochaRegressionModel)
    MochaRegressionWrapper.getparams(dms.mochaNet)
end
"""
`setparams!(dms::MochaRegressionModel,para::Array{Float64,1` sets the weights to `para`.
"""
function setparams!(dms::MochaRegressionModel,para::Array{Float64,1})
    MochaRegressionWrapper.setparams!(dms.mochaNet,para)
end

"""
`fetchnparams(dms::MochaRegressionModel)` returns the number of parameters of the model.
"""
function fetchnparams(dms::MochaRegressionModel)
  MochaRegressionWrapper.getnparams(dms.mochaNet)
end

#-------------------
# Initialisation functions
#-------------------
"""
Xavier Initialisation
"""
function init_xavier(dms::MochaRegressionModel)
  MochaRegressionWrapper.init_xavier(dms.mochaNet)
end
"""
Initialisation based on fan-in.
"""
function init_simple_fanin(dms::MochaRegressionModel)
  MochaRegressionWrapper.init_simple_fanin(dms.mochaNet)
end
"""
Gaussian Initialisation
"""
function init_gaussian(dms::MochaRegressionModel, initvar::Float64)
  MochaRegressionWrapper.init_gaussian(dms.mochaNet, initvar)
end
"""
Uniform Initialisation
"""
function init_uniform(dms::MochaRegressionModel, initvar::Float64)
  MochaRegressionWrapper.init_uniform(dms.mochaNet, initvar)
end

#-------------------
# evaluation functions
#-------------------
"""
Evaluate loss based on predictive probabilities
"""
function evaluate(dms::MochaRegressionModel,
  pred::Array{Float64,2}
  )
  err = pred[:]-dms.ground_truth
  loss = -0.5*dot(err,err)
  @dict(loss)
end
"""
Evaluate loss at weights `x`.
"""
function evaluate(dms::MochaRegressionModel,
  x::Vector{Float64}
  )
  loss = MochaRegressionWrapper.evaluateTestNN(dms.mochaNet, x, dms.batchsize)
  @dict(loss)
end

"""
Get predictions given weights `x`
"""
function getpredictions(dms::MochaRegressionModel,
  x::Vector{Float64}
  )
  predictions = MochaRegressionWrapper.getpredictions(dms.mochaNet, x, dms.batchsize)
  @dict(predictions)
end

end


text:
"""
    only(xs)

Like `first(xs)`, but asserts that the iterator `xs` only had one element.
"""
function only(xs)
    val, state = iterate(xs)
    @assert iterate(xs, state) == nothing
    return val
end

function better_readline(stream = stdin)
    #TODO doing this better still requires depending on the REPL.Terminals module
    if !isopen(stream)
        Base.reseteof(stream)
        @assert(isopen(stream))
    end
    if Sys.iswindows()
        # Apparently in windows it is already pretty nice
        return readline(stream)
    else
        #TODO: Put in raw mode, drop control characters etc
        return readline(stream)
    end
end


text:
mutable struct csi_t
    ref_count::Int64
    status::Int32
    finished::UInt32
    hooks::Ptr{Nothing}
end

mutable struct Interpreter <: GraphicsDevice
    ptr::Ptr{csi_t}

    function Interpreter(ptr::Ptr{Nothing})
        self = new(ptr)
        finalizer(destroy, self)
        self
    end
end

function destroy(csi::Interpreter)
    if csi.ptr == C_NULL
        return
    end
    ccall((:cairo_script_interpreter_destroy,_jl_libcsi), Nothing, (Ptr{Nothing},), csi.ptr)
    csi.ptr = C_NULL
    nothing
end

function Interpreter()
    ptr = ccall((:cairo_script_interpreter_create,_jl_libcsi),
                Ptr{Nothing},())
    Interpreter(ptr)
end

mutable struct cl <: GraphicsDevice
    s::CairoSurface
end

mutable struct InterpreterHooks <: GraphicsDevice
	closure::Ptr{Nothing}
	surface_create::Ptr{Nothing}
	surface_destroy::Ptr{Nothing}
	context_create::Ptr{Nothing}
	context_destroy::Ptr{Nothing}
	show_page::Ptr{Nothing}
	copy_page::Ptr{Nothing}
	create_source_image::Ptr{Nothing}

	function InterpreterHooks(;
                closure = C_NULL,
                surface_create = C_NULL,
                surface_destroy = C_NULL,
                context_create = C_NULL,
                context_destroy = C_NULL,
                show_page = C_NULL,
                copy_page = C_NULL,
                create_source_image = C_NULL)

		#self = new(cl(CairoSurface(Ptr{Nothing}(0))),C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,C_NULL)
        self = new(closure,surface_create,surface_destroy,context_create,context_destroy,show_page,copy_page,create_source_image)
		self
	end

    #function InterpreterHooks()
end

# typedef enum _csi_status {
#     CSI_STATUS_SUCCESS = CAIRO_STATUS_SUCCESS,
#     CSI_STATUS_NO_MEMORY = CAIRO_STATUS_NO_MEMORY,
#     CSI_STATUS_INVALID_RESTORE = CAIRO_STATUS_INVALID_RESTORE,
#     CSI_STATUS_INVALID_POP_GROUP = CAIRO_STATUS_INVALID_POP_GROUP,
#     CSI_STATUS_NO_CURRENT_POINT = CAIRO_STATUS_NO_CURRENT_POINT,
#     CSI_STATUS_INVALID_MATRIX = CAIRO_STATUS_INVALID_MATRIX,
#     CSI_STATUS_INVALID_STATUS = CAIRO_STATUS_INVALID_STATUS,
#     CSI_STATUS_NULL_POINTER = CAIRO_STATUS_NULL_POINTER,
#     CSI_STATUS_INVALID_STRING = CAIRO_STATUS_INVALID_STRING,
#     CSI_STATUS_INVALID_PATH_DATA = CAIRO_STATUS_INVALID_PATH_DATA,
#     CSI_STATUS_READ_ERROR = CAIRO_STATUS_READ_ERROR,
#     CSI_STATUS_WRITE_ERROR = CAIRO_STATUS_WRITE_ERROR,
#     CSI_STATUS_SURFACE_FINISHED = CAIRO_STATUS_SURFACE_FINISHED,
#     CSI_STATUS_SURFACE_TYPE_MISMATCH = CAIRO_STATUS_SURFACE_TYPE_MISMATCH,
#     CSI_STATUS_PATTERN_TYPE_MISMATCH = CAIRO_STATUS_PATTERN_TYPE_MISMATCH,
#     CSI_STATUS_INVALID_CONTENT = CAIRO_STATUS_INVALID_CONTENT,
#     CSI_STATUS_INVALID_FORMAT = CAIRO_STATUS_INVALID_FORMAT,
#     CSI_STATUS_INVALID_VISUAL = CAIRO_STATUS_INVALID_VISUAL,
#     CSI_STATUS_FILE_NOT_FOUND = CAIRO_STATUS_FILE_NOT_FOUND,
#     CSI_STATUS_INVALID_DASH = CAIRO_STATUS_INVALID_DASH,
#     CSI_STATUS_INVALID_DSC_COMMENT = CAIRO_STATUS_INVALID_DSC_COMMENT,
#     CSI_STATUS_INVALID_INDEX = CAIRO_STATUS_INVALID_INDEX,
#     CSI_STATUS_CLIP_NOT_REPRESENTABLE = CAIRO_STATUS_CLIP_NOT_REPRESENTABLE,
#     CSI_STATUS_TEMP_FILE_ERROR = CAIRO_STATUS_TEMP_FILE_ERROR,
#     CSI_STATUS_INVALID_STRIDE = CAIRO_STATUS_INVALID_STRIDE,
#     CSI_STATUS_FONT_TYPE_MISMATCH = CAIRO_STATUS_FONT_TYPE_MISMATCH,
#     CSI_STATUS_USER_FONT_IMMUTABLE = CAIRO_STATUS_USER_FONT_IMMUTABLE,
#     CSI_STATUS_USER_FONT_ERROR = CAIRO_STATUS_USER_FONT_ERROR,
#     CSI_STATUS_NEGATIVE_COUNT = CAIRO_STATUS_NEGATIVE_COUNT,
#     CSI_STATUS_INVALID_CLUSTERS = CAIRO_STATUS_INVALID_CLUSTERS,
#     CSI_STATUS_INVALID_SLANT = CAIRO_STATUS_INVALID_SLANT,
#     CSI_STATUS_INVALID_WEIGHT = CAIRO_STATUS_INVALID_WEIGHT,

#     /* cairo-script-interpreter specific errors */

#     CSI_STATUS_INVALID_SCRIPT,
#     CSI_STATUS_SCRIPT_INVALID_TYPE,
#     CSI_STATUS_SCRIPT_INVALID_INDEX,
#     CSI_STATUS_SCRIPT_UNDEFINED_NAME,
#     CSI_STATUS_INTERPRETER_FINISHED,

#     _CSI_STATUS_SCRIPT_LAST_ERROR,
#     CSI_INT_STATUS_UNSUPPORTED
# } csi_status_t;


text:
using Random

struct AliasSampler{IndexType <: Integer, ProbabilityParamType <: Real}
    # Some alias_indices will be left uninitialized after the
    # constructor finishes. This is not a bug!
    alias_indices::Vector{IndexType}
    nonalias_probs::Vector{ProbabilityParamType}
end

function AliasSampler(IndexType::Type,
                weights::AbstractVector{ProbabilityParamType}
            ) where{ProbabilityParamType <: Real}
    @assert IndexType <: Integer

    n = length(weights)
    s = sum(weights)
    scf = n/s
    sc_probabilities = (weights*scf)::Vector{ProbabilityParamType}

    smalls = Vector{IndexType}(undef, n)
    smalls_idx = 0
    larges = Vector{IndexType}(undef, n)
    larges_idx = 0

    for idx in 1:n
        if sc_probabilities[idx] >= 1.0
            larges_idx += 1
            larges[larges_idx] = idx
        else
            smalls_idx += 1
            smalls[smalls_idx] = idx
        end
    end

    aliases = Vector{IndexType}(undef, n)

    while larges_idx > 0 && smalls_idx > 0
        sm_p_idx = smalls[smalls_idx]
        lg_p_idx = larges[larges_idx]
        aliases[sm_p_idx] = lg_p_idx

        # This is slightly better numerically than the more obvious: sc_probabilities[lg_p_idx] += sc_probabilities[sm_p_idx] - 1.0
        sc_probabilities[lg_p_idx] -= 1.0
        sc_probabilities[lg_p_idx] += sc_probabilities[sm_p_idx]

        if sc_probabilities[lg_p_idx] < 1.0
            smalls[smalls_idx] = lg_p_idx
            larges_idx -= 1
        else
            smalls_idx -= 1
        end
    end

    while larges_idx > 0
        sc_probabilities[larges[larges_idx]] = 2.0
        larges_idx -= 1
    end

    while smalls_idx > 0
        sc_probabilities[smalls[smalls_idx]] = 2.0
        smalls_idx -= 1
    end

    return AliasSampler(aliases, sc_probabilities)

end

function a_sample(
                        alias_sampler::AliasSampler{IndexType, ProbabilityParamType},
                        rng=Random.GLOBAL_RNG
                 )::IndexType where{IndexType <: Integer, ProbabilityParamType <: Real}

    idx = rand(rng, 1:length(alias_sampler.alias_indices))
    if alias_sampler.nonalias_probs[idx] >= 1.0
        return idx
    end
    if rand(rng) < alias_sampler.nonalias_probs[idx]
        return idx
    else
        return alias_sampler.alias_indices[idx]
    end
end


text:
using ArgParse
using EventStreamMatrix
using EventStreamGLM
using YAML
using BSplines

import Dates: now, Millisecond
import Random: seed!
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--outfolder"
            help="Folder to write out the simulated datasets"
            default="/tmp/"
        "--name"
            help="Simulated dataset file name prefix"
            default="sim"
        "--nsims"
            help="Number of simulated datasets to create"
            default=1
            arg_type=Int
        "--maxtime"
            help="Maximum amount of time to simulate for"
            default = 10.0
            arg_type=Float64
        "--intensity_ub"
            help="Upper bound for self-point intensity"
            default=10.0
            arg_type=Float64
        "--seed"
            help="RNG seed"
            default=534
            arg_type=Int
        "--config"
            help="Path to YAML configuration file"
            default="./config.yaml"
            arg_type=String
    end
    return parse_args(s)
end

function main()
    args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in args
        println("  $arg  =>  $val")
    end
    conf = YAML.load_file(abspath(args["config"]))
    for (k, v) in pairs(conf)
        println("$k -> $v")
    end
    breakpoints = conf["breakpoints"]
    order = conf["order"]
    labels = ["source_$i" for i in 1:conf["inputnum"]] :: Vector{String}
    basis = BSplineBasis(order, breakpoints)
    # generate input events, one of ("homogpois")
    eventstream = Tuple{Float64, Int}[]
    if conf["input_type"] == "homogpois"
        for (j, intensity) in enumerate(conf["input_intensity"])
            P = HomogPoissonProcess(intensity, args["maxtime"])
            points = rand(P)
            evns = [(t, j) for t in points] :: Vector{Tuple{Float64, Int}}
            eventstream = vcat(eventstream, evns)
        end
    end
    sort!(eventstream) # ascending temporal order
    ## Create kernels
    other_kernels = [Spline(basis, q) for q in conf["other_coefs"]]
    self_kernel = Spline(basis, conf["self_coefs"])

    proc = EventStreamProcess(eventstream, labels, args["maxtime"],basis, support(basis)[2], conf["lambda_0"], other_kernels, self_kernel)
    # Sample
    for j in 1:args["nsims"]
        seed!(args["seed"]+j)
        fname = joinpath(args["outfolder"], "$(args["name"])_$(j).yaml")
        println("Simulating $fname")
        simtime_start = now()
        outpoints = rand(proc)
        out_data = Dict{String, Any}("nsources" => [length(labels)], "output" => outpoints)
        for (j,l) in enumerate(labels)
            out_data[l] = [e[1] for e in eventstream if e[2] == j]
        end
        simtime_end = now()
        # Locate some of the configuation data in each dataset.
        out_data["maxtime"] = [args["maxtime"]]
        out_data["true_other_coefs"] = conf["other_coefs"]
        out_data["true_self_coefs"] = conf["self_coefs"]
        out_data["breakpoints"] = conf["breakpoints"]
        out_data["order"] = conf["order"]
        out_data["simtime"] = Millisecond(simtime_start - simtime_end).value
        YAML.write_file(fname, out_data)
    end
end

main()


text:
# === IMPORTED FUNCTIONS AND TYPES === #

import Distributions.rand!
import Distributions._rand!
import Distributions.logpdf
import Distributions.ContinuousUnivariateDistribution
import Distributions.ContinuousMultivariateDistribution

# === IMPORTED DISTRIBUTIONS FROM DISTRIBUTIONS.JL === #

import Distributions.Dirichlet
import Distributions.Poisson
import Distributions.Chisq
import Distributions.Normal
import Distributions.Categorical

# === HELPER FUNCTIONS TO SET HYPERPARAMS === #

"""Specifies RateGamma distribution by first two moments."""
function specify_gamma(_mean::Real, _var::Real)
    β = _mean / _var
    α = _mean * β
    RateGamma(α, β)
end

# === GAMMA DISTRIBUTION === #

"""
RateGamma distribution, parameterized by shape and rate.

p(x | α, β) = (β^α / Γ(α)) ⋅ x^(α + 1) ⋅ exp{-βx}

This is a light wrapper to Distributions.jl which
uses the alternative parameterization based on
shape and scale to characterize the gamma distribution.
"""
struct RateGamma
    α::Float64
    β::Float64
end

logpdf(g::RateGamma, x::Float64) = logpdf(Distributions.Gamma(g.α, 1 / g.β), x)
rand(g::RateGamma) = rand(Distributions.Gamma(g.α, 1 / g.β))
mean(g::RateGamma) = g.α / g.β
var(g::RateGamma) = g.α / (g.β * g.β)

posterior(count_var::Real, prior::RateGamma) = (
    RateGamma(prior.α + count_var, prior.β + 1)
)

# === SYMMETRIC DIRICHLET === #

struct SymmetricDirichlet <: ContinuousMultivariateDistribution
    conc::Float64
    dim::Int64
    log_normalizer::Float64
end

function SymmetricDirichlet(conc::Float64, dim::Int64)
    SymmetricDirichlet(
        conc,
        dim,
        (dim * lgamma(conc)) - lgamma(dim * conc)  # log normalizer.
    )
end

length(g::SymmetricDirichlet) = g.dim

function _rand!(
        r::Random.AbstractRNG,
        g::SymmetricDirichlet,
        x::AbstractVector{Float64}
    )

    rand!(Distributions.Gamma(g.conc, 1.0), x)
    x ./= sum(x)
end

posterior(counts::AbstractVector, prior::SymmetricDirichlet) = (
    Dirichlet(counts .+ prior.conc)
)

logpdf(g::SymmetricDirichlet, p::AbstractVector) = (
    mapreduce(log, +, p) * (g.conc - 1) - g.log_normalizer
)

# === SCALED INVERSE CHI SQUARED === #

struct ScaledInvChisq
    ν::Float64
    s2::Float64
end

function logpdf(g::ScaledInvChisq, θ::Float64)
    # See pg 576 of Gelman et al., Bayesian Data Analysis.
    hν = 0.5 * g.ν
    return (
        hν * log(hν)
        - lgamma(hν)
        + hν * log(g.s2)
        - (hν + 1) * log(θ)
        - hν * g.s2 / θ
    )
end

rand(g::ScaledInvChisq) = g.ν * g.s2 / rand(Chisq(g.ν))

# === NORMAL INVERSE CHI SQUARED === #

"""
A normal-inverse-chi-squared distribution is a convienent
prior on the mean and variance parameters for a univariate
Gaussian. It can be seen as a reparameterization of the
normal inverse gamma distribution.

If:

    σ² ~ Inv-χ²(ν, s)
    μ | σ² ~ Normal(m, σ² / k)

Then (μ, σ²) follows a normal-inverse-chi-squared distribution
with parameters (k, m, ν, s)

Parameters
----------
k : Number of psuedo-observations of m. Note that k > 0.
m : Prior on the mean parameter.
ν : Number of psuedo-observations of s. Note that ν > 0.
s2 : Prior on the variance parameter. Note that s2 > 0.

References
----------
 - Pgs 42-43, 67-69 of Gelman et al., Bayesian Data Analysis
 - Pg 134 of Murphy, Machine Learning: A Probabilistic Perspective.
"""
struct NormalInvChisq
    k::Float64
    m::Float64
    ν::Float64
    s2::Float64
end

function rand(g::NormalInvChisq)
    σ2 = rand(ScaledInvChisq(g.ν, g.s2))
    μ = g.m + randn() * sqrt(σ2 / g.k)
    return μ, σ2
end

function posterior(
        n::Int64,
        sum_x::Float64,
        sumsq_x::Float64,
        prior::NormalInvChisq,
    )

    (n == 0) && return prior

    k = prior.k + n
    ν = prior.ν + n
    m = (sum_x / k) + (prior.k * prior.m / k)
    s2 = (
        prior.ν * prior.s2
        + (sumsq_x - sum_x * sum_x / n)  # == sum_i (xᵢ - xbar)^2
        + (prior.k * n * (sum_x / n - prior.m)^2) / k
    ) / ν

    return NormalInvChisq(k, m, ν, s2)
end

function logpdf(g::NormalInvChisq, μ::Float64, σ2::Float64)
    return (
        logpdf(ScaledInvChisq(g.ν, g.s2), σ2)
        + logpdf(Normal(g.m, sqrt(σ2 / g.k)), μ)
    )
end


text:
struct CanonicalDict{T}
    dict::Dict{T, TeXExpr}
end

CanonicalDict{T}() where T = CanonicalDict(Dict{T, TeXExpr}())

Base.setindex!(d::CanonicalDict, val::TeXExpr, key) = (d.dict[key] = val)
Base.setindex!(d::CanonicalDict{Char}, val::TeXExpr, key::String) = (d[first(key)] = val)

function Base.getindex(d::CanonicalDict, key)
    val = d.dict[key]

    if val isa Char
        return val
    else
        return copy(val)
    end
end

Base.getindex(d::CanonicalDict{Char}, key::String) = d[first(key)]

Base.haskey(d::CanonicalDict, key) = haskey(d.dict, key)
Base.haskey(d::CanonicalDict{Char}, key::String) = haskey(d, first(key))

function Base.get(d::CanonicalDict, key, default)
    haskey(d, key) && return d[key]
    return default
end

# Each symbol or command has a unique canonical representation, either
# as a TeXExpr
const symbol_to_canonical = CanonicalDict{Char}()
const command_to_canonical = CanonicalDict{String}()

function canonical_expr(char::Char)
    haskey(symbol_to_canonical, char) && return symbol_to_canonical[char]
    return TeXExpr(:char, char)
end

canonical_expr(command::String) = get(command_to_canonical, command, nothing)

# Symbols missing from the REPL completion data
latex_symbols[raw"\neq"] = "≠"

function get_symbol_char(command)
    if !haskey(latex_symbols, command)
        @warn "unkown command $command"
        return '?'
    end

    return first(latex_symbols[command])
end

# Numbers
for char in join(0:9)
    symbol_to_canonical[char] = TeXExpr(:digit, char)
end

##
## Special commands
##

command_to_canonical[raw"\frac"] = TeXExpr(:argument_gatherer, [:frac, 2])
command_to_canonical[raw"\sqrt"] = TeXExpr(:argument_gatherer, [:sqrt, 1])

##
## Commands from the commands_data.jl file
##

for symbol in spaced_symbols
    symbol_expr = TeXExpr(:symbol, symbol)
    symbol_to_canonical[symbol] = TeXExpr(:spaced, symbol_expr)
end

# Special case for hyphen that must be replaced by a minus sign
symbol_to_canonical['-'] = TeXExpr(:spaced, TeXExpr(:symbol, '−'))

for command in spaced_commands
    symbol = get_symbol_char(command)
    symbol_expr = TeXExpr(:symbol, symbol)
    symbol_to_canonical[symbol] = command_to_canonical[command] = TeXExpr(:spaced, symbol_expr)
end

for command in underover_commands
    symbol = get_symbol_char(command)
    symbol_expr = TeXExpr(:symbol, symbol)
    symbol_to_canonical[symbol] = command_to_canonical[command] = TeXExpr(:underover, Any[symbol_expr, nothing, nothing])
end

for func in underover_functions
    command = "\\" * func
    command_to_canonical[command] = TeXExpr(:underover, Any[TeXExpr(:function, func), nothing, nothing])
end

for command in integral_commands
    symbol = get_symbol_char(command)
    symbol_expr = TeXExpr(:symbol, symbol)
    symbol_to_canonical[symbol] = command_to_canonical[command] = TeXExpr(:integral, Any[symbol_expr, nothing, nothing])
end

for func in generic_functions
    command = "\\" * func
    command_to_canonical[command] = TeXExpr(:function, func)
end

for (command, width) in space_commands
    command_to_canonical[command] = TeXExpr(:space, width)
end

for (symbol, width) in space_symbols
    symbol_to_canonical[symbol] = TeXExpr(:space, width)
end

for command in combining_accents
    combining_char = get_symbol_char(command)
    symbol_expr = TeXExpr(:symbol, combining_char)
    command_to_canonical[command] = TeXExpr(:argument_gatherer, [:combining_accent, 2, symbol_expr])
end

for symbol in punctuation_symbols
    symbol = first(symbol)
    symbol_to_canonical[symbol] = TeXExpr(:punctuation, symbol)
end

for symbol in delimiter_symbols
    symbol = first(symbol)
    symbol_to_canonical[symbol] = TeXExpr(:delimiter, symbol)
end

for name in font_names
    command = "\\math$name"
    command_to_canonical[command] = TeXExpr(:argument_gatherer, [:font, 2, Symbol(name)])
end

##
## Default behavior
##
# We put it at the end to avoid overwritting existing commands

for (command, symbol) in latex_symbols
    symbol = first(symbol)  # Convert String to Char

    if haskey(command_to_canonical, command) || haskey(symbol_to_canonical, symbol)
        continue
    end

    symbol_expr = TeXExpr(:symbol, [symbol])
    symbol_to_canonical[symbol] = symbol_expr
    command_to_canonical[command] = symbol_expr
end


text:
# ------------------------------------------------------------------------------
#
#  Gmsh Julia tutorial 4
#
#  Holes in surfaces, annotations, entity colors
#
# ------------------------------------------------------------------------------
import gmsh

gmsh.initialize()

gmsh.model.add("t4")

cm = 1e-02
e1 = 4.5 * cm; e2 = 6 * cm / 2; e3 =  5 * cm / 2
h1 = 5 * cm; h2 = 10 * cm; h3 = 5 * cm; h4 = 2 * cm; h5 = 4.5 * cm
R1 = 1 * cm; R2 = 1.5 * cm; r = 1 * cm
Lc1 = 0.01
Lc2 = 0.003

function hypot(a, b)
    return sqrt(a * a + b * b)
end

ccos = (-h5*R1 + e2 * hypot(h5, hypot(e2, R1))) / (h5*h5 + e2*e2)
ssin = sqrt(1 - ccos*ccos)

# We start by defining some points and some lines. To make the code shorter we
# can redefine a namespace:
factory = gmsh.model.geo
factory.addPoint(-e1-e2, 0    , 0, Lc1, 1)
factory.addPoint(-e1-e2, h1   , 0, Lc1, 2)
factory.addPoint(-e3-r , h1   , 0, Lc2, 3)
factory.addPoint(-e3-r , h1+r , 0, Lc2, 4)
factory.addPoint(-e3   , h1+r , 0, Lc2, 5)
factory.addPoint(-e3   , h1+h2, 0, Lc1, 6)
factory.addPoint( e3   , h1+h2, 0, Lc1, 7)
factory.addPoint( e3   , h1+r , 0, Lc2, 8)
factory.addPoint( e3+r , h1+r , 0, Lc2, 9)
factory.addPoint( e3+r , h1   , 0, Lc2, 10)
factory.addPoint( e1+e2, h1   , 0, Lc1, 11)
factory.addPoint( e1+e2, 0    , 0, Lc1, 12)
factory.addPoint( e2   , 0    , 0, Lc1, 13)

factory.addPoint( R1 / ssin, h5+R1*ccos, 0, Lc2, 14)
factory.addPoint( 0        , h5        , 0, Lc2, 15)
factory.addPoint(-R1 / ssin, h5+R1*ccos, 0, Lc2, 16)
factory.addPoint(-e2       , 0.0       , 0, Lc1, 17)

factory.addPoint(-R2 , h1+h3   , 0, Lc2, 18)
factory.addPoint(-R2 , h1+h3+h4, 0, Lc2, 19)
factory.addPoint( 0  , h1+h3+h4, 0, Lc2, 20)
factory.addPoint( R2 , h1+h3+h4, 0, Lc2, 21)
factory.addPoint( R2 , h1+h3   , 0, Lc2, 22)
factory.addPoint( 0  , h1+h3   , 0, Lc2, 23)

factory.addPoint( 0, h1+h3+h4+R2, 0, Lc2, 24)
factory.addPoint( 0, h1+h3-R2,    0, Lc2, 25)

factory.addLine(1 , 17, 1)
factory.addLine(17, 16, 2)

# Gmsh provides other curve primitives than straight lines: splines, B-splines,
# circle arcs, ellipse arcs, etc. Here we define a new circle arc, starting at
# point 14 and ending at point 16, with the circle's center being the point 15:
factory.addCircleArc(14,15,16, 3)

# Note that, in Gmsh, circle arcs should always be smaller than Pi. The
# OpenCASCADE geometry kernel does not have this limitation.

# We can then define additional lines and circles, as well as a new surface:
factory.addLine(14,13, 4)
factory.addLine(13,12, 5)
factory.addLine(12,11, 6)
factory.addLine(11,10, 7)
factory.addCircleArc(8,9,10, 8)
factory.addLine(8,7, 9)
factory.addLine(7,6, 10)
factory.addLine(6,5, 11)
factory.addCircleArc(3,4,5, 12)
factory.addLine(3,2, 13)
factory.addLine(2,1, 14)
factory.addLine(18,19, 15)
factory.addCircleArc(21,20,24, 16)
factory.addCircleArc(24,20,19, 17)
factory.addCircleArc(18,23,25, 18)
factory.addCircleArc(25,23,22, 19)
factory.addLine(21,22, 20)

# But we still need to define the exterior surface. Since this surface has a
# hole, its definition now requires two curves loops:
factory.addCurveLoop([17,-15,18,19,-20,16], 21)
factory.addPlaneSurface([21], 22)
factory.addCurveLoop([11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10], 23)

factory.addPlaneSurface([23,21], 24)

# As a general rule, if a surface has N holes, it is defined by N+1 curve loops:
# the first loop defines the exterior boundary; the other loops define the
# boundaries of the holes.
factory.synchronize()

# Finally, we can add some comments by creating a post-processing view
# containing some strings:
v = gmsh.view.add("comments")

# Add a text string in window coordinates, 10 pixels from the left and 10 pixels
# from the bottom:
gmsh.view.addListDataString(v, [10, -10], ["Created with Gmsh"])

# Add a text string in model coordinates centered at (X,Y,Z) = (0, 0.11, 0),
# with some style attributes:
gmsh.view.addListDataString(v, [0, 0.11, 0], ["Hole"],
                            ["Align", "Center", "Font", "Helvetica"])

# These annotations are handled by a list-based post-processing view. For
# large post-processing datasets, that contain actual field values defined on
# a mesh, you should use model-based post-processing views instead, which
# allow to efficiently store continuous or discontinuous scalar, vector and
# tensor fields, or arbitrary polynomial order.

# Views and geometrical entities can be made to respond to double-click
# events, here to print some messages to the console:
gmsh.option.setString("View[0].DoubleClickedCommand",
                      "Printf('View[0] has been double-clicked!');")
gmsh.option.setString(
    "Geometry.DoubleClickedLineCommand",
    "Printf('Curve %g has been double-clicked!', Geometry.DoubleClickedEntityTag);")

# We can also change the color of some entities:
gmsh.model.setColor([(2, 22)], 127, 127, 127)
gmsh.model.setColor([(2, 24)], 160, 32, 240)
gmsh.model.setColor([(1, i) for i in 1:14], 255, 0, 0)
gmsh.model.setColor([(1, i) for i in 15:20], 255, 255, 0)

gmsh.model.mesh.generate(2)

gmsh.write("t4.msh")

gmsh.finalize()


text:
"""
    forwardsub(L,b)

Solve the lower-triangular linear system with matrix `L` and
right-hand side vector `b`.
"""
function forwardsub(L,b)
    n = size(L,1)
    x = zeros(n)
    x[1] = b[1]/L[1,1]
    for i in 2:n
        s = sum( L[i,j]*x[j] for j in 1:i-1 )
        x[i] = ( b[i] - s ) / L[i,i]
    end
    return x
end

"""
    backsub(U,b)

Solve the upper-triangular linear system with matrix `U` and
right-hand side vector `b`.
"""
function backsub(U,b)
    n = size(U,1)
    x = zeros(n)
    x[n] = b[n]/U[n,n]
    for i in n-1:-1:1
        s = sum( U[i,j]*x[j] for j in i+1:n )
        x[i] = ( b[i] - s ) / U[i,i]
    end
    return x
end

"""
    lufact(A)

Compute the LU factorization of square matrix `A`, returning the
factors.
"""
function lufact(A)
    n = size(A,1)        # detect the dimensions from the input
    L = diagm(ones(n))   # ones on main diagonal, zeros elsewhere
    U = zeros(n,n)
    Aₖ = float(copy(A))  # make a working copy

    # Reduction by outer products
    for k in 1:n-1
        U[k,:] = Aₖ[k,:]
        L[:,k] = Aₖ[:,k]/U[k,k]
        Aₖ -= L[:,k]*U[k,:]'
    end
    U[n,n] = Aₖ[n,n]
    return LowerTriangular(L),UpperTriangular(U)
end

"""
    plufact(A)

Compute the PLU factorization of square matrix `A`, returning the
triangular factors and a row permutation vector.
"""
function plufact(A)
    n = size(A,1)
    L = zeros(n,n)
    U = zeros(n,n)
    p = fill(0,n)
    Aₖ = float(copy(A))

    # Reduction by outer products
    for k in 1:n-1
        p[k] = argmax(abs.(Aₖ[:,k]))
        U[k,:] = Aₖ[p[k],:]
        L[:,k] = Aₖ[:,k]/U[k,k]
        Aₖ -= L[:,k]*U[k,:]'
    end
    p[n] = argmax(abs.(Aₖ[:,n]))
    U[n,n] = Aₖ[p[n],n]
    L[:,n] = Aₖ[:,n]/U[n,n]
    return LowerTriangular(L[p,:]),U,p
end


text:
function read_records(baseName::String)
	#Helper function to read back in saved records

	return read_file(string(baseName, ".real")) + 1im*read_file(string(baseName, ".imag"))
end

function read_file(fileName::String)

	#Get the size of the file and open
	fileSize = filesize(fileName)
	FID = open(fileName, "r")

	#Read out the first three dimensions (recordLength x numWaveforms x numSegments)
	sizes = read(FID, Int32, 3)

	#Infer the number or records from the file size
	numRecords = ifloor((fileSize-3*sizeof(Int32))/prod(sizes)/sizeof(Float32))

	a = mmap_array(Float32, tuple([sizes, numRecords]...), FID)
	close(FID)
	return a
end