text:
using Test

@testset "transfinite rectangle" begin
    fname = tempname() * ".msh"
    nx, ny = 10, 7
    coefx, coefy = 1.2, 1.2
    @gmsh_do begin
        @addPoint begin
            -1.0, -1.0, 0.0, 0.0, 1
            1.0, -1.0, 0.0, 0.0, 2
            1.0, 1.0, 0.0, 0.0, 3
            -1.0, 1.0, 0.0, 0.0, 4
        end
        @addLine begin
            1, 2
            2, 3
            4, 3
            1, 4
        end
        gmsh.model.geo.addCurveLoop([1, 2, -3, -4], 1)
        gmsh.model.geo.addPlaneSurface([1], 1)
        gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])
        @setTransfiniteCurve begin
            1, nx+1, "Progression", coefx
            3, nx+1, "Progression", coefx
            2, ny+1, "Progression", coefy
            4, ny+1, "Progression", coefy
        end
        gmsh.model.geo.mesh.setRecombine(2, 1)
        gmsh.model.geo.synchronize()
        gmsh.model.mesh.generate(2)
        gmsh.write(fname)
    end

    els = @gmsh_open fname begin
        gmsh.model.mesh.getElements(2)
    end
    @test els[1][1] == 3
    @test length(els[2][1]) == nx * ny
    rm(fname)
end

@testset "adaptive mesh with field" begin
    fname = tempname() * ".msh"
    @gmsh_do begin
        factory = gmsh.model.geo
        @addPoint begin
            0.0, 0.0, 0.0, 0.0, 1
            1.0, 0.0, 0.0, 0.0, 2
            1.0, 1.0, 0.0, 0.0, 3
            0.0, 1.0, 0.0, 0.0, 4
        end
        @addLine begin
            1, 2, 1
            2, 3, 2
            3, 4, 3
            4, 1, 4
        end
        factory.addCurveLoop([1, 2, 3, 4], 1)
        factory.addPlaneSurface([1], 1)
        @addField 1 "Distance" begin
            "EdgesList", [2]
            "NNodesByEdge", 100
        end
        @addField 2 "Threshold" begin
            "IField", 1
            "LcMin", 0.005
            "LcMax", 0.15
            "DistMin", 0.15
            "DistMax", 0.50
        end
        @addField 3 "MathEval" begin
            "F", "0.01*(1.0+50.*(y-x*x)*(y-x*x) + (1-x)*(1-x))"
        end
        @addField 4 "Distance" begin
            "NodesList", [4]
        end
        @addField 5 "MathEval" begin
            "F", "F4^3 + 0.005"
        end
        @addField 6 "Min" begin
            "FieldsList", [2, 3, 5]
        end
        gmsh.model.mesh.field.setAsBackgroundMesh(6)
        gmsh.model.geo.synchronize()
        gmsh.model.mesh.generate(2)
        gmsh.write(fname)
    end

    els = @gmsh_open fname begin
        gmsh.model.mesh.getElements(2)
    end
    @test els[1][1] == 2 # 3-node triangle
    @test !isempty(els[2][1])
    rm(fname)
end

@testset "option of terminal output" begin
    _stdout = stdout
    rd, wr = redirect_stdout()
    @gmsh_do begin
        factory = gmsh.model.geo
        @addOption begin
            "General.Terminal", 1
            "Mesh.CharacteristicLengthMax", 0.4
            "Mesh.CharacteristicLengthMin", 0.1
        end

        @addPoint begin
            0.0, 0.0, 0.0, 0.1, 1
            0.0, 1.0, 0.0, 0.3, 2
            1.0, 1.0, 0.0, 0.2, 3
            1.0, 0.0, 0.0, 0.4, 4
        end
        @addLine begin
            1, 2, 1
            2, 3, 2
            3, 4, 3
            4, 1, 4
        end
        factory.addCurveLoop([1, 2, 3, 4], 1)
        factory.addPlaneSurface([1], 1)
        gmsh.model.geo.synchronize()
        gmsh.model.mesh.generate(2)
        gmsh.write("temp.msh")
    end
    redirect_stdout(_stdout)
    close(wr)
    output = read(rd, String)
    close(rd)
    @test occursin("Done meshing", output)
end


text:
# CUDA errors

const driver_error_descriptions = (Int=>ASCIIString)[
	0 => "Success",
	1 => "Invalid value",
	2 => "Out of memory",
	3 => "Driver not initialized",
	4 => "Driver being shutdown",
	5 => "Profiler disabled",
	6 => "Profiler not initialized",
	7 => "Profiler already started",
	8 => "Profiler already stopped",
	100 => "No CUDA-capable device",
	101 => "Invalid device ordinal",
	200 => "Invalid kernel image",
	201 => "Invalid context",
	202 => "Context already current",
	205 => "Map operation failed",
	206 => "Unmap operation failed",
	207 => "Array mapped",
	208 => "Resource already mapped",
	209 => "No kernel image available/suitable for GPU",
	210 => "Resource already acquired",
	211 => "Resource not mapped",
	212 => "Resource not mapped as array",
	213 => "Resource not mapped as pointer",
	214 => "Uncorrectable ECC error detected",
	215 => "Unsupported limit",
	216 => "Context already in use",
	217 => "Peer access not supported",
	300 => "Invalid kernel source",
	301 => "File not found",
	302 => "Shared object symbol not found",
	303 => "Shared object initialization failed",
	304 => "OS call failed",
	400 => "Invalid handle",
	500 => "Named symbol not found",
	600 => "Not ready",
	700 => "Kernel launch failed",
	701 => "Launch out of resources",
	702 => "Launch timeout",
	703 => "Incompatible texturing mode",
	704 => "Peer access already enabled",
	705 => "Peer access not enabled",
	708 => "Primary context already active",
	709 => "Context destroyed",
	710 => "Assertion triggered failure",
	711 => "Too many peers",
	712 => "Host memory already registered",
	713 => "Host memory not registered",
	800 => "Operation not permitted",
	801 => "Operation not supported",
	999 => "Unknown error"
]

immutable CuDriverError
	code::Int
end

description(err::CuDriverError) = driver_error_descriptions[err.code]

ext: jl


text:
module TimeWarp

using StatsBase
using Parameters

# reexport to let user specify dtw(⋅,⋅,distance)
using Reexport
@reexport using Distances

export Sequence,
       SequenceArray,
       seq_to_array,
       dtw,
       DTWDistance,
       DTWMethod,
       ClassicDTW,
       FastDTW,
       distpath,
       dtw_cost_matrix,
       fastdtw

include("utils.jl")
include("sequence.jl")
include("distance_interface.jl")
include("dtw.jl")
include("windowed_matrix.jl")
include("fastdtw.jl")

end # module


text:
using Base.Test
using EGR
using Redis

println("TestRedisRuns")

client = RedisConnection()

run(pipeline(`redis-cli keys "Test*"`, `xargs redis-cli del`))

createOracleOutputLevel = 1
numEquivalentPasses = 1
algOutputLevel = 2
maxOutputNum=5

myREfunction(problem, opts, sd, expIndices,n) = returnIfExists(client, problem, opts, sd, expIndices,n; outputLevel = 0)
myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = writeFunction(client, problem,  opts, sd,k, gnum, origWant, fromOutputsFunction; outputLevel = 0)

thisOracle=Oracles[1]

(gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) = thisOracle

maxG  = Int(round(numEquivalentPasses*numTP))

expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

myOutputOpts =  OutputOpts(myOutputter,expIndices)

thisSD = SGsd( "SG.a=1")

longKey = thisDataName*":"*LossFunctionString*":"*string(L2reg)*":"thisSD.stepString*":"*"const"*":0"
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=20, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=10, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts,thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=10, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts,thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=30, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=100000, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=100001, outputLevel=algOutputLevel)

alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=30, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

myOpts= Opts(zeros(numVars,1); stepSizePower=0, maxG=30, outputLevel=algOutputLevel)
alg(thisProblem(Task(() -> getRandom(numTP, gradientOracle)))	, myOpts, thisSD, myOutputOpts, myWriteFunction)
println(" gnum:     "*string([parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]))
println(" origWant: "*string([parse(Int64,s) for s = lrange(client, longKey*":origWant",0, -1)]))

println("TestRedisRuns successful!")
println()


text:
module FlightMechanics

using Markdown
using LinearAlgebra

export
# Air Constants
GAMMA_AIR, R_AIR, RHO0, P0, T0, A0,
# Earth constants
GRAVITY_ACCEL, STD_GRAV_PARAMETER, GRAV_CONSTANT, ROT_VELOCITY, EARTH_MASS,
EARTH_MEAN_RADIUS,
# Ellipsoids
Ellipsoid,
Clarke1866, Clarke1880,
International, Bessel, Everest, ModifiedEverest,
AustralianNational, SouthAmerican1969, Airy, ModifiedAiry, Hough,
Fischer1960SouthAsia, Fischer1960Mercury, Fischer1968,
WGS60, WGS66, WGS72, WGS84,
# Conversion constants
LB2KG, KG2LB, FT2M, M2FT, IN2M, M2IN, SLUG2KG, SLUGFT2_2_KGM2, SLUGFT3_2_KGM3,
KT2MS, MS2KT, DEG2RAD, RAD2DEG, HP2WAT, WAT2HP, LBF2N, PA2PSF, PSF2PA,
RANK2KEL, KEL2RANK,
# Coordinates
body2hor, rot_matrix_body2hor, hor2body, rot_matrix_hor2body,
wind2hor, hor2wind, body2wind, wind2body,
ecef2hor, rot_matrix_ecef2hor, hor2ecef, rot_matrix_hor2ecef,
body2ecef, rot_matrix_body2ecef, ecef2body, rot_matrix_ecef2body,
quaternion2euler, euler2quaternion,
llh2ecef, ecef2llh,
# atmosphere
atmosphere_isa, atmosphere_f16,
# anemometry
qc2cas, qc2tas, qc2eas,
tas2eas, eas2tas, cas2eas, eas2cas, cas2tas, tas2cas,
tas_alpha_beta_from_uvw,
incompressible_qinf, compressible_qinf,
# Mechanics
rigid_body_velocity, rigid_body_acceleration, steiner_inertia,
# Flight Mechanics
coordinated_turn_bank, climb_theta, turn_rate_angular_velocity,
body_angular_velocity_to_euler_angles_rates,
euler_angles_rates_to_body_angular_velocity,
body_angular_velocity_to_quaternion_rates,
uvw_to_tasαβ,
tas_α_β_dot_from_uvw_dot, uvw_dot_from_tas_α_β,
# Models
six_dof_euler_fixed_mass,
six_dof_quaternion_fixed_mass,
six_dof_aero_euler_fixed_mass,
six_dof_ecef_quaternion_fixed_mass

include("constants.jl")
include("atmosphere.jl")
include("coordinates.jl")
include("anemometry.jl")
include("mechanics.jl")
include("flight_mechanics.jl")

# 6 DOF dynamic models
include("dynamics/sixdof_euler_fixed_mass.jl")
include("dynamics/sixdof_quaternion_fixed_mass.jl")
include("dynamics/sixdof_ecef_quaternion_fixed_mass.jl")
include("dynamics/sixdof_aero_euler_fixed_mass.jl")

# Models
include("models/Models.jl")
# Aircrafts
include("aircrafts/Aircrafts.jl")

end # module


text:
# coding: utf-8

# In[1]:

#bitonic sort is my special feature in this repository!!!
#bitonic sort is cool as well as difficult to understand
#if you look for detailed explanation, you should check out geeksforgeeks
# https://www.geeksforgeeks.org/bitonic-sort/
#its python version is written in recursion
#which i found it really difficult to understand what the heck it was doing
#so i spent the whole day scratching my head to finish this non-recursive version
#a more straight forward explanation is the visualization from wikipedia
# https://upload.wikimedia.org/wikipedia/commons/thumb/b/bd/BitonicSort1.svg/843px-BitonicSort1.svg.png

#enjoy~~~

# In[2]:

function bitonic_sort(arr)

    if length(arr)==1

        return arr

    end

    i=1

    #primary iteration
    while 2^(i-1)<length(arr)/2+1

        j=2^i

        direction=get_direction(j,length(arr))

        #secondary iteration
        while j>=2

            #double loop helps us to find out who is comparing with who for each round
            #the sorting order is given by another function called get_direction
            for k in 1:j:length(arr)

                #k+Int32(j/2)-1 is crucial
                #make sure it doesnt go out of bound
                for l in k:k+Int32(j/2)-1

                    if (arr[l]>arr[l+Int32(j/2)]) && (direction[l]==1)

                        arr[l],arr[l+Int32(j/2)]=arr[l+Int32(j/2)],arr[l]

                    end

                    if (arr[l]<arr[l+Int32(j/2)]) && (direction[l]==-1)

                        arr[l],arr[l+Int32(j/2)]=arr[l+Int32(j/2)],arr[l]

                    end

               end

            end

            #Int32 is important
            #list indices dont support float
            #when u do /, we get float by default
            j=Int32(j/2)

        end

        i+=1

    end

    return arr

end

# In[3]:

#the annoying part of bitonic sort is that
#some parts of the sequence are sorted by descending order
#other parts of the sequence are sorted by ascending order
#for each primary iteration, we create a sorting order sequence for each element
#for each secondary iteration inside a primary iteration
#the sorting order for each element is the same
function get_direction(iteration_number,len)

    direction=[0 for _ in 1:len]

    ind=[i for i in 1:iteration_number:len]

    for i in ind[1:2:length(ind)]

        direction[i]=1

    end

    for i in ind[2:2:length(ind)]

        direction[i]=-1

    end

    for i in 1:length(direction)

        if direction[i]==0

            direction[i]=direction[i-1]

        end

    end

    return direction

end

# In[4]:

#bitonic sort only works on a list with 2^n elements
for _ in 1:100

    test_arr=rand(2^rand(1:10))

    if !(bitonic_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: bitonic sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: bitonic sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: bitonic sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
abstract type BehaviorModel end
abstract type AbstractMLRewardModel end
abstract type AbstractMLDynamicsModel end

mutable struct MLMDP{S, A, DModel<:AbstractMLDynamicsModel, RModel<:AbstractMLRewardModel}  <: MDP{S, A}
    dmodel::DModel
    rmodel::RModel
    discount::Float64
end

mutable struct MLPOMDP{S, A, O, DModel<:AbstractMLDynamicsModel, RModel<:AbstractMLRewardModel}  <: POMDP{S, A, O}
  dmodel::DModel
  rmodel::RModel
  discount::Float64
end

struct OriginalRewardModel <: AbstractMLRewardModel
	r_crash::Float64
	accel_cost::Float64
	decel_cost::Float64
	invalid_cost::Float64
	lineride_cost::Float64
	lanechange_cost::Float64
end

struct IDMMOBILModel <: AbstractMLDynamicsModel
	nb_cars::Int
    phys_param::PhysicalParam

	BEHAVIORS::Array{BehaviorModel,1}
	NB_PHENOTYPES::Int

	encounter_prob::Float64
	accels::Array{Int,1}
end

function IDMMOBILModel(nb_cars, phys_param; encounter_prob=0.5, accels=Int[-3,-2,-1,0,1])
    BEHAVIORS = IDMMOBILBehavior[IDMMOBILBehavior(x[1],x[2],x[3],idx) for (idx,x) in enumerate(Iterators.product(["cautious","normal","aggressive"],[phys_param.v_slow+0.5;phys_param.v_med;phys_param.v_fast],[phys_param.l_car]))]
    return IDMMOBILModel(nb_cars, phys_param, BEHAVIORS, length(BEHAVIORS), encounter_prob, accels)
end

# TODO for performance, parameterize this by BehaviorModel
struct CarState
    x::Float64
    y::Float64
	vel::Float64 #v_x
	lane_change::Float64 # ydot # in units of LANES PER SECOND
	behavior::BehaviorModel
    id::Int # car id to track from state to state - ego is ALWAYS 1
end

function ==(a::CarState,b::CarState)
    return a.x==b.x && a.y==b.y && a.vel==b.vel && a.lane_change == b.lane_change && a.id == b.id && a.behavior==b.behavior
end
Base.hash(a::CarState, h::UInt64=zero(UInt64)) = hash(a.vel, hash(a.x, hash(a.y, hash(a.lane_change, hash(a.behavior, hash(a.id, h))))))
"Return a representation that will produce a valid object if executed"
Base.repr(c::CarState) = "CarState($(c.x),$(c.y),$(c.vel),$(c.lane_change),$(c.behavior),$(c.id))"

mutable struct MLState
    x::Float64 # total distance traveled by the ego
    t::Float64 # total time of the simulation
	cars::Array{CarState,1} #NOTE ego car is first car
    Control_Signal::Bool  # If we use control Logic
    terminal::Nullable{Any} # SHOULD BE Nullable{Symbol} if not null, this is a terminal state, see below

    #Target_achieve ::Bool
end
# more constructors at bottom

#=
Terminal states: Each terminal state is not considered different, ther terminal states are
    :crash
    :lane
    :brake
    :distance
=#

function ==(a::MLState, b::MLState)
    if isnull(a.terminal) && isnull(b.terminal) # neither terminal
        return a.x == b.x && a.t == b.t && a.cars == b.cars
    elseif !isnull(a.terminal) && !isnull(b.terminal) # both terminal
        return get(a.terminal) == get(b.terminal)
    else # one is terminal
        return false
    end
end
function Base.hash(a::MLState, h::UInt64=zero(UInt64))
    if isnull(a.terminal)
        return hash(a.x, hash(a.t, hash(a.cars,h)))
    else
        return hash(get(a.terminal), h)
    end
end

struct MLAction
    acc::Float64
    lane_change::Float64 # ydot
end

#MLAction(acc::Float64,lane_change::Float64) = MLAction(acc,lane_change,false)
MLAction() = MLAction(0.0,0.0)
==(a::MLAction,b::MLAction) = (a.acc==b.acc) && (a.lane_change==b.lane_change)
Base.hash(a::MLAction,h::UInt64=zero(UInt64)) = hash(a.acc,hash(a.lane_change,h))
function MLAction(x::Array{Float64,1})
	assert(length(x)==2)
	lane_change = abs(x[2]) <= 0.3? 0: sign(x[2])
	return MLAction(x[1],lane_change)
end
vec(a::MLAction) = Float64[a.acc;a.lane_change]

const OriginalMDP=MLMDP{MLState, MLAction, IDMMOBILModel, OriginalRewardModel}

struct ActionSpace
	actions::Vector{MLAction}
end

struct CarPhysicalState
    x::Float64
    y::Float64
    vel::Float64
    lane_change::Float64
    id::Int
end
const CarStateObs = CarPhysicalState

==(a::CarPhysicalState, b::CarPhysicalState) = (a.x == b.x) && (a.y == b.y) && (a.vel == b.vel) && (a.lane_change == b.lane_change) && (a.id == b.id)
Base.hash(a::CarPhysicalState, h::UInt64=zero(UInt64)) = hash(a.x, hash(a.y, hash(a.vel, (hash(a.lane_change, hash(a.id,h))))))
CarPhysicalState(cs::CarState) = CarPhysicalState(cs.x, cs.y, cs.vel, cs.lane_change, cs.id)
function CarState(cps::CarPhysicalState, behavior::BehaviorModel)
    return CarState(cps.x, cps.y, cps.vel, cps.lane_change, behavior, cps.id)
end

struct MLPhysicalState
    x::Float64
    t::Float64
    cars::Array{CarPhysicalState,1}
    terminal::Nullable{Any} # Should be Nullable{Symbol}
end
const MLObs = MLPhysicalState

MLPhysicalState(s::MLState) = MLPhysicalState(s.x, s.t, CarPhysicalState[CarPhysicalState(cs) for cs in s.cars], s.terminal)

function ==(a::MLPhysicalState, b::MLPhysicalState)
    if isnull(a.terminal) && isnull(b.terminal) # neither terminal
        return a.x == b.x && a.t == b.t && a.cars == b.cars
    elseif !isnull(a.terminal) && !isnull(b.terminal) # both terminal
        return get(a.terminal) == get(b.terminal)
    else # one is terminal
        return false
    end
end
function Base.hash(a::MLPhysicalState, h::UInt64=zero(UInt64))
    if isnull(a.terminal)
        return hash(a.x, hash(a.t, hash(a.cars,h)))
    else
        return hash(get(a.terminal), h)
    end
end

MLState(ps::MLPhysicalState, cars::Vector{CarState}) = MLState(ps.x, ps.t, cars, ps.terminal)
MLState(s::MLState, cars::Vector{CarState}) = MLState(s.x, s.t, cars, s.Control_Signal, s.terminal)
MLState(x::Float64, t::Float64, cars::Vector{CarState}, Control_Signal::Bool) = MLState(x, t, cars, Control_Signal, nothing)

#############################################

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
function get_robust_path_two(start_node, end_node, p, q, c, d, Gamma_u, Gamma_v, origin, destination)
  @assert length(start_node)==length(end_node)
  @assert length(start_node)==length(p)
  @assert length(start_node)==length(q)
  @assert length(start_node)==length(c)
  @assert length(start_node)==length(d)

  no_arc = length(start_node)

  # Step 0
  diag = p.*q + q.*c + q.*d
  idx = sortperm(diag, rev=true)
  # Theta = Array((Float64, Float64),1)
  # Theta = Array{Tuple{Float64, Float64}}(1)
  Theta = [(0.0,0.0)]

  for t=1:no_arc
    ij = idx[t]

    # Step 1
    Theta = [Theta; (0, p[ij]*d[ij]+q[ij]*d[ij])]

    # Step 2
    Ay  = p[ij]*d[ij] + q[ij]*d[ij]
    Bx  = q[ij]*c[ij]
    Cxy = p[ij]*d[ij] + q[ij]*c[ij] + q[ij]*d[ij]
    Dx  = q[ij]*c[ij] + q[ij]*d[ij]
    Ey  = p[ij]*d[ij]

    for s=t+1:no_arc
      kl = idx[s]

      AAy  = p[kl]*d[kl] + q[kl]*d[kl]
      BBx  = q[kl]*c[kl]
      CCxy = p[kl]*d[kl] + q[kl]*c[kl] + q[kl]*d[kl]
      DDx  = q[kl]*c[kl] + q[kl]*d[kl]
      EEy  = p[kl]*d[kl]

      # Pattern 1
      if Ay>AAy && Bx>BBx && Ey <EEy
        Theta = [Theta; (BBx, Ay); (Cxy-EEy, EEy)]
      end

      # Pattern 2
      if Ay>AAy && Bx>BBx && Ey>EEy && Dx>DDx
        Theta = [Theta; (BBx, Ay); (Dx, EEy)]
      end

      # Pattern 3
      if Ay>AAy && Bx>BBx && Dx<DDx && Ey>EEy
        Theta = [Theta; (BBx, Ay); (Dx, CCxy-Dx)]
      end

      # Pattern 4
      if EEy<Ay && Ay<AAy && Bx>BBx && Ey<EEy
        Theta = [Theta; (CCxy-Ay,Ay); (Cxy-EEy, EEy)]
      end

      # Pattern 5
      if EEy<Ay && Ay<AAy && Bx>BBx && Ey>EEy && Dx > DDx
        Theta = [Theta; (CCxy-Ay,Ay); (Dx, EEy)]
      end

      # Pattern 6
      if EEy<Ay && Ay<AAy && Bx>BBx && Ey>EEy && Dx < DDx
        Theta = [Theta; (CCxy-Ay,Ay); (Dx, CCxy-Dx)]
      end

      # Pattern 7
      if Bx>DDx && Ay<EEy
        Theta = [Theta; (DDx, Ay); (Bx, EEy)]
      end

      # Pattern 8
      if Bx<BBx && Ey<EEy
        Theta = [Theta; (BBx, Cxy-BBx); (Cxy-EEy, EEy)]
      end

      # Pattern 9
      if Bx<BBx && Ey>EEy && Dx>DDx
        Theta = [Theta; (BBx, Cxy-BBx); (Dx, EEy)]
      end

      # Pattern 10
      if Bx<BBx && Ey>EEy && Dx<DDx
        Theta = [Theta; (BBx, Cxy-BBx); (Dx, CCxy-Dx)]
      end

      # Pattern 11
      if Ey>AAy && Dx<BBx
        Theta = [Theta; (Dx, AAy); (BBx, Ey)]
      end
    end
  end
  # The set Theta contains all necessary (theta_u, theta_v)

  theta_u = zeros(size(Theta))
  theta_v = zeros(size(Theta))
  z_front = zeros(size(Theta))
  for i=1:length(Theta)
    theta_u[i] = Theta[i][1]
    theta_v[i] = Theta[i][2]
    z_front[i] = Gamma_u*theta_u[i] + Gamma_v*theta_v[i]
  end

  z_idx = sortperm(z_front, rev=true)

  best_obj = Inf
  best_path = []
  best_x = []

  for i=1:length(Theta)
    k = z_idx[i]
    if Gamma_u*theta_u[k] + Gamma_v*theta_v[k] >= best_obj
        break
    end

    rho_mu = zeros(no_arc)
    for ij=1:no_arc
        qc = q[ij]*c[ij]
        pd = p[ij]*d[ij]
        qd = q[ij]*d[ij]

        Condition1 = theta_u[k]>=qc && theta_v[k]>=pd+qd
        Condition2 = theta_u[k]<=qc && theta_v[k]>=pd+qd
        Condition3 = pd<=theta_v[k] && theta_v[k]<=pd+qd && theta_u[k]+theta_v[k]<=pd+qc+qd
        Condition4 = pd<=theta_v[k] && theta_v[k]<=pd+qd && theta_u[k]+theta_v[k]>=pd+qc+qd
        Condition5 = theta_u[k]<=qc+qd && theta_v[k]<=pd
        Condition6 = theta_u[k]>=qc+qd && theta_v[k]<=pd

        if Condition1 || Condition4
          rho_mu[ij] = 0
        elseif Condition2
          rho_mu[ij] = qc-theta_u[k]
        elseif Condition3 || Condition5
          rho_mu[ij] = pd+qc+qd-theta_u[k]-theta_v[k]
        elseif Condition6
          rho_mu[ij] = pd-theta_v[k]
        end
    end

    link_cost = p.*c + rho_mu
    (pp, xx, ll) = get_shortest_path(start_node, end_node, link_cost, origin, destination)
    # current_obj = Gamma_u*theta_u[k] + Gamma_v*theta_v[k] + dot(link_cost, xx)
    current_obj = Gamma_u*theta_u[k] + Gamma_v*theta_v[k] + ll

    if current_obj < best_obj
        best_obj = current_obj
        best_path = pp
        best_x = xx
    end

  end

  return (best_path, best_x, best_obj)

end


text:
result = begin
    a = 10
    b = 20
    a + b
end

println(result)

result = (a = 10; b = 20; c = 30; a+b+c)

println(result)

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
function Plots._apply_recipe(d::Dict, plm::PairwiseListMatrix; kw...)

    get!(d, :yflip, true)
    get!(d, :linetype, :heatmap)

    # args
    (plm, )
end

function arcdiagram(plm::PairwiseListMatrix)

  n = plm.nelements
  x = 1:n

  plt = scatter(x, zeros(n), m=(10,:orange), leg=false)

  zmin,zmax = extrema(plm)
  grad = ColorGradient(:bluesreds)
  curvecolor(v) = getColorZ(grad, (v-zmin)/(zmax-zmin))

  for i in 1:(n-1)
    for j in i:n
                                                       # i-j
      curve = BezierCurve(P2[(x[i],0), ((x[i]+x[j])/2, mat[i,j]), (x[j], 0)])
      plot!(curve_points(curve), line = (curvecolor(mat[i,j]), 0.5, 2))
    end
  end

  plt
end


text:
using OmegaCore, OmegaGrad, Distributions, Test

function zygotetest_1()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Normal(0, 1)
  z(ω) = x(ω) + y(ω)
  ω = defω()
  x_ = x(ω)
  y_ = y(ω)
  @test grad(z, ω, ZygoteGrad) = (1.0, 1.0)
end

function zygotetest_2()
  x = [i ~ Normal(0, 1) for i = 1:5]
  y(ω) = sum(x(ω))
  ω = defω()
  y_ = y(ω)
  @test grad(y, ω, ZygoteGrad) = (1.0,1.0,1.0,1.0,1.0)
end

function zygotetest_3()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Normal(0, 1)
  z(ω) = x(ω) * y(ω)
  ω = defω()
  x_ = x(ω)
  y_ = y(ω)
  @test grad(z, ω, ZygoteGrad) = (ω.data[[2]][2], ω.data[[1]][2])
end

function zygotetest_4()
  x = 1 ~ Normal(0, 1)
  z(ω) = x(ω) * x(ω) * x(ω)
  ω = defω()
  x_ = x(ω)
  @test grad(z, ω, ZygoteGrad) = (3.0 * (ω.data[[1]][2] ^ 2),)
end

function zygotetest_5()
  x = 1 ~ Normal(0, 1)
  z(ω) = sin(x(ω))
  ω = defω()
  x_ = x(ω)
  @test grad(z, ω, ZygoteGrad) = (cos(ω.data[[1]][2]),)
end

@testset "Gradients" begin
  zygotetest_1()
  zygotetest_2()
  zygotetest_3()
  zygotetest_4()
  zygotetest_5()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Test

using Jecco
using Jecco.AdS5_3_1

@testset "All tests:" begin

    @time @testset "derivative operators tests:" begin include("deriv_operators_tests.jl") end

    @time @testset "grid tests:" begin include("grid_tests.jl") end

    @time @testset "output tests:" begin include("output_tests.jl") end

    @time @testset "filtering tests:" begin include("filtering_tests.jl") end

    @time @testset "AdS5_3_1 inner grid coefficients tests:" begin
        include("AdS5_3_1_equation_inner_coeff.jl")
    end

    @time @testset "AdS5_3_1 outer grid coefficients tests:" begin
        include("AdS5_3_1_equation_outer_coeff.jl")
    end

    @time @testset "AdS5_3_1 homogeneous black brane tests:" begin
        include("AdS5_3_1_homog_Bbrane.jl")
    end

end

nothing


text:
#!/usr/bin/env julia
export @include, @genNew, @updateAll

primitive_types = Dict(
    "bool"      => "Bool",
    "int8"      => "Int8",
    "byte"      => "Int8",
    "uint8"     => "UInt8",
    "char"      => "UInt8",
    "int16"     => "Int16",
    "uint16"    => "UInt16",
    "int32"     => "Number",
    "uint32"    => "Unsigned",
    "int64"     => "Number",
    "uint64"    => "Unsigned",
    "float32"   => "Number",
    "float64"   => "Number",
    "string"    => "String",
    "bool[]"    => "Union{Vector{Bool}, Array{Bool}}",
    "int8[]"    => "Union{Vector{Int8}, Array{Int8}}",
    "byte[]"    => "Union{Vector{Int8}, Array{Int8}}",
    "uint8[]"   => "Union{Vector{UInt8}, Array{UInt8}}",
    "char[]"    => "Union{Vector{UInt8}, Array{UInt8}}",
    "int16[]"   => "Union{Vector{Int16}i, Array{UInt16}}",
    "uint16[]"  => "Union{Vector{UInt16}, Array{UInt16}}",
    "int32[]"   => "Union{Vector{<:Number}, Array{<:Number}}",
    "uint32[]"  => "Union{Vector{<:Number}, Array{<:Number}}",
    "int64[]"   => "Union{Vector{<:Number}, Array{<:Number}}",
    "uint64[]"  => "Union{Vector{<:Number}, Array{<:Number}}",
    "float32[]" => "Union{Vector{<:Number}, Array{<:Number}}",
    "float64[]" => "Union{Vector{<:Number}, Array{<:Number}}",
    "string[]"  => "Union{Vector{String}, Array{String}}",
    "time"      => "Union{Cxx.CxxCore.CppPtr, Cxx.CxxCore.CppValue, Cxx.CxxCore.CppRef}",
    "duration"  => "Union{Cxx.CxxCore.CppPtr, Cxx.CxxCore.CppValue, Cxx.CxxCore.CppRef}",
    "Header"    => "Union{Cxx.CxxCore.CppPtr, Cxx.CxxCore.CppValue, Cxx.CxxCore.CppRef}",
)

"""
    include(expr)

Inlude msg and srv files from generated ROS packages.
"""
macro include(expr)
    n = replace(string(expr), " " => "")
    n = replace(n, "(" => "")
    n = replace(n, ")" => "")
    toinclude = split(n,':')
    if length(toinclude) != 2
        throw(ErrorException("Error while trying to read @rosinclude $n"))
    end
    pkg = string(toinclude[1])
    for file in split(toinclude[2],',')
        f = joinpath(@__DIR__,"gen", "$(pkg)_$(file).jl")
        if isfile(f)
            include(f)
        else
            throw(ErrorException("Error while trying to include $f. Try running the @ROS.genNew macro."))
        end
    end
end

function scout(pkg, f)
    package_dir = package.getPath(pkg)
    # NOTE not finding the package path is
    # not handled because the function
    # is guaranteed to be called when the
    # package is found
    if length(package_dir) > 0
        additional_folder = ""
        if occursin("/src/", package_dir) # the package was found in a catkin workspace, and not inside the ROS installation (/opt/ros/...)
            additional_folder = "devel"
            found_ws_root = false
            ws_root = package_dir
            while !found_ws_root
                ws_root = normpath(ws_root, "..")
                found_ws_root = !occursin("src", ws_root)
            end
            package_dir = joinpath(ws_root, additional_folder, "include", pkg, "$(f).h")
        else
            package_dir = joinpath(normpath(package_dir, "../.."), additional_folder, "include", pkg, "$(f).h")
        end
    end
    return package_dir
end

function startGen(pkg, file, type, fullpath, overwrite)
    jl_fname = "$(pkg)_$(file).jl"
    target = joinpath(@__DIR__,"gen",jl_fname)
    if ! overwrite && isfile(target)
        println(target," found but already generated. Skipping.")
    else
        thefile = joinpath(fullpath, file*"."*type)
        includes, types, fields, prefix, suffix, constants = fileAnalysis(thefile, pkg, type)
        generated = finishGen(target, pkg, file, type, includes, types, fields, prefix, suffix, constants)
        open(target, "w") do f
            write(f, generated)
        end
    end
end

function finishGen(jl_fname, pkg, name, type, includes, types, fields, prefix, suffix, constants)
    header_location = ""
    if type == "action"
        header_location = scout(pkg,name*"Action")
    else
        header_location = scout(pkg,name)
    end
    cxxinclude(header_location)
    meta = "# This file was auto-generated using typegen.jl.\n"
    meta *= "include(\"__ROSTYPE.jl\")\n"
    meta *= "cxxinclude(\"$header_location\")\n"
    for inc in includes
        inc = replace(inc, "[]"=>"")
        meta *= "cxxinclude(\"$inc\")\n"
        jl_inc = String.(split(inc, '/'))
        jl_inc1 = last(jl_inc)
        jl_inc1 = replace(jl_inc1, ".hpp"=>"")
        jl_inc1 = replace(jl_inc1, ".h"=>"")
        jl_inc2 = jl_inc[length(jl_inc)-1]
        meta *= "include(\"$(jl_inc2)_$(jl_inc1).jl\")\n"
    end
    meta *= "\n\n"

    # Minimize meta-calls
    secondary_types::Dict{String, String} = Dict()
    for s in suffix
        if !isempty(s) && !haskey(secondary_types, s)
            ptr_value = string(eval(Meta.parse("typeof(@cxxnew $pkg::$name::$s())")))
            ptr = string(eval(Meta.parse("cxxt\"$pkg::$name::$s*\"")))
            const_ptr = string(eval(Meta.parse("cxxt\"const $pkg::$name*\"")))
            secondary_types[s] = "Union{$ptr_value, $ptr, $const_ptr}"
        end
    end

    ptr_value = string(eval(Meta.parse("typeof(@cxxnew $pkg::$name())")))
    ptr = string(eval(Meta.parse("cxxt\"$pkg::$name*\"")))
    const_ptr = string(eval(Meta.parse("cxxt\"const $pkg::$name*\"")))
    main_type = "Union{$ptr_value, $ptr, $const_ptr}"

    for s in unique(suffix)
        if !isempty(s)
            meta *= "const $(pkg)_$(name)_$s = cxxt\"$pkg::$name::$s\"\n"
            meta *= "$(pkg)_$(name)_$s() = @cxxnew $pkg::$name::$s()\n\n"
        end
    end
    meta *= "const $(pkg)_$(name) = cxxt\"$pkg::$name\"\n"
    meta *= "$(pkg)_$name() = @cxxnew $pkg::$name()\n\n"

    for constant in constants
        meta *= "function getWrap(c::$main_type, f::Val{:$(constant[2])})\n"
        meta *= "    return  "
        if constant[1] == "Number"
            meta *= "parse(Float64, \"$(constant[3])\")\n"
        elseif constant[1] == "Unsigned"
            meta *= "parse(UInt64, \"$(constant[3])\")\n"
        elseif constant[1] == "String"
            meta *= "\"$(constant[3])\"\n"
        else
            meta *= "parse($(constant[1]), \"$(constant[3])\")\n"
        end
        meta *= "end\n\n"
    end

    for i in 1:length(types)
        if !isempty(suffix[i])
            meta *= "function getWrap(c::$(secondary_types[suffix[i]]), f::Val{:$(fields[i])})\n"
        else
            meta *= "function getWrap(c::$(main_type), f::Val{:$(fields[i])})\n"
        end
        meta *= "    "
        # Primitive type
        if types[i] in values(primitive_types)
            # Special handling for string
            if types[i] == "String"
                meta *= "unsafe_string(@cxx c->$(fields[i]))\n"
            # Vector of primitives
            elseif occursin("Vector{", types[i])
                meta *= "unsafe_wrap(DenseArray, @cxx c->$(fields[i]))\n"
            # Primitive but not built-in types
            elseif startswith(types[i], "Union{")
                meta *= "@cxx &c->$(fields[i])\n"
            # Rest of built-in types
            else
                meta *= "@cxx c->$(fields[i])\n"
            end
        # Vector of misc message type
        elseif types[i] == "Union{Vector, Array}"
            meta *= "unsafe_wrap(DenseArray, @cxx c->$(fields[i]))\n"
        # Misc message type
        else
            meta *= "@cxx &c->$(fields[i])\n"
        end
        meta *= "end\n\n"
    end

    for s in unique(suffix)
        if !isempty(s)
            meta *= "function getWrap(c::$(secondary_types[s]), f)\n"
            meta *= "    throw(ErrorException(\"$pkg::$(name)$s does not have a field of \$(f)!\"))\nend\n\n"

            meta *= "function Base.getproperty(c::$(secondary_types[s]), f::Symbol)\n"
            meta *= "    getWrap(c, Val(f))\n"
            meta *= "end\n\n"
        end
    end

    for p in prefix
        meta *= "function getWrap(c::$main_type, f::Val{:$p})\n"
        meta *= "    @cxx &c->$p\n"
        meta *= "end\n\n"
    end

    meta *= "function getWrap(c::$main_type, f)\n"
    meta *= "    throw(ErrorException(\"$pkg::$name does not have a field of \$(f)!\"))\nend\n\n"

    meta *= "function Base.getproperty(c::$main_type, f::Symbol)\n"
    meta *= "    getWrap(c, Val(f))\n"
    meta *= "end\n\n"

    for i in 1:length(types)
        if !isempty(suffix[i])
            meta *= "function setWrap!(c::$(secondary_types[suffix[i]]), f::Val{:$(fields[i])}, value::$(types[i]))\n"
        else
            meta *= "function setWrap!(c::$(main_type), f::Val{:$(fields[i])}, value::$(types[i]))\n"
        end
        meta *= "    "
        # Primitive type
        if types[i] in values(primitive_types)
            # Special handling for string
            if types[i] == "String"
                meta *= "icxx\" \$c->$(fields[i]) = \$(pointer(value));\"\n"
            elseif occursin("Vector{", types[i])
                meta *= "icxx\" \$c->$(fields[i]).clear();\"\n"
                meta *= "    for i in 1:length(value)\n"
                meta *= "        icxx\" \$c->$(fields[i]).push_back(\$(value[i]));\"\n"
                meta *= "    end\n"
            # Primitive but not built-in types
            elseif startswith(types[i], "Union{")
                meta *= "icxx\" \$c->$(fields[i]) = \$value;\"\n"
            # Rest of built-in types
            else
                meta *= "icxx\" \$c->$(fields[i]) = \$value;\"\n"
            end
        # Vector of misc message type
        elseif types[i] == "Union{Vector, Array}"
            meta *= "icxx\" \$c->$(fields[i]).clear();\"\n"
            meta *= "    for i in 1:length(value)\n"
            meta *= "        icxx\" \$c->$(fields[i]).push_back(*\$(value[i]));\"\n"
            meta *= "    end\n"
        # Misc message type
        else
            meta *= "icxx\" \$c->$(fields[i]) = \$value;\"\n"
        end
        meta*= "end\n\n"
    end

    for s in unique(suffix)
        if !isempty(s)
            meta *= "function setWrap!(c::$(secondary_types[s]), f, value)\n"
            meta *= "    throw(ErrorException(\"$pkg::$name does not have a field of \$(f)!\"))\nend\n\n"

            meta *= "function Base.setproperty!(c::$(secondary_types[s]), f::Symbol, value)\n"

            meta *= "    setWrap!(c, Val(f), value)\n"
            meta *= "end\n\n"
        end
    end

    meta *= "function setWrap!(c::$main_type, f, value)\n"
    meta *= "    throw(ErrorException(\"$pkg::$name does not have a field of \$(f)!\"))\nend\n\n"

    meta *= "function Base.setproperty!(c::$main_type, f::Symbol, value)\n"

    meta *= "    setWrap!(c, Val(f), value)\n"
    meta *= "end\n\n"

    return meta
end

function fileAnalysis(file, parent_pkg, t)
    p = ["", "request", "response", "goal", "result", "feedback"]
    s = ["", "Request", "Response", "Goal", "Result", "Feedback", "ActionGoal", "ActionResult", "ActionFeedback"]
    includes::Vector{String} = []
    types::Vector{String} = []
    fields::Vector{String} = []
    # Constants contain type,field,value
    constants::Vector{Tuple{String, String, String,String,String}} = []
    # The prefix vector is used for services and actions
    # to include one of the values of p
    prefix::Vector{String} = []
    # The suffix vector is used for services and actions
    # to include one of the values of s
    suffix::Vector{String} = []
    pindex = 1
    if t == "srv"
        pindex += 1
        push!(prefix, p[pindex])
    elseif t == "action"
        pindex += 3
        push!(prefix, p[pindex])
    end
    for line in eachline(file)
        line = strip(replace(line, r"#.*" => ""))
        # Not a comment
        if !isempty(line)
            # Change of service/action definition
            if line == "---"
                pindex += 1
                push!(prefix, p[pindex])
            # Not a constant
            elseif !('=' in line)
                type,field = String.(split(line))
                if type == "Header"
                    type = "std_msgs/Header"
                end
                if haskey(primitive_types, type)
                    type = primitive_types[type]
                elseif haskey(primitive_types, replace(type, r"\[(.*)\]" => ""))
                    type = primitive_types[replace(type, r"\[(.*)\]" => "")]
                else
                    if !occursin("/", type)
                        type = parent_pkg * "/" * type
                    end
                    pkg, f = String.(split(type, '/'))
                    push!(includes, scout(pkg,f))
                    type = replace(type, '/' => '.')
                end
                if occursin(".", type)
                    if endswith(type, r"\[(.*)\]")
                        type = "Union{Vector, Array}"
                    else
                        type = "Union{Cxx.CxxCore.CppPtr, Cxx.CxxCore.CppValue, Cxx.CxxCore.CppRef}"
                    end
                end
                push!(types, type)
                push!(fields, field)
                push!(suffix, s[pindex])
                #  if t == "action"
                    #  push!(types, type)
                    #  push!(fields, field)
                    #  push!(suffix, "Action"*s[pindex])
                #  end
            else
                typefield, value = String.(split(line, '='))
                type, field = String.(split(typefield))
                push!(constants, (type,field,value,p[pindex],s[pindex]))
            end
        end
    end
    return includes, types, fields, prefix, suffix, constants
end

function msgGenerator(overwrite=false)
    println("Generating ROS types...")
    #  types_of_interest = ["msg", "srv", "action"]
    types_of_interest = ["msg", "srv"]
    packages = package.getAll()
    for p in packages
        pp = package.getPath(p)
        for t in types_of_interest
            d = joinpath(pp, t)
            if isdir(d)
                for file in readdir(d)
                    suffix_pos = findlast(isequal('.'),file)
                    if suffix_pos != nothing && file[nextind(file,suffix_pos):end] == t
                        startGen(p, file[1:prevind(file, suffix_pos)], t, d, overwrite)
                    end
                end
            end
        end
    end
end

"""
    genNew()

Initiate the generation of Julia code for all newly found msg and srv ROS packages.
"""
macro genNew()
    msgGenerator(false)
end

"""
    updateAll()

Initiate the generation of Julia code for all msg and srv ROS packages (new and already existing ones).
"""
macro updateAll()
    msgGenerator(true)
end


text:
module jMOR

import Base: .*
using PyCall
using JLD, HDF5
using Base.SparseArrays
unshift!(PyVector(pyimport("sys")["path"]), "")

export rbfile, fem, jldimport, @femimport, load

# import python fem solver #
macro femimport(eq, eqname...)
    eval(:(@pyimport $eq as fem))

    len = length(eqname)

    if len > 0 && (len != 2 || eqname[1] != :as)
        throw(ArgumentError("The correct use is @femimport eq or @femimport eq as eqname"))
    elseif len == 2
        name = string(eqname[2])
    else
        name = string(eq)
    end

    try
        jldimport(name, "r+")
    catch
        jldimport(name, "w")
    end

    clean()
end

# import jld file with data from the RB problem #
function jldimport(name::AbstractString, access::AbstractString)
    global rbfile = jldopen(name * ".jld", access)

    # initialize all the variables
    if access == "r+"
        for var in map(x -> Symbol(x), names(rbfile))
            q = quote
                global $var = 0
                export $var
            end

            eval(q)
        end
    end

    nothing
end

export params, S, ns, sigma, sigma2, N, N2, Nh, Afem, bfem, U, U2, Xh, r

function clean()
    global params = 0
    global S = 0
    global ns = 0
    global sigma = 0
    global sigma2 = 0
    global N = 0
    global N2 = 0
    global Nh = 0
    global Afem = 0
    global bfem = 0
    global U = 0
    global U2 = 0
    global Xh = 0
    global r  = 0

    nothing
end

# return variable from the data #
macro load(var, inds...)
    len = length(inds)
    name = string(var)

    if len == 0
        q = quote
            if $var == 0
                global $var = read(rbfile, $name)
            end

            $var
        end
    else
        #temp code until getindex in JLD module is fixed for arrays of inds
        i1 = inds[1]

        if len == 1
            q = quote
                if $var == 0
                    global $var = read(rbfile, $name)
                end

                $var[$i1]
            end
        else
            i2 = inds[2]

            q = quote
                if $var == 0
                    global $var = read(rbfile, $name)
                end

                $var[$i1, $i2]
            end
        end
    end
end

# load variable readable from REPL #
function load(var::Symbol)
    name = string(var)

    q = quote
        if $var == 0
            global $var = read(rbfile, $name)
            export $var
        end

        $var
    end

    eval(q)
end

# pod related functions #

export poderror, pod, podsub, offsol, offtime

# compute pod for the selected parameters #
ParamsUnion = Union{Matrix{Vector{Float64}}, FloatRange{Float64}, Vector{Vector{Float64}}}
function pod(params::ParamsUnion, error::Float64 = 1.0e-6)
    offsol(params)

    pod(error)
end

function pod(params::Union{Matrix{Vector{Float64}}, FloatRange{Float64}, Vector{Vector{Float64}}}, N::Int64)
    offsol(params)

    pod(N)
end

function pod(error = 1.0e-6, S::Matrix{Float64} = @load S)
    pod(0, S)

    o_delete(rbfile, "N")
    rbfile["N"] = poderror(error)

    o_delete(rbfile, "N2")
    rbfile["N2"] = poderror(error, @load sigma2)

    nothing
end

function pod(N1::Int64, S::Matrix{Float64} = @load S)
    @load r
    @load ns
    @load Nh

    temp = fem.Xh()
    Xh = spzeros(Nh, Nh)

    sparsepytojl(Xh, temp)

    rbfile["Xh"] = Xh

    U = eigfact(Symmetric(S' * Xh * S), ns-r+1:ns)

    sigma = sqrt(abs(U[:values][end:-1:1]))
    U = S * U[:vectors][:, end:-1:1] * spdiagm(1 ./ sigma)

    rbfile["U"] = U
    rbfile["sigma"] = sigma
    rbfile["N"] = N1

    F = svdfact(S)

    rbfile["U2"] = F[:U][:, 1:r]
    rbfile["sigma2"] = F[:S][1:r]
    rbfile["N2"] = N1

    affinedim = fem.affinedim()

    if affinedim[1] == -1
        temp = fem.Afem()
        Afem = spzeros(Nh, Nh)

        sparsepytojl(Afem, temp)

        rbfile["Afem"] = Afem
    end

    if affinedim[2] == -1
        rbfile["bfem"] = fem.bfem()
    end

    Ut = U'

    for i = 1:affinedim[1]
        temp = fem.Afem(i)
        Afem = spzeros(Nh, Nh)

        sparsepytojl(Afem, temp)

        rbfile["Afem$i"] = (Ut * Afem * U)'
    end

    for i = 1:affinedim[2]
        rbfile["bfem$i"] = Ut * fem.bfem(i)
    end

    nothing
end

# compute the snapshots and prepare the jld file with the information #
function offsol(params1::Vector{Vector{Float64}})
    if exists(rbfile, "params") # check if there exists snapshots in params1 already computed in the jld
        @load S
        @load params
        @load ns

        ind = find(x -> length(find(y -> y == x , params)) != 0, params1) # indexes of the snapshots already in the jld from params1
        deleteat!(params1, ind) # delete those parameters in ind

        # compute and merge the new snapshots and the related info
        S1 = offcoef(params1)
        ns1 = length(params1)

        global S = [S S1]
        global ns += ns1
        global params = [params; params1]

        # delete the old info in the jld
        o_delete(rbfile, "S")
        o_delete(rbfile, "ns")
        o_delete(rbfile, "params")
        o_delete(rbfile, "r")
    else
        global S = offcoef(params1)
        global ns = length(params1)
        global params = params1
        global Nh = fem.Nh()

        rbfile["Nh"] = Nh
    end

    rbfile["S"] = S
    rbfile["ns"] = ns
    rbfile["params"] = params
    rbfile["r"] = rank(S)

    nothing
end

# this offsol function change the structure of params to be compatible with the function above #
offsol(params::Union{Matrix{Vector{Float64}}, FloatRange{Float64}}) = offsol(paramsvec(params))
# vectorize parameters #
paramsvec(params::Matrix{Vector{Float64}}) = vec(params)
paramsvec(params::Union{FloatRange{Float64}, Float64}) = map(x -> [x], params)

# call the fem solver for a parameter and returns its solution as a function #
function offsol(param::Vector{Float64})
    fem.offsol(param)

    x -> fem.offeval(x)
end

# call the fem solver for the parameters, returns the coefficients of the
# solutions in each column #
function offcoef(params::Vector{Vector{Float64}})
    Nh = fem.Nh()
    ns = length(params)

    sols = Matrix{Float64}(Nh, ns)

    for (ind, p) in enumerate(params)
        sols[:, ind] = offcoef(p)
    end

    sols
end

# call the fem solver for a singular param #
offcoef(param::Vector{Float64}) = fem.offcoef(param)

# compute for which N the truncation of pod gives the error bound indicated #
function poderror(error::Float64 = 1.0e-6, sigma::Vector{Float64} = @load sigma)
    sigma2 = sigma.^2
    s1 = sum(sigma2[2:end])

    r = @load r
    N = 1

    while s1 > error && N < r
        N += 1
        s1 -= sigma2[N]
    end

    N
end

# compute pod error for N elements #
function poderror(N::Int64, sigma::Vector{Float64} = @load sigma)
	  sigma2 = sigma.^2

    sum(sigma2[N+1:end])
end

# auxiliar function to convert scipy sparse matrix to julia sparse structure #
function sparsepytojl(A::SparseMatrixCSC{Float64,Int64}, o::PyObject)
    ind = o[:nonzero]()
    indx = ind[1] + 1
    indy = ind[2] + 1
    odata = filter(x -> x != 0, o[:data])

    for i in 1:length(indx)
        A[indx[i], indy[i]] = odata[i]
    end
end

# compute pod for a subset of parameters #
function podsub(inds::Vector{Int64}, name::AbstractString)
    params1 = @load(params, inds)
    S1 = @load(S, :, inds)
    @load(Afem)

    jldnew(name)

    pod(params1, S1, Afem)
end

# test how much time will take to finish the offline stage #
offtime(params::Union{Matrix{Vector{Float64}}, FloatRange{Float64}}) = offtime(paramsvec(params))
function offtime(params::Vector{Vector{Float64}})
	  offcoef(rand(params))

	  tic()
	  offcoef(rand(params))
	  time = toc();

    ns = length(params)
	  time = time * ns

    timeprint(time)
end

function timeprint(time::Float64)
    if time <= 60
		    print(time)
		    print(" sec")
	  elseif 60 < time && time <= 3600
		    print(time/60)
		    print(" min")
	  else
		    print(time/3600)
		    print(" hr")
	  end
end

#################################################################

# errors related functions #

export errors, poderrors, singularvalues, offontime

# compare the offline and online times #
offontime(param::Float64) = offontime(paramsvec(param))
function offontime(param::Vector{Float64})
    println("Offline computations:")
    tic()
    offcoef(param)
    off = toc()
    println("")

    println("Online computations:")
    tic()
    oncoef(param)
    on = toc()
    println("")

    diff = off - on

    if off >= on
        speedup = round(Int, 100 * diff / on)
        print("There is a speed up of $speedup%")
    else
        print("Online computations are slower than the offline ones with a difference of $diff sec")
    end
end

# compute the l2 errors from 1 to N for a parameter #
errors(param::Float64, N::Int64) = errors(paramsvec(param), N)
function errors(param::Vector{Float64}, N::Int64, UXh::Bool = true, normXh::Bool = true)
    if normXh
        Xh = @load Xh
        normfun = x -> sqrt(dot(x, Xh*x))
    else
        normfun = norm
    end

    off = offcoef(param)

    Afem = @load Afem
    bfem = fem.bfem(param)
    U = UXh ? rbfile["U"][:, 1:N] : rbfile["U2"][:, 1:N]

    errors = Float64[]

    for i = 1:N
        RB = U[:, 1:i]
        RBt = RB'

        on = RB * (((RBt * Afem * RB)') \ (RBt * bfem))

        push!(errors, normfun(off - on))
    end

    errors
end

# compute the pod errors from 1 to N #
poderrors(N::Int64, sigma::Vector{Float64} = @load sigma) = map(x -> poderror(x, sigma), 1:N)

# return the first N singular values
singularvalues(N::Int64, sigma::Vector{Float64} = @load sigma) = sigma[1:N]

# online stage related functions

export onbasis, onload, oncoef, onsol, ontofem

# return i-th base of the reduced basis #
function onbasis(i::Int64)
    coef = vec(rbfile["U"][:,i])

    fem.offsolcoef(coef)

    x -> fem.offeval(x)
end

.*(a::Vector{Vector{Float64}}, b::Vector{Float64}) = map(i -> a[i] * b[i], 1:length(b))
.*(a::Vector{SparseMatrixCSC{Float64,Int64}}, b::Vector{Float64}) = map(i -> a[i] * b[i], 1:length(b))

function prepare()
   @load Afem
   nothing
end

# compute the online coefficients for param #

@generated function oncoef(param::Vector{Float64})
    affinedim = fem.affinedim()

    RB = @load U
    RBt = RB'

    if affinedim[1] == -1 # Afem independent
        Apde = @load Afem
        ARB = factorize((RBt * Apde * RB)')
        if affinedim[2] == 0 # bfem not affine
            return :($ARB \ ($RBt * fem.bfem(param)))
        else # bfem affine
            bfems = Vector{Float64}[]

            for i = 1:affinedim[2]
                push!(bfems, read(rbfile, "bfem$i"))
            end

            return :($ARB \ sum($bfems .* fem.bcoefs(param)))
        end
    elseif affinedim[1] == 0 # Afem not affine
        if affinedim[2] == -1 # bfem independent
            bRB = RBt * (@load bfem)

            return :(($RBt * fem.Afem(param) * $RB)' \ $bRB)
        elseif affinedim[2] == 0 # bfem not affine
            return :(($RBt * fem.Afem(param) * $RB)' \ ($RBt * fem.bfem(param)))
        else # bfem affine
            bfems = Vector{Float64}[]

            for i = 1:affinedim[2]
                push!(bfems, read(rbfile, "bfem$i"))
            end

            return :(($RBt * fem.Afem(param) * $RB)' \ sum($bfems .* fem.bcoefs(param)))
        end
    else # Afem affine
        Afems = Matrix{Float64}[] # change this

        for i = 1:affinedim[1]
            push!(Afems, read(rbfile, "Afem$i"))
        end

        if affined[2] == -1 # bfem independent
            bRB = RBt * (@load bfem)

            return :(sum($Afems .* fem.Acoefs(param)) \ $bRB)
        elseif affined[2] == 0 # bfem not affine
            return :(sum($Afems .* fem.Acoefs(param)) \ ($RBt * fem.bfem(param)))
        else
            bfems = Vector{Float64}[]

            for i = 1:affinedim[2]
                push!(bfems, read(rbfile, "bfem$i"))
            end

            return :(sum($Afems .* fem.Acoefs(param)) \ sum($bfems .* fem.bcoefs(param)))
        end
    end
end

# return the online solution as a function #
function onsol(param::Vector{Float64})
    fem.offsolcoef(ontofem(oncoef(param)))

    x -> fem.offeval(x)
end

# convert the reduced basis coefficients to the fem coefficients #
@generated function ontofem(param::Vector{Float64})
    RB = @load U

    return :(RB * oncoef(param))
end

#=

# Solve the equation for a param and put in evals
# the full evaluations of some prevaluated values
# of the basis functions

export offsol!, onsol!, offbase, onbase, solve!, oncoef, test

test() = @load Afem

offsol!(param::Vector{Float64}, evals::Vector{Float64}, values::Vector{Vector{Float64}}) = solve!(param, evals, values, offcoef)

function onsol!(param::Vector{Float64}, evals::Vector{Float64}, values::Vector{Vector{Float64}})
    #@load Afem
    #@load U

    solve!(param, evals, values, oncoef)
end

function solve!(param::Vector{Float64}, evals::Vector{Float64}, values::Vector{Vector{Float64}}, coefun::Function)
    coef = coefun(param)

    for (ind, vals) in enumerate(values)
        evals[ind] = dot(coef, vals)
    end

    nothing
end

# Evaluate the base in x

@generated function offbase(x::Vector{Float64})
    fem.buildbase()

    :(fem.offbase(x))
end

@generated function onbase(x::Vector{Float64})
    if !fem.isRB()
        RBt = (@load U)'
        fem.buildRB(RBt)
    end

    :(fem.onbase(x))
end

=#

####################################################################################

# EIM functions

export offeim, oneim, eim, eimevals, eimevalsgrad

function eimevals(g::Function, paramsx::Vector{Vector{Float64}}, paramsmu::Vector{Vector{Float64}})
    evals = Vector{Float64}[]
    gs_x = map(mu -> (x -> g(x, mu)), paramsmu)

    for g_x in gs_x
        push!(evals, map(g_x, paramsx))
    end

    hcat(evals...)
end

function eimevalsgrad(g::Function, paramsx::Vector{Vector{Float64}}, paramsmu::Vector{Vector{Float64}})
    results = Array{DiffBase.DiffResult{1, Float64, Tuple{Array{Float64, 1}}}, 1}[]
    gs_x = map(mu -> (x -> g(x, mu)), paramsmu)

    result = DiffBase.GradientResult(paramsx[1])

    for g_x in gs_x
        push!(results, map(x -> ForwardDiff.gradient!(copy(result), g_x, x), paramsx))
    end

    evals = hcat(map(x -> map(DiffBase.value, x), results)...)
    grad1 = hcat(map(x -> map(y -> DiffBase.gradient(y)[1], x), results)...)
    grad2 = hcat(map(x -> map(y -> DiffBase.gradient(y)[2], x), results)...)

    (evals, grad1, grad2)
end

function offeim(evals::Matrix{Float64}, Q::Matrix{Float64}, J::Vector{Int64}, Jmu::Vector{Int64}, iter_max::Int64 = 2^11, tol::Float64 = 200eps())
    M = size(Q)[2]
    error_M = tol + 1

    siz = size(evals)

    Q = hcat(Q, Matrix{Float64}(siz[1], iter_max - M))

    minmaxeval = [0., 0.]
    function findmax(values::Matrix{Float64})
        imax = indmax(values)
        imin = indmin(values)
        minmaxeval[1] = values[imax]
        minmaxeval[2] = abs(values[imin])

        imax = (indmax(minmaxeval) == 1) ? imax : imin

        ind2sub(siz, imax)
    end

    Q_M = Q[:, 1:M]
    allres = evals - Q_M * (Q_M[J, :] \ evals[J, :])

    ind_mu_M = findmax(allres)

    error_M = abs(allres[ind_mu_M[1], ind_mu_M[2]])

    r = allres[:, ind_mu_M[2]]

    count = 0

    while M < iter_max && count < 8
        tic()
        M += 1

        #i_M = indmax(abs(r))
        i_M = ind_mu_M[1]
        rho_M = r / r[i_M]

        Q[:, M] = rho_M
        push!(J, i_M)
        push!(Jmu, ind_mu_M[2])

        Q_M = Q[:, 1:M]

        allres = evals - Q_M * (Q_M[J, :] \ evals[J, :])
        #=
        temp = Q_M[J, :] \ evals[J, :]
        maxvalue = 0
        for i in 1:siz[2]
            allres = evals[:, i] - Q_M * temp[:, i]
        end
        =#

        ind_mu_M = findmax(allres)

        error_M = abs(allres[ind_mu_M[1], ind_mu_M[2]])

        r = allres[:, ind_mu_M[2]]

        if error_M <= tol
            count += 1
        end

        toc()
        println(error_M)
    end

    (Q[:,1:M], J, Jmu, error_M)
end
function offeim(evals::Matrix{Float64}, iter_max::Int64 = 1500, tol::Float64 = 1e-11)
    M = 0
    error_M = tol + 1

    siz = size(evals)
    iter_max = min(iter_max, siz[2])

    Q = Matrix{Float64}(siz[1], iter_max)
    J = Int64[]
    Jmu = Int64[]

    minmaxeval = [0., 0.]
    function findmax(values::Matrix{Float64})
        imax = indmax(values)
        imin = indmin(values)
        minmaxeval[1] = values[imax]
        minmaxeval[2] = abs(values[imin])

        imax = (indmax(minmaxeval) == 1) ? imax : imin

        ind2sub(siz, imax)
    end

    ind_mu_M = findmax(evals)

    r = evals[:, ind_mu_M[2]]

    count = 0

    while M < iter_max && count < 1
        tic()
        M += 1

        #i_M = indmax(abs(r))
        i_M = ind_mu_M[1]
        rho_M = r / r[i_M]

        Q[:, M] = rho_M
        push!(J, i_M)
        push!(Jmu, ind_mu_M[2])

        Q_M = Q[:, 1:M]

        allres = evals - Q_M * (Q_M[J, :] \ evals[J, :])
        #=
        temp = Q_M[J, :] \ evals[J, :]
        maxvalue = 0
        for i in 1:siz[2]
            allres = evals[:, i] - Q_M * temp[:, i]
        end
        =#

        ind_mu_M = findmax(allres)

        error_M = abs(allres[ind_mu_M[1], ind_mu_M[2]])

        r = allres[:, ind_mu_M[2]]

        if error_M <= tol
            count += 1
        end

        toc()
        println(error_M)
    end

    (Q[:,1:M], J, Jmu, error_M)
end

function offeim(evals::Matrix{Float64}, grad1::Matrix{Float64}, grad2::Matrix{Float64}, iter_max::Int64 = 2^11, tol::Float64 = 200eps())
    M = 0
    error_M = tol + 1
    error_Mgrad1 = 0
    error_Mgrad2 = 0

    siz = size(evals)

    Q = Matrix{Float64}(siz[1], iter_max)
    Qgrad1 = Matrix{Float64}(siz[1], iter_max)
    Qgrad2 = Matrix{Float64}(siz[1], iter_max)
    J = Int64[]

    minmaxeval = [0., 0.]
    function findmax(values::Matrix{Float64})
        imax = indmax(values)
        imin = indmin(values)
        minmaxeval[1] = values[imax]
        minmaxeval[2] = abs(values[imin])

        imax = (indmax(minmaxeval) == 1) ? imax : imin

        ind2sub(siz, imax)
    end

    ind_mu_M = findmax(evals)

    r = evals[:, ind_mu_M[2]]
    rgrad1 = grad1[:, ind_mu_M[2]]
    rgrad2 = grad2[:, ind_mu_M[2]]

    rref = r

    count = 0

    while M < iter_max && count < 8
        tic()
        M += 1

        i_M = indmax(abs(rref))
        rho_M = r / rref[i_M]
        rho_Mgrad1 = rgrad1 / rref[i_M]
        rho_Mgrad2 = rgrad2 / rref[i_M]

        Q[:, M] = rho_M
        Qgrad1[:, M] = rho_Mgrad1
        Qgrad2[:, M] = rho_Mgrad2
        push!(J, i_M)

        Q_M = Q[:, 1:M]
        Q_Mgrad1 = Qgrad1[:, 1:M]
        Q_Mgrad2 = Qgrad2[:, 1:M]

        coefs = Q_M[J, :] \ evals[J, :]
        allres = evals - Q_M * coefs
        allresgrad1 = grad1 - Q_Mgrad1 * coefs
        allresgrad2 = grad2 - Q_Mgrad2 * coefs
        #=
        temp = Q_M[J, :] \ evals[J, :]
        maxvalue = 0
        for i in 1:siz[2]
            allres = evals[:, i] - Q_M * temp[:, i]
        end
        =#

        ind_mu_M = findmax(allres)
        ind_mu_Mgrad1 = findmax(allresgrad1)
        ind_mu_Mgrad2 = findmax(allresgrad2)

        error_M = abs(allres[ind_mu_M[1], ind_mu_M[2]])
        error_Mgrad1 = abs(allresgrad1[ind_mu_Mgrad1[1], ind_mu_Mgrad1[2]])
        error_Mgrad2 = abs(allresgrad2[ind_mu_Mgrad2[1], ind_mu_Mgrad2[2]])

        max123 = indmax([error_M, error_Mgrad1, error_Mgrad2])
        indref = ind_mu_M[2]

        if max123 == 1
            rref = r
        elseif max123 == 2
            indref = ind_mu_Mgrad1[2]
            rref = rgrad1
        else
            indref = ind_mu_Mgrad2[2]
            rref = rgrad2
        end

        r = allres[:, indref]
        rgrad1 = grad1[:, indref]
        rgrad2 = grad2[:, indref]

        if max123 == 1
            rref = r
        elseif max123 == 2
            indref = ind_mu_Mgrad1[2]
            rref = rgrad1
        else
            indref = ind_mu_Mgrad2[2]
            rref = rgrad2
        end

        if error_M <= tol && error_Mgrad1 <= tol && error_Mgrad2 <= tol
            count += 1
        end

        toc()
        println(error_M)
        println(error_Mgrad1)
        println(error_Mgrad2)
        println("ok")
    end

    (Q[:,1:M], Qgrad1[:,1:M], Qgrad2[:,1:M], J, error_M, error_Mgrad1, error_Mgrad2)
end

function offeim(g::Function, paramsx::Vector{Vector{Float64}}, paramsmu::Vector{Vector{Float64}}, iter_max::Int64 = 2^14, tol::Float64 = 200eps(), preevaluate::Bool = true)
    M = 0
    error_M = tol + 1

    siz = (length(paramsx), length(paramsmu))

    Q = Matrix{Float64}(siz[1], iter_max)
    J = Int64[]
    mus = Vector{Float64}[]

    if preevaluate
        evals = Vector{Float64}[]
        gs_x = map(mu -> (x -> g(x, mu)), paramsmu)

        for g_x in gs_x
            push!(evals, map(g_x, paramsx))
        end

        evals = hcat(evals...)

        ind_mu_M = ind2sub(siz, indmax(abs(evals)))

        r = evals[:, ind_mu_M[2]]

        count = 0

        while M < iter_max && count < 8
            tic()
            M += 1

            i_M = indmax(abs(r))
            rho_M = r / r[i_M]

            Q[:, M] = rho_M
            push!(J, i_M)

            Q_M = Q[:, 1:M]

            allres = evals - Q_M * (Q_M[J, :] \ evals[J, :])

            ind_mu_M = ind2sub(siz, indmax(abs(allres)))

            error_M = abs(allres[ind_mu_M[1], ind_mu_M[2]])

            r = allres[:, ind_mu_M[2]]

            if error_M <= tol
                count += 1
            end

            toc()
            println(error_M)
        end
    else
        mu_M = argmax(g, paramsx, paramsmu)[2]

        r = map(x -> g(x, mu_M), paramsx)

        count = 0

        while M < iter_max && count < 8
            tic()
            M += 1

            i_M = indmax(abs(r))
            rho_M = r / r[i_M]

            Q[:, M] = rho_M
            push!(J, i_M)

            Q_M = Q[:, 1:M]

            (error_M, mu_M) = argmax(g, Q_M, J, paramsx, paramsmu)

            gmu = map(x -> g(x, mu_M), paramsx)

            r = gmu - Q_M*(Q_M[J, :] \ gmu[J])

            if error_M <= tol
                count += 1
            end

            toc()
            println(error_M)
        end
    end

    return (Q[:,1:M], mus, J, error_M)
end

function argmax(g::Function, paramsx::Vector{Vector{Float64}}, paramsmu::Vector{Vector{Float64}})
    fx = map(mu -> (x -> g(x, mu)), paramsmu)

    maxval = 0
    maxmu = 0

    for (i, f) in enumerate(fx)
        m = maximum(abs(map(f, paramsx)))

        if m > maxval
            maxval = m
            maxmu = paramsmu[i]
        end
    end

    return (maxval, maxmu)
end

function argmax(g::Function, Q::Matrix{Float64}, J::Vector{Int64}, paramsx::Vector{Vector{Float64}}, paramsmu::Vector{Vector{Float64}})
    fx = map(mu -> (x -> g(x, mu)), paramsmu)

    len = length(J)

    maxval = 0
    maxmu = 0

    Q1 = factorize(Q[J, :])

    for (i, f) in enumerate(fx)
        gmu = map(f, paramsx)
        m = maximum(abs(gmu - Q*(Q1 \ gmu[J])))

        if m > maxval
            maxval = m
            maxmu = paramsmu[i]
        end
    end

    return (maxval, maxmu)
end

function oneim(mu::Vector{Float64}, g::Function, Q::Matrix{Float64}, J::Vector{Int64}, paramsx::Vector{Vector{Float64}})
    gJmu = map(x -> g(x, mu), paramsx[J])

    Q[J, :] \ gJmu
end

oneim(mu::Vector{Float64}, g::Function, QQ_J::Base.LinAlg.LU{Float64,Array{Float64,2}}, paramsx_J::Vector{Vector{Float64}}) = QQ_J \ map(x -> g(x, mu), paramsx_J)

function oneim(mu1, mu2, g::Function, QQ_J::Base.LinAlg.LU{Float64,Array{Float64,2}}, paramsx_J::Vector{Vector{Float64}})
    mu = [mu1, mu2]
    QQ_J \ map(x -> g(x, mu), paramsx_J)
end

function eim(g::Function, mus::Vector{Vector{Float64}}, xs::Vector{Vector{Float64}})
    basis = Function[]
    coefs = [0.0]

    push!(basis, x -> g(x, mus[1]) / g(xs[1], mus[1]))

    function Ig_mu(x::Vector{Float64}, M::Int64)
        result = 0

        for (j, rho) in enumerate(basis[1:M])
            result += coefs[j] * rho(x)
        end

        result
    end

    for (i, mu) in enumerate(mus[2:end])
        tic()
        B_M = Matrix{Float64}(i, i)
        g_M = map(x -> g(x, mu), xs[1:i])

        for (j, rho) in enumerate(basis)
            B_M[:, j] = map(rho, xs[1:i])
        end

        coefs = B_M \ g_M

        push!(basis, x -> (g(x, mu) - Ig_mu(x, i)) / (g(xs[i+1], mu) - Ig_mu(xs[i+1], i)))
        toc()
    end

    basis
end

#######################################################################################

# Unfinished functions

#=
function greedy(params::Vector{Vector{Float64}}, mu1::Vector{Float64}, errormax::Float64 = 1e-6, Nmax::Int64 = 2000)
    Nh = fem.Nh()

    temp = fem.Afem()
    Afem = spzeros(Nh, Nh)
    sparsepytojl(Afem, temp)
    temp = 0

    c = 1 / 1.e-5

    len = length(params)
    V = Matrix{Float64}(fem.Nh(), Nmax)
    bfems = Matrix{Float64}(Nh, len)
    sampling = Vector{Float64}[]

    tic()
    for i in 1:len
        bfems[:, i] = fem.bfem(params[i])
    end
    toc()

    muN = mu1
    N, error = 0, errormax + 1

    function gramschmidt(sol::Vector{Float64}, N::Int64)
        if N != 0
            sol = sol - V[:, 1:N] * (V[:, 1:N])' * sol
        end

        sol = sol / norm(sol)
        nothing
    end
    function argmax(N::Int64)
        err = 0
        mu = 0

        VN = V[:, 1:N]
        VNt = VN'
        onsols = (VNt * Afem * VN)' \ (VNt * bfems)
        res = c * (bfems - Afem * VN * onsols)

        for i in 1:size(res)[2]
            normres = norm(res[:, i])

            if normres >= err
                err = normres
                mu = param
            end
        end

        (err, mu)
    end
    while N < Nmax && error > errormax
        tic()
        N += 1

        sol = offcoef(muN)
        gramschmidt(sol, N-1)

        V[:, N] = sol
        push!(sampling, muN)

        error, muN = argmax(N)
        println(error)
        toc()
    end
end

=#

end


text:
# specific stuff for Oct trees

mutable struct OctTreeNode{T<:AbstractPoint3D} <: SpatialTreeNode
    id::Int64
    r::Float64
    midx::Float64
    midy::Float64
    midz::Float64
    is_empty::Bool
    is_divided::Bool
    point::T
    lxlylz::OctTreeNode{T}
    lxhylz::OctTreeNode{T}
    hxlylz::OctTreeNode{T}
    hxhylz::OctTreeNode{T}
    lxlyhz::OctTreeNode{T}
    lxhyhz::OctTreeNode{T}
    hxlyhz::OctTreeNode{T}
    hxhyhz::OctTreeNode{T}
    function OctTreeNode{T}(r::Number, midx::Number, midy::Number, midz::Number) where T <:AbstractPoint3D
        n = new(0, r, midx, midy, midz, true, false, T())
        n.lxlylz = n
        n.lxhylz = n
        n.hxlylz = n
        n.hxhylz = n
        n.lxlyhz = n
        n.lxhyhz = n
        n.hxlyhz = n
        n.hxhyhz = n
        n
    end
end

OctTreeNode(r::Number, midx::Number, midy::Number, midz::Number, ::Type{T}) where T <:AbstractPoint3D = OctTreeNode{T}(r, midx, midy, midz)
OctTreeNode(::Type{T}) where T <:AbstractPoint3D = OctTreeNode(0.5, 0.5, 0.5, 0.5, T)
OctTreeNode() = OctTreeNode(Point3D)

mutable struct OctTree{T<:AbstractPoint3D} <: SpatialTree
    head::OctTreeNode{T}
    number_of_nodes_used::Int64
    nodes::Array{OctTreeNode, 1}
    faststack::Array{OctTreeNode{T}, 1}
    function OctTree{T}(r::Number, midx::Number, midy::Number, midz::Number, n::Int64=100000) where T <:AbstractPoint3D
        nodes = OctTreeNode[OctTreeNode(r, midx, midy, midz, T) for i in 1:n]
        new(nodes[1], 1, nodes, [OctTreeNode(T) for i in 1:10000])
    end
end

OctTree(r::Number, midx::Number, midy::Number, midz::Number, ::Type{T};n=1000) where T <:AbstractPoint3D = OctTree{T}(r, midx, midy, midz, n)
OctTree(::Type{T};n=10000) where T <: AbstractPoint3D = OctTree(0.5, 0.5, 0.5, 0.5, T; n=n)
OctTree(n::Int64) = OctTree(Point3D; n=n)
OctTree() = OctTree(Point3D)
Base.eltype(::OctTree{T}) where T <:AbstractPoint3D = T

function initnode!(q::OctTreeNode{T}, r::Number, midx::Number, midy::Number, midz::Number) where T <:AbstractPoint3D
    q.r = r
    q.midx = midx
    q.midy = midy
    q.midz = midz
    q.is_empty = true
    q.is_divided = false
    q.id = 0
end

function divide!(h::OctTree{T}, q::OctTreeNode{T}) where T<:AbstractPoint3D
    # make sure we have enough nodes
    if length(h.nodes) - h.number_of_nodes_used < 8
        new_size = length(h.nodes)+(length(h.nodes) >>> 1)
        sizehint!(h.nodes, new_size)
        for i in 1:(new_size-length(h.nodes))
            push!(h.nodes, OctTreeNode(T))
        end
    end

    # populate new nodes
    q.lxlylz = h.nodes[h.number_of_nodes_used+1]
    q.lxhylz = h.nodes[h.number_of_nodes_used+2]
    q.hxlylz = h.nodes[h.number_of_nodes_used+3]
    q.hxhylz = h.nodes[h.number_of_nodes_used+4]
    q.lxlyhz = h.nodes[h.number_of_nodes_used+5]
    q.lxhyhz = h.nodes[h.number_of_nodes_used+6]
    q.hxlyhz = h.nodes[h.number_of_nodes_used+7]
    q.hxhyhz = h.nodes[h.number_of_nodes_used+8]

    # set new nodes properties (dimensions etc.)
    r2 = q.r/2.0
    initnode!(q.lxlylz, r2, q.midx-r2, q.midy-r2, q.midz-r2)
    initnode!(q.lxhylz, r2, q.midx-r2, q.midy+r2, q.midz-r2)
    initnode!(q.hxlylz, r2, q.midx+r2, q.midy-r2, q.midz-r2)
    initnode!(q.hxhylz, r2, q.midx+r2, q.midy+r2, q.midz-r2)
    initnode!(q.lxlyhz, r2, q.midx-r2, q.midy-r2, q.midz+r2)
    initnode!(q.lxhyhz, r2, q.midx-r2, q.midy+r2, q.midz+r2)
    initnode!(q.hxlyhz, r2, q.midx+r2, q.midy-r2, q.midz+r2)
    initnode!(q.hxhyhz, r2, q.midx+r2, q.midy+r2, q.midz+r2)

    # update tree and parent node
    h.number_of_nodes_used += 8
    q.is_divided = true
    if !q.is_empty
        # move point in parent node to child node
        sq = _getsubnode(q, q.point)
        sq.is_empty = false
        q.is_empty = true
        sq.point = q.point
        q.point = T()
    end
    q
end

function _getsubnode(q::OctTreeNode{T}, point::T) where T <:AbstractPoint3D
    x=getx(point)
    y=gety(point)
    z=getz(point)
    if x<q.midx
        # lx
        if y<q.midy
            # ly
            z<q.midz && return q.lxlylz
            return q.lxlyhz
        else
            # hy
            z<q.midz && return q.lxhylz
            return q.lxhyhz
        end
    else
        # hx
        if y<q.midy
            # ly
            z<q.midz && return q.hxlylz
            return q.hxlyhz
        else
            # hy
            z<q.midz && return q.hxhylz
            return q.hxhyhz
        end
    end
end

function map(t::OctTree{T}, cond_data) where T <: AbstractPoint3D
    curr_stack_ix = 1
    t.faststack[1] = t.head
    @inbounds while curr_stack_ix > 0
        q = t.faststack[curr_stack_ix]
        curr_stack_ix -= 1
        if !stop_cond(q, cond_data) && q.is_divided
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxlylz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxlyhz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxhylz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.lxhyhz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.hxlylz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.hxlyhz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.hxhylz
            curr_stack_ix += 1
            t.faststack[curr_stack_ix] = q.hxhyhz
        end
    end
end


text:
using Mocha

backend = GPUBackend()
init(backend)

include("mnist-layers.jl")

train_net = Net("train-net", backend, [train_dl, non_data_layers...])
test_net = Net("test-net", backend, [test_dl, non_data_layers...])

adam_instance = Adam()
params = make_solver_parameters(adam_instance;
    max_iter=500_000,
    regu_coef=1e-4,
    lr_policy=LRPolicy.Fixed(1e-3))
solver = Solver(adam_instance, params)
add_coffee_break(solver, TrainingSummary(), every_n_iter=10000)
add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=10000)

base_dir = "snapshots_tues"
setup_coffee_lounge(solver, save_into="$base_dir/statistics.jld", every_n_iter=10000)
add_coffee_break(solver, Snapshot(base_dir), every_n_iter=10000)

solve(solver, train_net)

destroy(train_net)
destroy(test_net)

shutdown(backend)

1


text:
const cc = CrawlCorpus(docs_location, true)

_next_doc_id = 1
function allocate_docids(n)
    global _next_doc_id
    docids = _next_doc_id:(_next_doc_id+n-1)
    _next_doc_id += n
    docids
end

function assign_unique_docids(crps::Corpus)
    ndocs = length(crps.documents)

    # allocate doc ids from one master processor
    docids = remotecall_fetch(1, allocate_docids, ndocs)

    cdocs = documents(crps)

    # create a lookup map
    doc_to_id = Dict{String,Int}()
    id_to_doc = Dict{Int,String}()
    for docidx in 1:length(docids)
        d_id = docids[docidx]
        d_name = cdocs[docidx].metadata.name

        doc_to_id[d_name] = d_id
        id_to_doc[d_id] = d_name
    end

    invidx = inverse_index(crps)
    # update the doc ids in inverted index to global doc ids
    for (n,v) in invidx
        for docidx in 1:length(v)
            docid = v[docidx]
            doc = cdocs[docid]
            path = doc.metadata.name
            v[docidx] = doc_to_id[path]
        end
    end

    # return the id-to-doc map to be used during search
    id_to_doc
end

function store_index(invidx::Dict, id_to_doc::Dict{Int,String}, sername::String)
    path = joinpath(part_idx_location, sername)
    as_serialized(invidx, path)

    path = joinpath(id_to_doc_location, sername)
    as_serialized(id_to_doc, path)
end

function mkstr(ascdata::Array{Uint8,1})
    const l::Int = length(ascdata)
    for i in 1:l
        (ascdata[i] >= 0x80) && (ascdata[i] = ' ')
    end
    ASCIIString(ascdata)
    #UTF8String(ascdata)
end

function index_archive_file(archive::URI)
    fname = basename(archive.path)
    # since we may end up encountering archive files files from different segments in the crawl corpus with possibly same name,
    # create a unique id as hash of the complete URL.
    uniqid = bytes2hex(AWS.Crypto.md5(archive.path))
    sername = "$(uniqid).jser"

    # restartability aid: check if archive is indexed already and not repeat it
    if isfile(joinpath(part_idx_location, sername))
        println("skipping $archive already indexed as $(sername)")
        return fname
    end

    t1 = time()
    # read all entries from the archive
    f = open(cc, archive)
    entries = read_entries(cc, f, "text/")
    close(f)
    println("\t\tentry created at $(time()-t1)secs")

    # create a corpus of StringDocuments
    docs = GenericDocument[]
    for entry in entries
        # ignoring anything otherthan ascii for this example
        ascstr = mkstr(entry.data)
        sd = StringDocument(ascstr)
        sd.metadata = TextAnalysis.DocumentMetadata()
        sd.metadata.name = entry.uri
        push!(docs, sd)
    end
    crps = Corpus(docs)
    println("\t\tcorpus created at $(time()-t1)secs")

    # process the corpus
    crps = preprocess(crps)
    println("\t\tpreprocessed at $(time()-t1)secs")
    update_inverse_index!(crps)
    invidx = inverse_index(crps)

    # in the corpus, doc ids are allocated starting from 1
    # since we are processing in chunks, we need to make them unique across all runs and all processors
    docids = assign_unique_docids(crps)
    println("\tprocessing time $(time()-t1)secs")

    t1 = time()
    store_index(invidx, docids, sername)
    println("\tstoring time $(time()-t1)secs")

    # remove the archive file to save space. can be commented if disk space is not an issue
    clear_cache(cc, archive)
    fname
end

##
# execute index creation parallely for n_archives archive files
function create_index(n_archives::Int)
    arc_file_uris = archives(cc, n_archives)

    println("working with $(length(arc_file_uris)) archive files")
    arcs_indexed = @parallel append! for arc in arc_file_uris
        index_archive_file(arc)
        [string(arc)]
    end
    println("archives indexed: $(length(arcs_indexed))")
end


text:
module testMod;

export test;

function test(x::Int64)
   sum([i for i=1:x])
end

end

max_stars_repo_name: Alfheim/node-julia

max_issues_repo_name: Alfheim/node-julia

max_forks_repo_name: Alfheim/node-julia


text:
module Foo

using Hyperspecialize

import Bar

@concretize AlsoNotAType Set{Type}([Int64])
@concretize Float64 Set{Type}([UInt16])

# Make sure that other modules can also change concretizations
function __init__()
  @widen Bar.Float64 Set{Type}([UInt16, Float32])
end

end #module


text:
using BinaryBuilder

# Collection of sources required to build MbedTLS
# mbedtls release: 2.12.0
sources = [
    "https://github.com/ARMmbed/mbedtls.git" =>
    "294ad8d725c159400d10599cb26a4579a2d39956",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mbedtls/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DUSE_SHARED_MBEDTLS_LIBRARY=On
make -j${nproc} && make install
if [ $target == "x86_64-w64-mingw32" ]; then
    cp $prefix/lib/*.dll $prefix/bin/.
elif [ $target == "i686-w64-mingw32" ]; then
    cp $prefix/lib/*.dll $prefix/bin/.
else
    cd $prefix/lib; for f in $(find . -name '*.so'); do strip $f ; done
fi
"""

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

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libmbedcrypto", :libmbedcrypto),
    LibraryProduct(prefix, "libmbedtls", :libmbedtls),
    LibraryProduct(prefix, "libmbedx509", :libmbedx509)
]

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

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

max_stars_repo_name: invenia/MbedTLSBuilder

max_issues_repo_name: invenia/MbedTLSBuilder

max_forks_repo_name: invenia/MbedTLSBuilder


text:
using MatrixFactorizations, LinearAlgebra, Random, ArrayLayouts, Test
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!, require_one_based_indexing, checksquare
import MatrixFactorizations: QRCompactWYQLayout

struct MyMatrix <: LayoutMatrix{Float64}
    A::Matrix{Float64}
end

Base.getindex(A::MyMatrix, k::Int, j::Int) = A.A[k,j]
Base.setindex!(A::MyMatrix, v, k::Int, j::Int) = setindex!(A.A, v, k, j)
Base.size(A::MyMatrix) = size(A.A)
Base.strides(A::MyMatrix) = strides(A.A)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()

struct MyVector <: LayoutVector{Float64}
    A::Vector{Float64}
end

Base.getindex(A::MyVector, k::Int) = A.A[k]
Base.setindex!(A::MyVector, v, k::Int) = setindex!(A.A, v, k)
Base.size(A::MyVector) = size(A.A)
Base.strides(A::MyVector) = strides(A.A)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyVector) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyVector}) = DenseColumnMajor()

include("test_ul.jl")

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

Random.seed!(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

# helper functions to unambiguously recover explicit forms of an implicit QL Q
squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); lmul!(Q, Matrix{eltype(Q)}(I, sq, sq)))
rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)

@testset "QR" begin
    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        raw_a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        raw_a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
        asym = raw_a' + raw_a                  # symmetric indefinite
        apd  = raw_a' * raw_a                 # symmetric positive-definite
        ε = εa = eps(abs(float(one(eltya))))

        @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
            raw_b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
            εb = eps(abs(float(one(eltyb))))
            ε = max(εa, εb)
            tab = promote_type(eltya, eltyb)

            @testset "QR decomposition of a Number" begin
                α = rand(eltyb)
                aα = fill(α, 1, 1)
                @test qrunblocked(α).Q * qrunblocked(α).R ≈ qrunblocked(aα).Q * qrunblocked(aα).R
                @test abs(qrunblocked(α).Q[1,1]) ≈ one(eltyb)
            end

            for (a, b) in ((raw_a, raw_b),
                (view(raw_a, 1:n-1, 1:n-1), view(raw_b, 1:n-1, 1)))
                a_1 = size(a, 1)
                @testset "QR decomposition (without pivoting)" begin
                    qra   = @inferred qrunblocked(a)
                    @inferred qrunblocked(a)
                    q, r = qra.Q, qra.R
                    @test_throws ErrorException qra.Z
                    @test q'*squareQ(q) ≈ Matrix(I, a_1, a_1)
                    @test q*squareQ(q)' ≈ Matrix(I, a_1, a_1)
                    @test q'*Matrix(1.0I, a_1, a_1)' ≈ squareQ(q)'
                    @test squareQ(q)'q ≈ Matrix(I, a_1, a_1)
                    @test Matrix(1.0I, a_1, a_1)'q' ≈ squareQ(q)'
                    @test q*r ≈ a
                    @test a*(qra\b) ≈ b atol=3000ε
                    @test Array(qra) ≈ a
                    sq = size(q.factors, 2)
                    @test *(Matrix{eltyb}(I, sq, sq), adjoint(q)) * squareQ(q) ≈ Matrix(I, sq, sq) atol=5000ε
                    if eltya != Int
                        @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab}, q)
                        ac = copy(a)
                        @test qrunblocked!(a[:, 1:5])\b == qrunblocked!(view(ac, :, 1:5))\b
                    end
                    qrstring = sprint((t, s) -> show(t, "text/plain", s), qra)
                    rstring  = sprint((t, s) -> show(t, "text/plain", s), r)
                    qstring  = sprint((t, s) -> show(t, "text/plain", s), q)
                    @test qrstring == "$(summary(qra))\nQ factor:\n$qstring\nR factor:\n$rstring"
                end
                @testset "Thin QR decomposition (without pivoting)" begin
                    qra   = @inferred qrunblocked(a[:, 1:n1], Val(false))
                    @inferred qrunblocked(a[:, 1:n1], Val(false))
                    q,r   = qra.Q, qra.R
                    @test_throws ErrorException qra.Z
                    @test q'*squareQ(q) ≈ Matrix(I, a_1, a_1)
                    @test q'*rectangularQ(q) ≈ Matrix(I, a_1, n1)
                    @test q*r ≈ a[:, 1:n1]
                    @test q*b[1:n1] ≈ rectangularQ(q)*b[1:n1] atol=100ε
                    @test q*b ≈ squareQ(q)*b atol=100ε
                    @test_throws DimensionMismatch q*b[1:n1 + 1]
                    @test_throws DimensionMismatch b[1:n1 + 1]*q'
                    sq = size(q.factors, 1)
                    @test *(LowerTriangular(Matrix{eltyb}(I, sq, sq)), adjoint(q))*squareQ(q) ≈ Matrix(I, a_1, a_1) atol=5000ε
                    if eltya != Int
                        @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab},q)
                    end
                end
            end
            if eltya != Int
                @testset "Matmul with QR factorizations" begin
                    a = raw_a

                    qra = qrunblocked(a[:,1:n1], Val(false))
                    q, r = qra.Q, qra.R
                    @test rmul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
                    @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
                    @test rmul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
                    @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),adjoint(q))
                    @test_throws ErrorException size(q,-1)
                    @test_throws DimensionMismatch q * Matrix{Int8}(I, n+4, n+4)
                end
            end
        end
    end

    @testset "transpose errors" begin
        @test_throws MethodError transpose(qrunblocked(randn(3,3)))
        @test_throws MethodError adjoint(qrunblocked(randn(3,3)))
        @test_throws MethodError transpose(qrunblocked(randn(3,3), Val(false)))
        @test_throws MethodError adjoint(qrunblocked(randn(3,3), Val(false)))
        @test_throws MethodError transpose(qrunblocked(big.(randn(3,3))))
        @test_throws MethodError adjoint(qrunblocked(big.(randn(3,3))))
    end

    @testset "Issue 7304" begin
        A = [-√.5 -√.5; -√.5 √.5]
        Q = rectangularQ(qrunblocked(A).Q)
        @test norm(A-Q) < eps()
    end

    @testset "qrunblocked on AbstractVector" begin
        vl = [3.0, 4.0]
        for Tl in (Float32, Float64)
            for T in (Tl, Complex{Tl})
                v = convert(Vector{T}, vl)
                nv, nm = qrunblocked(v)
                @test nv*nm ≈ v
            end
        end
    end

    @testset "Issue 24589. Promotion of rational matrices" begin
        A = rand(1//1:5//5, 4,3)
        @test first(qrunblocked(A)) == first(qrunblocked(float(A)))
    end

    @testset "Issue Test Factorization fallbacks for rectangular problems" begin
        A  = randn(3,2)
        Ac = copy(A')
        b  = randn(3)
        b0 = copy(b)
        c  = randn(2)
        @test A \b ≈ ldiv!(c, qrunblocked(A ), b)
    end

    @testset "Wide QR" begin
        A = randn(3,5)
        Q,R = qrunblocked(A)
        @test Q*R ≈ A
    end

    @testset "lmul!/rmul!" begin
        A = randn(100,100)
        Q,R = qrunblocked(A)
        x = randn(100)
        b = randn(100,2)
        @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
        @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
        @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
        @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
        c = randn(2,100)
        @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
        @test rmul!(copy(c), Q') ≈ c*Matrix(Q')

        A = randn(100,103)
        Q,R = qrunblocked(A)
        x = randn(100)
        b = randn(100,2)
        @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
        @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
        @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
        @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
        c = randn(2,100)
        @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
        @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
    end

    @testset "LinearAlgebra.QR -> MatrixFactorizations.QR" begin
        A = randn(10,10)
        F̃ = LinearAlgebra.qrfactUnblocked!(copy(A))
        F = MatrixFactorizations.QR(F̃)
        Q,R = F
        @test Q*R ≈ A
        Q̃,_ = F̃
        @test MatrixFactorizations.QRPackedQ(Q̃)*R ≈ A
    end

    @testset "Ambiguity" begin
        A = randn(10,10)
        b = rand(10)
        B = randn(10,3)
        F = qrunblocked(A)
        @test ldiv!(F, MyVector(copy(b))) ≈ A \ b
        @test ldiv!(F, MyMatrix(copy(B))) ≈ A \ B
    end
end

@testset "QL" begin
    @testset "Compare with QR" begin
        n = 10
        A = randn(n,n)
        Q̃, R = @inferred qr(A[end:-1:1,end:-1:1])
        Q, L = @inferred ql(A)
        Q̄, L̄ = MatrixFactorizations.generic_qlfactUnblocked!(copy(A))
        @test Q.factors ≈ Q̄.factors
        @test Q.τ ≈ Q̄.τ
        @test R[end:-1:1,end:-1:1] ≈ L ≈ L̄
        @test Q̃[end:-1:1,end:-1:1] ≈ Q ≈ Q̄

        A = randn(n,n+2)
        Q̃, R = qr(A[end:-1:1,end:-1:1])
        Q, L = ql(A)
        Q̄, L̄ = MatrixFactorizations.generic_qlfactUnblocked!(copy(A))
        @test Q.factors ≈ Q̄.factors
        @test Q.τ ≈ Q̄.τ
        @test R[end:-1:1,end:-1:1] ≈ L ≈ L̄
        @test Q̃[end:-1:1,end:-1:1] ≈ Q ≈ Q̄

        A = randn(n+2,n)
        Q̃, R = qr(A[end:-1:1,end:-1:1])
        Q, L = ql(A)
        Q̄, L̄ = MatrixFactorizations.generic_qlfactUnblocked!(copy(A))
        @test Q.factors ≈ Q̄.factors
        @test Q.τ ≈ Q̄.τ
        @test R[end:-1:1,end:-1:1] ≈ L ≈ L̄
        @test Q̃[end:-1:1,end:-1:1] ≈ Q ≈ Q̄
    end

    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        raw_a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        raw_a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
        asym = raw_a' + raw_a                  # symmetric indefinite
        apd  = raw_a' * raw_a                 # symmetric positive-definite
        ε = εa = eps(abs(float(one(eltya))))

        @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
            raw_b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
            εb = eps(abs(float(one(eltyb))))
            ε = max(εa, εb)
            tab = promote_type(eltya, eltyb)

            @testset "QL decomposition of a Number" begin
                α = rand(eltyb)
                aα = fill(α, 1, 1)
                @test ql(α).Q * ql(α).L ≈ ql(aα).Q * ql(aα).L
                @test abs(ql(α).Q[1,1]) ≈ one(eltyb)
            end

            for (a, b) in ((raw_a, raw_b),
                (view(raw_a, 1:n-1, 1:n-1), view(raw_b, 1:n-1, 1)))
                a_1 = size(a, 1)
                @testset "QL decomposition (without pivoting)" begin
                    qla   = @inferred ql(a)
                    @inferred ql(a)
                    q, l  = qla.Q, qla.L
                    @test_throws ErrorException qla.Z
                    @test q'*squareQ(q) ≈ Matrix(I, a_1, a_1)
                    @test q*squareQ(q)' ≈ Matrix(I, a_1, a_1)
                    @test q'*Matrix(1.0I, a_1, a_1)' ≈ squareQ(q)'
                    @test squareQ(q)'q ≈ Matrix(I, a_1, a_1)
                    @test Matrix(1.0I, a_1, a_1)'q' ≈ squareQ(q)'
                    @test q*l ≈ a
                    @test a*(qla\b) ≈ b atol=3000ε
                    @test Array(qla) ≈ a
                    sq = size(q.factors, 2)
                    @test *(Matrix{eltyb}(I, sq, sq), adjoint(q)) * squareQ(q) ≈ Matrix(I, sq, sq) atol=5000ε
                    if eltya != Int
                        @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab}, q)
                        ac = copy(a)
                        @test ql!(a[:, 1:5])\b == ql!(view(ac, :, 1:5))\b
                    end
                    qlstring = sprint((t, s) -> show(t, "text/plain", s), qla)
                    rstring  = sprint((t, s) -> show(t, "text/plain", s), l)
                    qstring  = sprint((t, s) -> show(t, "text/plain", s), q)
                    @test qlstring == "$(summary(qla))\nQ factor:\n$qstring\nL factor:\n$rstring"
                end
                @testset "Thin QL decomposition (without pivoting)" begin
                    qla   = @inferred ql(a[:, 1:n1], Val(false))
                    @inferred ql(a[:, 1:n1], Val(false))
                    q,l   = qla.Q, qla.L
                    @test_throws ErrorException qla.Z
                    @test q'*squareQ(q) ≈ Matrix(I, a_1, a_1)
                    @test q'*rectangularQ(q) ≈ Matrix(I, a_1, n1)
                    @test q*l ≈ a[:, 1:n1]
                    @test q*b[1:n1] ≈ rectangularQ(q)*b[1:n1] atol=100ε
                    @test q*b ≈ squareQ(q)*b atol=100ε
                    @test_throws DimensionMismatch q*b[1:n1 + 1]
                    @test_throws DimensionMismatch b[1:n1 + 1]*q'
                    sq = size(q.factors, 1)
                    @test *(LowerTriangular(Matrix{eltyb}(I, sq, sq)), adjoint(q))*squareQ(q) ≈ Matrix(I, a_1, a_1) atol=5000ε
                    if eltya != Int
                        @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab},q)
                    end
                end
            end
            if eltya != Int
                @testset "Matmul with QL factorizations" begin
                    a = raw_a

                    qla = ql(a[:,1:n1], Val(false))
                    q, l = qla.Q, qla.L
                    @test rmul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
                    @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
                    @test rmul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
                    @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),adjoint(q))
                    @test_throws ErrorException size(q,-1)
                    @test_throws DimensionMismatch q * Matrix{Int8}(I, n+4, n+4)
                end
            end
        end
    end

    @testset "transpose errors" begin
        @test_throws MethodError transpose(ql(randn(3,3)))
        @test_throws MethodError adjoint(ql(randn(3,3)))
        @test_throws MethodError transpose(ql(randn(3,3), Val(false)))
        @test_throws MethodError adjoint(ql(randn(3,3), Val(false)))
        @test_throws MethodError transpose(ql(big.(randn(3,3))))
        @test_throws MethodError adjoint(ql(big.(randn(3,3))))
    end

    @testset "Issue 7304" begin
        A = [-√.5 -√.5; -√.5 √.5]
        Q = rectangularQ(ql(A).Q)
        @test norm(A+Q) < eps()
    end

    @testset "ql on AbstractVector" begin
        vl = [3.0, 4.0]
        for Tl in (Float32, Float64)
            for T in (Tl, Complex{Tl})
                v = convert(Vector{T}, vl)
                nv, nm = ql(v)
                @test nv*nm ≈ v
            end
        end
    end

    @testset "Issue 24589. Promotion of rational matrices" begin
        A = rand(1//1:5//5, 4,3)
        @test first(ql(A)) == first(ql(float(A)))
    end

    @testset "Issue Test Factorization fallbacks for rectangular problems" begin
        A  = randn(3,2)
        Ac = copy(A')
        b  = randn(3)
        b0 = copy(b)
        c  = randn(2)
        @test_broken A \b ≈ ldiv!(c, ql(A ), b)
    end

    @testset "Wide QL" begin
        A = randn(3,5)
        Q,L = ql(A)
        @test Q*L ≈ A
    end

    @testset "lmul!/rmul!" begin
        A = randn(100,100)
        Q,L = ql(A)
        x = randn(100)
        b = randn(100,2)
        @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
        @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
        @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
        @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
        c = randn(2,100)
        @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
        @test rmul!(copy(c), Q') ≈ c*Matrix(Q')

        A = randn(100,103)
        Q,L = ql(A)
        x = randn(100)
        b = randn(100,2)
        @test lmul!(Q, copy(x)) ≈ Matrix(Q)*x
        @test lmul!(Q, copy(b)) ≈ Matrix(Q)*b
        @test lmul!(Q', copy(x)) ≈ Matrix(Q)'*x
        @test lmul!(Q', copy(b)) ≈ Matrix(Q)'*b
        c = randn(2,100)
        @test rmul!(copy(c), Q) ≈ c*Matrix(Q)
        @test rmul!(copy(c), Q') ≈ c*Matrix(Q')
    end
end

# choleskyinv.jl
@testset "choleskyinv" begin
    etol = 1e-9
    n = 20
	# real matrices
    A = randn(n, n)
    P = A*A'
    Pi = inv(P)
    C = choleskyinv!(copy(Matrix(P)))
	@test(norm(C.c.L*C.c.U-P)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
    C = choleskyinv(P)
	@test(norm(C.c.L*C.c.U-P)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)

	# repeat the test for complex matrices
    A=randn(ComplexF64, n, n)
    P=A*A'
    Pi=inv(P)
	C=choleskyinv!(copy(Matrix(P)))
	@test(norm(C.c.L*C.c.U-P)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
    C = choleskyinv(P)
	@test(norm(C.c.L*C.c.U-P)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Pi)/√n < etol)
end


text:
#!/usr/bin/env julia

# Usage
#
# You can execute this file directly as long as the `julia` executable is
# in your PATH and if you are on a Linux/macOS system. If `julia` is not in
# your PATH or if you are on a Windows system, call Julia and explicitly
# provide this file as a command line parameter, e.g., `path/to/julia
# Trixi.jl/utils/install.jl`.

import Pkg

# Get Trixi root directory
trixi_root_dir = dirname(@__DIR__)

# Install Trixi dependencies
println("*"^80)
println("Installing dependencies for Trixi...")
Pkg.activate(trixi_root_dir)
Pkg.instantiate()

# Install Trixi2Img dependencies
println("*"^80)
println("Installing dependencies for Trixi2Img...")
Pkg.activate(joinpath(trixi_root_dir, "postprocessing", "pkg", "Trixi2Img"))
Pkg.instantiate()

# Install Trixi2Vtk dependencies
println("*"^80)
println("Installing dependencies for Trixi2Vtk...")
Pkg.activate(joinpath(trixi_root_dir, "postprocessing", "pkg", "Trixi2Vtk"))
Pkg.instantiate()

println("*"^80)
println("Done.")