- TODO
- Haz una cadena de
a
repetida 100 veces.
"a"^100
- Construye un vector de 100
a
s. Considera usar unacomprehension
.
['a' for i in 1:100] == fill('a', 100)
Una matriz se puede hacer
[1 for i in 1:3, j in 1:3]
- Cuanto es
2^100
?
big(2)^100
Nota: No usen Floats
para guardar monedas, usen punto fijo como DecFP.jl
4. Crea una cadena de la 'a'
a la 'z'
.
'a':'z'
|> collect
|> join
- Intenta sumar
"a" + "b"
y'a' + 'b'
. Ahora intenta restar'A' - '0'
. Intenta multiplicar cadenas. Que sale? Nota.Cuidado con eltype piracy
- definir funciones o metodos sobre tipos que no definieron ustedes. Esto en general es muy mal visto; evitar.
'a' - '0' # Espacio afin, como las horas
- Cual es el
Int64
mas grande? ElInt32
? Y todos los Ints nativos?
typemax(Int64) == 9223372036854775807
- Cual es el
Float32
mas pequeno? Cual es la division mas pequena entreFloat32
? - Como obtienes un
Inf
? Cuantos hay? Cuantos0
s hay?
1/0 # -1/0
-0.0 === 0.0
- Haz un arreglo con todos los numeros pares del 1 al 100.
2:2:100 |> collect
[i for i in 1:100 if iseven(i)]
- Haz un arreglo con todos los impares del 1 al 100.
[1:2:100...]
- Haz un arreglo con todos los numeros cuadrados del 1 al 100 (incluyendo 100) Nota: Definimos funciones en Julia de al menos 3 maneras:
# La definicion en bloque
function nombre(param1, param2)
if param1 < 3
return param1 + param2
end
param2
end
# La definicion "inline"
nombre(param1, param2) = param1 + param2
# definicion anonima
function (param1, param2)
param1 + param2
end
# Otra funcion anonima
x -> x + 1
# funcion parcial
# ==(3)
issquare(n) = isinteger(sqrt(n))
[i for i in 1:100 if issquare(i)]
(1:10) .^ 2
- Haz una matriz de 10x10 de solo booleanos.
rand(Bool, 10, 10)
- Haz un vector de 10 flotantes entre 0 y 1 aleatorios
[rand() for i in 1:10]
-
Define una funcion para:
- sumar 1
sq = (1:10) .^ 2 f(x) = x + 1 f.(sq)
- sumar todos los numeros del 1 al 100
function cien() cuenta = 0 for i in 1:100 cuenta += i end cuenta end
- sumar todos los pares del 1 al 100
sum(2:2:100)
- sumar todos los impares del 1 al 100
sum(1:2:99)
- Imprimr
"Cadena"
si le dan una cadena y"Numero"
si le dan un numero
function g(x) if isa(x, Number) println("NΓΊmero") elseif isa(x, String) println("Cadena") end end # Operador ternario "?" # condicional ? caso1 : caso2 g(x) = isa(x, Number) ? println("num") : println("cadena")
- Aplica la funcion de
sumar 1
a un arreglo.
f.(x) map(x -> x + 1, 1:100) # Los impares filter(isodd, 1:100) # lo mismo que una suma, pero el paralelo reduce(+, 1:100)
Nota: Esto se conoce como "Higer order functions" - funciones que aceptan funciones. Son MUY utiles, y vale la pena conocer:
map
,filter
,reduce
,findfirst
,count
,foldr
7. Encontrar el primer numero divisible entre 7 despues entre el 503 y el 520.findfirst(n -> n % 7 == 0, 503:520) 8. Encuentra el primer indice del numero igual a 42 en un arreglo de numeros aleatorios de tamanio 1M si todos estan entre el 1 y el 100. ```julia findfirst(==(42), rand(1:100, 1_000_000))
- Que hacen las funcioens
any
,all
? Como las puedes combinar con un generador? Y con tests?
any(i for i in (true, true, false)) @test all(iseven, sq)
- Que es
infix
? Que esprefix
? Como puedes definir una funcion parcial prefix que revise si algo es"igual a 3"
. Infix es como3 + 2
prefix es como `+(3, 2) - Contar el numero de
#
en la siguiente cadena: (define una funcion e ignora las comillas)
"?#?##?#?#???####??##?#??#?#?#?###??#?##??#?##??##?#?#???###??#???#??####??##?##???##?#??#?##????##?#"
count('#', str) count(str .== '#')
Nota: No hay
NULL
como en R de Julia. 12. Calcular la distancia de Hamming 1. La distancia de Hamming es la cuenta de cuantos caracteres difieron entre una cadena y otra (asumiendo que tienen la misma cantidad de caracteres):AAABBB
yAAACCC
tienen una distancia de Hamming de3
. ``julia function hamming(a, b) cuenta = 0 for i in 1:length(a) if a[i] != b[i] cuenta = cuenta + 1 end end cuenta end
Temario:
lstopo
y latencia- Instalacion de VSCode y el plugin
- Structs y parametricos
- Vectores, Algebra lineal, Principio de optimizaciones
- Instalacion via juliaup - USEN ESTO!!!
function dead() # Dead code optimization
#x = 1 # NO LO NECESITO
#x = 2 # TAMPOCO
x = 3
y = 1
x + y # Register allocation/saving
end
# Loop Invariant Code Motion -> "saca" una operacion del Loop
function licm(vector, n) # v = [1,2,3], n = 3
avg = 0
for i in vector
avg += i
end
avg/n
end
function branch(x) # branch elision / Branch compaction
if x < 5
if x < 10
return true
end
end
end
- Define un struct para un punto en 3D, con un campo por coordenada.
Hint: Se va a ver algo como
Usas
struct Point3D ... end
Int
s oFloat
s? Escoge cualquiera de los dos por ahora. - Define una funcion que calcule la norma del
Punto3D
. - Define una funcion para sumar
Punto3D
s. - Define una funcion para restar
Punto3D
s. - Define una funcion que cree un Punto en 3D de una tupla, y una de 3 numeros sueltos.
- Ahora define un Punto en
Int
s oFloat
s, el que no hayas usado. - Define una funcion que pueda sumar ambos tipos de
Punto3D
. - Escribe una suite de tests para verificar la funcionalidad de tu paquete
- Instala
PkgTemplates.jl
con]add PkgTemplates
(o bienusing Pkg; Pkg.add("PkgTemplates.jl")
). Busca en la documentacion como hacer un nuevo paquete. (Esto puede tardar un rato si no tienes perfil degit
. Si no es asi, usaPkg.generate("Punto3D")
, y genera untest/runtests.jl
a manita.) - Genera un nuevo paquete para
Punto3D
y mete los tests que escribiste anteriormente entests/runtests.jl
- Exporta tus funciones para hacer un
Punto3D
. - Activa un ambiente para
"Punto3D"
con]activate
- tu consola deberia verse azul y cambiar de nombre a"Punto3D>"
. Agrega el paquete deStaticArrays
dentro de este ambiente. Si lo has hecho correctamente, deberias ver que]st
incluye aStaticArrays
en el output. - Agrega funcionalidad a
Punto3D
para que acepteStaticArrays
y tests. - Benchmarkea la diferencia con
BenchmarkTools.jl
.
Repaso
- Que es la
estabilidad de tipos
? Por que es importante? Intenta evaluarsqrt(-9)
y explicar por que no da unComplex64
.
La estabilidad de tipos permite que Julia entienda la entrada y salida de tipos de datos y especialize el codigo de tu funcion con base en ello.
`sqrt(-9)` romperia esa estabilidad si se permitiera en regresar `Float64` en algunos inputs y `Complex64` en otros.
- Intenta definir
- un
struct
dePunto3D
con 3 camposx,y,z
. - Define
Base.:+(p1::Punto3D, p2::Punto3D)
coordenada a coordenada. - Ahora define una funcion que sume un arreglo:
julia function suma(v::Vector{Punto3D}) res = 0 for i in v res += i end res #... end
Intenta usar sumar el siguiente arreglo con tu funcion:julia v = [Punto3D(1,1,1) for i in 1:1_000_000]
- Ahora intenta hacer que
Punto3D
funcione para cualquierInt
yFloat64
nativo. LlamaloP3D{T}
. - Inspecciona los tipos y ensamblador de
Punto3D
vs una version que no tengan tipos losx,y,z
. 1. Recuerda que@code_warntype
te alerta si Julia no puede inferir (deducir) el tipo de datos que sale de tu funcion. Si falla, intenta usarCthulhu.jl
. 2.@code_llvm misuma(v)
va a decir cosas interesantes - que relaciones vez con@code_native misuma(v)
y@btime
deBenchmarkTools.jl
cuando no usas tipos en los campos dex,y,z
vs cuando si? Y cuando los parametrizas? - Bonus: Intenta usar
StructArrays
.
- La matriz de Strang se de fine como una matriz cuadrada que tiene un
2
en la diagonal,-1
, en las diagonales no centradas (off-diagonal) y 0 en todos los demas lugares - define una funcion para crear una matriz de Strang para un tamano
n
.
2 -1 0
-1 2 -1
0 -1 2
Solucion primera:
function strang(n)
mat = [0 for i in 1:n, j in 1:n]
for i in 1:n
for j in 1:n
if i == j
mat[i, j] = 2
end
end
end
for i in 1:n
for j in 1:n
if abs(i - j) == 1
mat[i, j] = -1
end
end
end
mat
end
- Ahora hazla generica usando las funciones
zeros
,one
. BONUS: Puedes definirla en una sola linea?
? :D - Define una funcion
suertuda(T, n)
dondeT
sea un tipo de dato numerico yn
sea el tamano cuadrado de la matriz, y todas las entradas valgan7
en ese tipo numerico. - Convierte la Matriz de Strang en una matriz
TriDiagonal
. - Genera una matriz de con ~10% de las entradas siendo 0, y el resto siendo un numero flotante entre 0 y 1 de tamanio 1000x1000 - no uses elementos cuyos indices difieran en mas de 10 unidades.
a = [rand() < 0.001 && abs(i - j) <= 10 ? 1 : 0 for i in 1:1000, j in 1:1000]
- Cuadrala y toma el tiempo con
@time
. - Conviertela en una matriz
SparseArray
y cuadrala, mide las diferencias de tiempos con@time
. - Cuantos ordenes de magnitued mas rapido fue uno que el otro?
- Multiplica
[1,2,3] .* [1 2 3]
y[1,2,3] .* [1,2,3]
. Por que son distintos? Cuales son los tipos de datos de salida? - Usa
@edit [1,2,3] .* [1,2,3]
y@edit [1,2,3] .* [1 2 3]
para verificar tus resultados. - Que es un
CartesianIndex
? Apoyate en el manual - Que es
IndexStyle
?
De ahora en adelante queremos usar `eachindex`.
- Como puedes accessar las factorizaciones de una matriz? Cuantas hay en
LinearAlgebra
? - Define una matriz de 2x2x2 con los numeros del 1 al 8.
- Que es un numero de condicion?
- Investiga el libro de Gilbert Strang de Algebra lineal para ingenieros, pues la rama de
Analisis numerico
yalgebra lineal numerico
es parte de esta disciplina.
- Siempre hay un poquito "mas" que investigar/hacer/calcular y ahi se rompe al abstraccion de tus datos. El chiste es que tu herramienta (Ya sea Julia o algo mas) te permite construir eso facilmente y que sea rapido.
Arreglos:
x = [1,2,3]
x[begin] == 1
x[end] == 3
x[3] # x[begin:end -1] es un idioma util
- Descarga
UnicodePlots.jl
. Usa 2 vectores aleatorios de 100 para intentar varias graficas. - Descarga
Plots, DataFrames, CSV, FileTrees
. Grafica alguna grafica de que hayas hecho conUnicodePlots
pero ahora enPlots
.
- Agrega una leyenda
- cambia los limites
- agrega un titulo
- guardalo en formato
.png
. - cambia el backend a
plotlyjs
e interactua un rato.
- Por que preferirias
Makie.jl
en vez dePlots.jl
?
- Corre un demo de
Makie
en 3D y juega con los parametros. - Visita la
Makie gallery
y encuentra algo cool.
- Que es un API? Como lo puedes usar para descargar un CSV?
Para nosotros, un API para descargar datos simplemente pide una llamada de funcion con HTTP.jl y permite descargar datos.
Se pueden ver com o
using HTTP
HTTP.get("https://datos.gob.mx/covid/salud/", key = "1234")
# using Downloads
NOTA: Gran acordeon de DataFrames.jl
- Interned Strings y
String15
- Trata de homogeneizar los tipos de datos de tus columnas: si tienes muchos
Union{X,Y,Z,A,B,...}
el compilador de Julia va a sufrir y ademas va a ser lento el performance en run time.
- Como puedes cargar un
CSV
en Julia? Intenta con el ejemplo del INEGI aqui. UnCSV
es un archivo de valores separados por comas, ieComma Separated Values
. - Ejemplo
FileTrees.jl
yGlob.jl
- define una funcion para contar todas las palabas separadas por
" "
en una cadena. - aplicala a todo un directorio de manera recursiva con
FileTrees.jl
. Que solo cuente las palabras
- Cuando convienen
DataFrames.jl
y cuando convieneSQL
(u otro tipo de base de datos) + la "ventaja Karaminski" en Julia.
SELECT *
FROM gente
WHERE gente >= 10
- SQL:
- permite guardar tus datos en casos de catastrofe
- optimiza la cantidad de veces que uno "va a disco"
- tiene medidas para prevenir corrupcion de datos
- medidas de seguridad para encriptar los archivos y restringir acceso NOTA: checar el curso de Andy Pavlo de CMU
- TODO
- Mide antes de optimizar! Leer este paper de COST
Puedes tener varios cores hasta que demuestres que sabes usar uno solo
- Demo: usando el profiler de VSCode
# Intenta correr este snippet de codigo dentro de VSCode
function profile_test(n)
for i = 1:n
A = randn(100,100,20)
m = maximum(A)
Am = mapslices(sum, A; dims=2)
B = A[:,:,5]
Bsort = mapslices(sort, B; dims=1)
b = rand(100)
C = B.*b
end
end
# compilacion
@profview profile_test(1)
# runtime puro
@profview profile_test(10)
- Al menos puedes usar
@time
/@btime
- Daniel y INEGI (DataFrames.jl + StasBase.jl)
- Brenda y optimizacion en redes/grafos (Graphs.jl antes LightGraphs.jl)
- Jesica y NLP
- Arturo y procesamiento de covid en batch/paralelo (Lecture de
parallel prefix sum
mas abajo yFileTrees.jl
)
Paralelismo a que nivel? NOTA: Breviario cultural - PaSH scripts de shell en paralelo Este es un gran lecture de donde se inspiran estas notas.
- Inter Word Parallelism, Super Word Parallelism
01 01 01 01 01 01 01 01
- Instruction pipelining +
ports
: muchas cosas pueden operar al mismo tiempo
v = v + 1
x = x + 1
y = y + 1
z = z + 1
- Atomics
x = 1
y = x + 1
NOTA: No se preocupen todavia por usar atomics
, mutex
, o locks
, syncronizacion
NOTA: Seguridad
- Multicore (usar muchos cores en una sola compu) : en 1 sola compu, usar varios cores. πππ Este usaremos! πππ
- Hyperthreading (no lo vamos a usar)
- GPU/aceleradores (arquitectura heterogenea) -> π Este usaremos! π
- Distribuido (muchas compus enchufadas) -> Muchas compus trabajando en el mismo problema
Distributed.jl
. Cual creen que nos vamos a enfocar? Multicore! π - Concurrencia 1:N: 1 trabajador puede
- Tiende a servir cuando esperas a los demas, ie, I/O y si el procesamiento es ligero.
- Contencion: que tanto distinto trabajadores/cores estan peleando por accesar un recurso.
- Paralelismo y teoria
- Ejemplo a mano
- Camino critico
- Pregunta: Cual es la diferencia entre
Camino critico
yseccion critica
? Respuesta: Laseccion critica
tiene que ver con cuando queremos garantizar que exista acceso unico a un recurso. Elcamino critico
es el cuello de botella serial de un algoritmo.
- Pregunta: Cual es la diferencia entre
- Span
- Paralelismo
- Cilk, Futhark y multiversiones
- "Task parallel, depth first, randomized work stealing scheduler"
- Task parallel: "fork join"
import Base.Threads.@spawn function fib(n::Int) if n < 2 return n end t = @spawn fib(n - 2) return fib(n - 1) + fetch(t) end
- depth first:
- randomized work stealing:
- scheduler:
Base.@threads
: no es "composable", ie, no puedo usar@threads
dentro de otro@threads
. Sirve cuanto solo quieres dividir la chamba (que es igual para todos) entre todos tus cores.ThreadsX.jl
:Folds.jl
:- usando y sin usar
Folds.jl
, trata de resolver los siguientes ejercicios en tu equipo:- Genera una cadena de ~ 1M de caracteres aleatoriamente escogidos entre
!@#$
y cuenta cuantos son iguales a#
. - Encuentra la primera secuencia igual a
!@#$
. Encuentra las entradas de indice igual a!@#$
- Agrega 1 a un vector de de enteros de tamano 1M, con 20 tamanos logaritmicamente espaciados entre 100 y 1M. Cuando gana la version serial? Cuando la paralela? Que tamanio debe tener el vector para que tu benchmarking deje de ser superlineal?
- Genera una cadena de ~ 1M de caracteres aleatoriamente escogidos entre
- usando y sin usar
- Cuando te conviene
pmap
y cuando@threads
? - Ver la siguiente lecture y encontrar 3 ejemplos en donde se aplican algoritmos tipo
parallel prefix sum
.
RESUMEN: No hay una sola respuesta de como resolver el paralelismo (aun). Esta presentacion del creador de Futhark es un buen resumen de por que es tan dificil hacer bien paralelismo.
- Julia comparte el mismo tipo de tecnologia que
IntelTBB/Cilk
yRayon
de Rust con el scheduler, cuando usas@spawn
. El scheduler deGo
usa un algoritmo debreadth first
, en vez dedepth first
, el cual optimiza por casos de mucho I/O***. - Julia comparte el mismo tipo de tecnologia (hoy) que
OpenMP
conBase.@threads
.
Recuerda que siempre puedes ver workshops/talleres de Julia en el canal de YouTube de Julia y leer la parte del manual que te interesa
- Como defines una expresion en Julia? Para que usarias
quote begin ... end
?
ex = :(x = 1)
ex = quote x = 1 end
- Define una expresion
x = 1
. Cambia el valor de1
a42
accesando la expresion. Evaluala.
ex = :(x = 1) # recuerda Meta.@dump x = 1
ex.args[2] = 42
ex
- Que es un macro? Que restricciones tiene? Cuando corre?
Un macro es un tipo de funcion especial que
1. recibe una `Expr`
2. la manipula
3. regresa otra `Expr`
4. Corre al tiempo del `parser`
- Que son los macros
_str
? Como usarias uno? Construge una expresion regular, un vector de bites, y una cadena literal con sus macros respectivos (recuerda que puedes usar?"pancho"
en el prompt para buscar a pancho en los docstrings de Julia) - Que es la higiene en los macros?
Los macros definen variables que se resuelven (o valen) en un `module` distinto al global, es suyo propio, para evitar el choque de nombres.
- (Dificil, salta esto y regresa al rato) Usa
MacroTools.jl
para definir un ejemplo de arboles de LindenMayer (No olvides usarMeta.@dump
y@macroexpand
).
- Un arbol de linden se ve asi:
- Nota: Nunca usar
Base.@pure
. - Entender que casi no vale la pena escribir macros en Julia general - usa funciones!
- BONUS: Usa
SnoopCompile.jl
yCthulhu.jl
para
- Define la siguiente funcion, y corre los siguientes macros:
f(x) = x^2
Meta.@dump f(3)
# ASTMeta.@lower f(3)
# Lowered@code_lowered f(3)
@code_warntype f(3)
@code_typed f(3)
# Typed@code_llvm f(3)
# Optimizando...- (Aqui sucede la magia de optimizaciones de LLVM)
@code_native f(3)
- Trata de dar ejemplos para usar los siguientes macros
@enter
, deDebugger.jl
@warn
- arroja una advertencia@show
@info
@log
@which
- ensenia que metodo aplica@which 3+ 3
@doc
- ensenia docstrings y los puede anadir a una funcion@elapsed
- cuanto tiempo ha transcurrido@inbounds
- quita elbounds checking
@simd
- ~trata de forzar la vectorizacion de tu codigo, usa en vezLoopVectorization.jl
y su@turbo
.@nexprs
- bendito sea Cthulhu que ya no escriboFortran
!@nloops
- algebra multidimensional lineal!@views
- para cuando usas slicesx[:] = a[:] + b[:]
del lado derecho de una igualdad, ahorra las copias.@.
- Broadcast toda una expresion sin tener que usar.*
todo el tiempo.@test_broken
,@test_throws
- Define una expresion de
for
y usapush!
para definir algo interesante
- Para que su codigo sea vea re-chulo, usen JuliaMono y pongan las
ligatures
en activo.
- Dudas y resolucion de problemas
Hoy vamos a seguir (tanto como podamos) el workshop de LightGraphs.jl (Ahora Graphs.jl) encontrado en este repositorio
NOTA: Hay documentacion que menciona LightGraphs.jl
- eso yo no existe, usen Graphs.jl
- deberia ser 100% compatible.
- Que es un nodo? Que es una arista? Que es un sucesor? Que es un vecino? Que es un sucesor? Que es una grafica?
Nodos son los "puntos"
vertices son las lineas que los conectan
Vecino - si 2 nodos estan conectados por 1 arista, son vecinos
Sucesor - Si 1 arista sale de un nodo y llega a otro, al que le llega la arista es el sucesor
Grafica* - {Nodes x Vertices}, en ingles le dicen "Graph = {Nodes x Edges}"
*Grafica simple:
1.5 - Instala Graphs.jl
con using Pkg; Pkg.add("Graphs.jl")
2. Define una grafica simple con 3 nodos y 3 vertices entre ellos.
g = SimpleGraph(0)
for i in 1:3
add_vertex!(g)
end
add_edge!(g, 1 => 2)
add_edge!(g, 2 => 3)
add_edge!(g, 3 => 1)
# Para graficar:
using Random, Cairo, Fontconfig, GraphPlot
g = complete_bipartite_graph(2, 2)
Random.seed!(42)
gplot(g, nodelabel=vertices(g))
- Muestra su matriz de adyacencia y su espectro de adjacencia
adjacency_matrix(g)
adjancency_spectrum(g)
- Grafica una grafica bipartita con tamanio
(2,2)
aleatoreamente generada.
g = complete_bipartite_graph(2, 2)
- Usa
Base.summarysize
para saber cuanto pesa unVector{Int}
vacio y unSet{Int}
vacio.
Base.summarysize(Int[]) == 40
Base.summarysize(Set{Int}()) == 336
- Compara el tiempo de creacion de una grafica
SimpleGraph
de tamanios10
a10000
en potencias de10
.
for i in (10, 100, 1000)
@btime SimpleGraph(i)
end
- Genera una grafica complete de tamanio
10
.
g = complete_graph(10)
- Que es una grafica dirigida? Genera de tamanio 13, y graficala.
Una grafica dirigida es una grafica en donde la direccion de las aristas importa: si sale una arista de 1 a 2, se dice que 1 es el antecesor de 2, y 2 es el sucesor de 1. Tambien puedes tener una arista que salga de 2 a 1.
g = SimpleDiGraph(5)
add_edge!(g, 1 => 2)
add_edge!(g, 1 => 3)
add_edge!(g, 3 => 1)
# tip: intenta agregar un nodo a si mismo
- Describe los campos de una grafica simple
SimpleDiGraph{Int}
. Por que hay redundancia?
Para favorecer ciertos casos/iteraciones algoritmicamente via despacho multiple -
- Que es mas facil: recorrer todos los sucesores de un nodo o todos sus sucesores?
Depende el sapo la pedrada.
`g.fadjlist` es para "forward adjacency list" y `g.badjlist` es para "backward adjacency list"
- Genera una grafica ciclica de tamanio 5. Plottea los casos de todos los ejemplos de
complete_*
.
g = complete_graph(5)
gplot(g, nodelabel = vertices(g))
g = complete_bipartite_graph(3,5)
gplot(g, nodelabel = vertices(g))
g = complete_digraph(5)
gplot(g, nodelabel = vertices(g))
g = complete_multipartite_graph(1:5)
gplot(g, nodelabel = vertices(g))
- Construye una grafica dirigida simple con la matriz
0 1 1 0
0 0 1 0
0 0 0 0
0 0 1 0
gplot(SimpleGraph([
0 1 1 0
0 0 1 0
0 0 0 0
0 0 1 0
]))
- Genera una grafica lineal de un iterador por medio de
SimpleGraphFromIterator
y graficala.
```julia
iter = (Edge(i, i+1) for i in ???)
```julia
iter = (Edge(i, i+1) for i in 1:4)
g = SimpleGraphFromIterator(iter)
gplot(g, nodelabels = vertices(g))
- Da ejemplos de las funciones
nv
,neighbors
, - Calcula los componentes conectados de una grafica
g
de tu interes. Ahora usa losstrongly_connected_components
en una grafica dirigida - Investiga si alguno de los algoritmos en
Graphs.Parallel.<TAB>
te sirven para agilizar un codigo.
- Observa todos los tipos de graficas en
smallgraph
. Recolecta sus aristas conedges(g) |> collect
.
gplot(smallgraph("house"))
gplot(smallgraph("petersen"))
gplot(smallgraph("tutte"))
gplot(smallgraph("sedgewickmaze"))
gplot(smallgraph(:karate))
- Que es un camino? Que es un camino Hamiltoniano? Que es una grafica aciclica? Que es un arbol? Que es un
arbol de expansion minimo
? Calcula usandoprimm_mst
Un camino es una sucesion de nodos conectados por aristas dirigidas en donde solo se sigue la direccion de la arista.
Una grafica aciclica es una grafica en donde todas los nodos nunca pueden visitarse a si mismos siguiendo todos los posibles caminos que salen de ellos mismos.
Un camino Hamiltoniano es un camino que puede visitar todos los nodos de una grafica conexa sin repetir nodos.
Un arbol es una grafica aciclica en donde 1 solo nodo (llamado raiz) puede visitar a todos los demas nodos por medio de caminos
- Encuentra como llamar los siguientes algoritmos en una grafica de tu interes:
- Dijkstra
- A*
- Prim
- Bellman-Ford
- Floyd-Warshall
- Carga un dataset con
GraphsIO.jl
ySNAPDatasets.jl
. Corre un algoritmo de tu interes y benchmarkealo - Define
g = smallgraph(:diamond)
y guardalo conGraphIO.savegraph
(busca en los docs como se usa) - Que es
GraphBLAS
? Que esSuiteSparseGraphBLAS
? Como puedes usarlo en Julia? NOTA: BLAS significa "Basic Linear Algebra Subroutines" y existen desde hace 6 anios y todos las usan.
# Trabajo de Dr. Tim Davis (experto en algebra lineal sparse) y Will Kimmerer es el que lo ha traido a Julia
using SuiteSparseGraphBLAS
Dia 08: Computo Simbolico, Ecuaciones diferenciales con DifferentialEquations.jl y DeepLearning con Flux.jl
- Ada Lovelace, tejidos, proto ensamblador
- encuesta
- Proyecto de Arturo, mencionar Arrow.jl
- Brenda?
- Noel y Laplace?
- Jesica y cargar datos?
- Learning from Data de Gilbert Strang
- Parallel Computing and Scientific Machine Learning
- Introduction to Computational Thinking
- Instala
Symbolics.jl
conusing Pkg; Pkg.add("Symbolics")
. - Define unas dos variables
x
yy
simbolicamente. Usa Julia para confirma las formulas de sumas y diferencia de cuadrados:
@vars x y
ex = (x - y) * (x + y)
simplify(ex, expand = true)
(x - y) * (x + y) == x^2 - (y^2)
(x + y) * (x * y) == x^2 + y^2 + 2x*y
simplify((x+y)^2, expand = true)
- Usa
Symbolics.jl
para simplificar la espresionsin^2(x) + cos^2(x)
.
ex = cos(x)^2 + sin(x)^2
simplify(ex)
- Define una funcion como
f(x) = 3x^2
y toma su derivada
@vars x
D = Differential(x)
expand_derivatives(3x^2)
- Como es distinto el computo simbolico de usar
BigInt
s?
Symbolics trata de representar el concepto matematico de un numero, no su implementacion per se.
- Usa
Latexify.jl
para imprimir una matriz simbolica y pegala en otro lugar usandoclipboard
.
@vars x
m = [1 cos(x); sin(x) -1]
latexify(m, clipboard = true)
- BONUS: Investiga como usaron
Symbolics.jl
para su tecnica de "sparsity detection".. Moraleja - combinar el mundo numerico y simbolico lleva a ventajas algoritmicas que los Fortraneros nunca sonarian. - Como difiere
Symbolics.jl
de SymPy? Y de Mathematica? Lee el manual para saber mas. - Si te late la teoria de categorias, checa Catlab.jl
- Instala
DifferentialEquations.jl
- se va a tardar un rato. - Compara tu solver de ecuaciones diferenciales favorito vs
DifferentialEquations.jl
en esta tabla del Dr. Chris Rackauckas.
- Que es el atractor de Lorenz? Resuelve el atractor de Lorenz con el solver
Tsit5
. Intenta utilizar el macro@ode_def
.
using DifferentialEquations
- BONUS plottea varias trayectorias para convencerte de su comportamiento caotico.
- TAREA: Si no te gustan mis preguntas, completa cualquiera de los siguientes tutoriales de SciML.
- Que es el modelo de Lotka-Volterra? Encuentra, resuelve, y plottea el model de Lotka-Volterra. Lo puedes hacer con emojis?
Si se puede, busca en el manual π
- Encuentra la mayor cantidad de Solvers estilo diferencias finitas.
- Que es un modelo "surrogado"? Corre uno de tu interes con Surrogates.jl y compara su perfromance y precision
- Escoge un tutorial de
Flux.jl
y siguelo lo mejor que puedas. Los puedes encontrar aqui - Que es la diferenciacion automatica? Como puedo derivar una catapulta?
- Picarle a la documentacion de Flux.jl y llevar un curso de Deep Learning. Volverse mesias tecnologico y comprar la casa en Bahamas.
- Instalar los siguientes paquetes:
using JuMP, BioSequences, FASTX, HiGHS, DataFrames
- Constancias y formulario
- Proyectos finales
- BONUS BONUS: Leer
La triple helice
de Richard Lewontin.
- Instala
BioSequences.jl
conusing Pkg; Pkg.add("BioSequences")
- Si quieres ver un repositorio enorme de problemas de secuencias genomicas, revisa rosalind.info
- Abre la documentacion de BioSequences.jl y contesta las siguientes preguntas
- Cuantos nucleotidos existen?
alphabet(DNA) |> length == 16
alphabet(RNA) |> length == 16
- Cuanto pesan en memoria
"a"
,'a'
,0
,true
? Muestra sus bits si es posible.
Base.summarysize.(["a", 'a', 0, true]) == [9, 4, 8, 1] # Recuerda que estoy mide en bytes! O sea, mulitplica por 8 y es el numero de bits
- Ahora busca la representacion en memoria de adenina de RNA y adenina de DNA
Base.summarysize(RNA_A) == Base.summarysize(DNA_A) == 8
- Convierte
'C'
en citosina de ADN.
convert(DNA_C, 'C')
- Cuanto pesa en memoria la leucina
AA_L
? Como esta representada? Como esta representado'a'
?
Base.summarysize(AA_L) == 1
# Si usamos @edit AminoAcid('L'), y buscamos la tabla... y encontramos que AA_L es 0x0a
bitsring(0x0a) == "00001010"
- Usa
LongDNA{4}(...)
para convertir la cadena"TTANC"
. Que cosas interesantes puedes decir de sus estructura? Y de sus supertipos? Que metodos puedes llamar sobre ese objeto?
LongDNA{4}("TTANC")
- Concatena dos cadenas de ADN
"ACGT"
juntas (Osea, que te resulte una cadena de ADN "ACGTACGT"`)
LongDNA{4}("ACGT") * LongDNA{4}("ACGT")
- Genera una secuencia de
"ACGT"
de tamanio100
.
repeat(LongDNA{4}("ACGT"), 25)
- Usa los macros
dna"..."
yrna"..."
para definir cadenas de tamanio 10. Define el tuyo para moder hacermidna"ACGT"
y que funciones igual queLongDNA{4}("ACGT")
dna"ACGTACGTAC"
rna"AAUUUGCUCA"
macro midna_str(x)
LongDNA{4}(x)
end
midna"ACGT"
- Usa los siguientes algoritmos comunes en tus ADNs:
pop!
,popfirst!
pop!(dna"ACGT") == dna"T"
popfirst!(dna"ACGT") == dna"A"
push!
,pushfirst!
push!(dna"ACGT", DNA_T) == dna"ACGTT"
pushfirst!(dna"ACGT", DNA_T) == dna"TACGT"
reverse
,reverse!
reverse(dna"ACGT") == dna"TGCA"
insert!
insert!(dna"ACGT", 2, DNA_T) == dna"ATCGT"
deleteat!
deleteat!(dna"ACGT", 2) == dna"AGT"
- Revisa con la secuencia de ADN de
seq = dna"ACGTTGCAAGCTAGAT"
- si es palindromica
ispalindormic(seq) == false
- obtener su complemento
complement(seq)
- es canonica
iscanonical(seq)
- Descarga la secuencia del covid19 con
FASTX.jl
de este sitio. Lee el archivo.
- encuentra si la secuencia de
"ACCC"
ocurre
query = ExactSearchQuery(aa"ACCC")
occursin(query, covid)
- encuentra la ultima secuencia de
"ACCC"
, la primera, y cuenta todas
query = ExactSearchQuery(aa"ACCC")
findlast(query, covid)
findfirst(query, covid)
count(query, covid) # TODO!
- haz una busqueda aproximada de
"CATCAT"
query = ApproximanteSearchQuery(aa"CATCAT")
findfirst(query, 0, covid)
- Bonus: Cuenta cuantos nucleotidos tipo
'C'
hay en la secuencia del Covid. Recuerda que puedes usar la iteracion estilo
count = 0
for in covidrna
if ...
...
end
end
count(==(AA_C), covid)
O usar metodos funcionales.
- Haz un benchmarking con
@btime
deBenchmarkTools.jl
para contar y comparalo con tu implementacion.
fakecovid = copy(covid)
for i in 1:1000:length(covid)
fakecovid[i] = AA_C
end
matches(fakecovid, covid) == mapreduce(==, +, fakecovid, matches)
- Ahora corre los analisis anteriores pero convirtiendo la secuencia del covid en un
String
convencional. Cual es la diferencia en consumo de memoria? Velocidad? Cuanto pesa uno con respecto al otro en memoria? (Hint: UsaBase.summarysize(x)
) - Inspecciona el ensamblador para ver si hay muchos uso de registros estilo SIMD/vectorizado. (Esos vectores se ven como
xmm
,ymm
,zmm
) - Define una nueva secuencia identica a la de covid, modificala en 10 lugares, y toma la distancia de Hamming. Haz un benchmarking de tu metodo y el de
BioSequences.jl
.
- Llenar el formulario
- Proyectos finales
- dudasquejascomentariossugerencias
- JULIACON 2022! Es gratis! Online! (Mis presentaciones)
- π Julia WATs???? π
- JuMP y optimizacion
Moraleja del curso:
- Nadie tiene tiempo infinito - tus herramientas deben valer la pena contra sus alternativas viables
- Aprender a medir, pedir ayuda, y resolver tus propias dudas con tus propias herramientas - tener "feedback loops" estrechos sinergiza mucho.
- Escoge el algoritmo toericamente optimo -> Minimiza el trabajo que debe hacer tu codigo
- Entiende como se representan las cosas en memoria y como se mueven - no gastes de mas.
- 1 solo core de tu laptop es muy poderoso (puede hacer varios GHz por segundo, o ~10^9 operaciones por segundo), no temas a exprimirle lo maximo y poder aproximar si algo se deberia poder hacer en menos tiempo.
- Discutir los solvers de DifferentialEquations.jl
- Meta:Vamos a resolver juntos el problema de la dieta en JuMP
- Importa las dependencias
using JuMP, DataFrames, HiGHS
- Carga los datos
foods = DataFrames.DataFrame(
[
"hamburger" 2.49 410 24 26 730
"chicken" 2.89 420 32 10 1190
"hot dog" 1.50 560 20 32 1800
"fries" 1.89 380 4 19 270
"macaroni" 2.09 320 12 10 930
"pizza" 1.99 320 15 12 820
"salad" 2.49 320 31 12 1230
"milk" 0.89 100 8 2.5 125
"ice cream" 1.59 330 8 10 180
],
["name", "cost", "calories", "protein", "fat", "sodium"],
)
- Establece las restricciones
limits = DataFrames.DataFrame(
[
"calories" 1800 2200
"protein" 91 Inf
"fat" 0 65
"sodium" 0 1779
],
["name", "min", "max"],
)
- Define un
model
o deJuMP
. Como tenemos un problema lineal, usaremos el optimizadorHiGHS
.
model = Model(HiGHS.Optimizer)
- Que otros solvers para problemas lineales existen? Cuales son gratis? Cual es el que mas conviene? Intenta instalar otro y resolver este problema de la misma manera.
- Define las variables de decision.
@variable(model, x[foods.name] >= 0)
- Que pasa si llamas
@variable model x[foods.name] >= 0
? Son equivalentes?
- Ahora puedes definir una funcion objetivo (que es lo que quieres que cumpla tu funcion principal - esto es lo que buscas "optimizar")
@objective(
model,
Min,
sum(food["cost"] * x[food["name"]] for food in eachrow(foods)),
)
- NOTA: Tambien puedes usar codigo en Julia para definir estas restricciones:
for limit in eachrow(limits)
intake = @expression(
model,
sum(food[limit["name"]] * x[food["name"]] for food in eachrow(foods)),
)
@constraint(model, limit.min <= intake <= limit.max)
end
- Imprime el modelo con
print(model)
- Ahora optimiza y describe el modelo con
optimize!(model)
solution_summary(model)
- Imprime la solucion optima:
for food in foods.name
println(food, " = ", value(x[food]))
end
- Trata de agregar un nuevo
@constraint
(restriccion) al modelo donde
$$
x_{milk} + x_{icecream} \le 6.0
$$
12. Instala NLopt.jl
y reemplaza a HiGHS.Optimizer
como solver.
Repasemos la pagina de los Supported Solvers.
- Cuantos de esos crees poder instalar? Que tan facil es instalarlos? Cuanto tiempo crees que gastarias en instalar eso para tu grupo de alumnos? Tus colaboradores?
- Alguien mas tiene el algoritmo optimo para contar la cantidad de numeros primos mas bajos que un
n
. - No vale la pena reescribir todos los codigos de
C/C++/Fortran/Rust
del mundo en Julia, ergo hay que hacer una interfaz para llamarlo desde Julia, conocido como"wrapper"
(envoltorio, paquete). - Pero esperar que los usuarios puedan instalar las dependencias de terceros (instalar un codigo en
C
, digamos) con todas las combinaciones de compiladores y sistemas operativos (Mac, Ubuntu, Windows, FreeBSD, Arch) es una pesadilla como desarrollador y como usuario (que si no hay instrucciones para tu compu de 32bits π?) - El sueΓ±o βοΈ es que solo tuvieras que hacer
using numerosprimos_jll
- Subir una receta para todos los sistemas operativos
- BinaryBuilder.jl consiste en
- Un ambiente virtualizado ("falso") en donde uno instala las dependencias
- Se generan recetas para todas las arquitecturas
- Se consiguen robots para verificar que las instalaciones en todos los sistemas operativos funcionan
- Suben los binarios a un repositorio central (
Yggdrasil.jl
, el arbol de la vida) - Con hacer
using mipaquete_jll
ya los usuarios lo descarga automagicamente.
- Demo con
Primes_jll
.
using primecount_jll
myprimecount(x) = @ccall libprimecount.primecount_pi(x::Clonglong)::Clonglong
myprimecount(100) == 25
@time myprimecount(1e8)
@time myprimecount(1e14)
- Vamos a
resolver juntosdejar de tarea el problema de Rocket Control
- Resuelve el problema
- Cambia el solver y vuelve a resolver el problema, compara sus tiempos de resolucion y consumo de recursos.
NO NECESITAS ESTO 99% del tiempo: Pero hoy queriamos pavonear π¦ Este es una version un poquito mas interactiva pedagogica que el articulo que traduje en ingles aca, originalmente por Jiahao Chen.
- Ordena los siguiente "pasos" de compilacion en el orden en el que corren:
@code_typed foo(x)
foo(x)
@code_native foo(x)
@time foo(x)
@code_llvm foo(x)
@code_lowered foo(x)
Solucion:
@time foo(x) # <- tiempo del "parse"
@code_lowered foo(x)
@code_typed foo(x) # <- Entre estos dos puntos sucede la recompilacion con funciones `generated`
@code_llvm foo(x) # <-
@code_native foo(x)
foo(x)
- Genera una instruccion de un bucle for
for
sobre la variablei
, de1:n
con cuerpo vacio.
- Llena el cuerpo del for con instrucciones para llenar un arreglo
arr
en las posicionesi
hastai:10
con la funcionf(i)
. (Hint: esto se va a ver algo asi comoex = :(for ...))
, despuespush!(ex.args[...], ...)
). - Ahora define una funcion que si toma un vector
v
de entrada, toma el promedio con los 2 numeros anteriores y los 2 numeros posteriores. Se dice que laventana
es de tamanio 2, y es unrolling average
. (Asume que el vector es de tamanio10
)
- Que optimizaciones podrias hacer si no sabes el tamanio de la ventana de antemano?
- Que es una funcion generada/
generated function
? Cuando corre? Que restricciones tiene?
- Necesitas el `@generated` frente a tu funcion
- Solo puedes operar sobre tipos, no valores
- Debes regresar una `Expr`.
- Solo pueden llamar funciones que ya fueron definidas
- No pueden mutar/observar estado global
BONUS: Leer el paper sobre el world age en Julia
3. Ugh... ok, haremos benchmarking sobre StagedFilters.jl
.
Tarea para Miguel: Yo me comprometo a hacer videocapsulas de los siguientes temas: TODO: GLOSARIO!
- Juliaup y manejando distintas versiones de Julia en una sola maquina
- FileTrees.jl y procesamiento de contar palabras en un directorio
- VSCode plugin: Instalacion, uso, testing, desarrollo
- Revise.jl y workflow
- Makie.jl
- Plots y cambiar de backend
- Programmacion funcional y benchmarking y
- Franklin.jl y blog
- Punto2D, tipos parametricos y performance
- PkgCompiler.jl
- JET.jl, Aqua.jl
- ThreadsX.jl y Folds.jl
- CUDA.jl y
- Matriz de Strang, algebra lineal numerica
- sparse arrays y multiplicacion de ellos
- Compartiendo un MWE en Discourse/Slack/Zulip/Forem?
- JuliaMono y ligatures
rand()
y features delRNG
- Ventaja de Julia: Graphs.jl y GraphBLAS y SparseGraphBLAS
- Jack Dongarra, BLAS benchmarks, Octavian.jl, LoopVectorization.jl
- BLAS subroutines via Julia
- Cuando intersectar/buscar con un vector es mas rapido que en un Set
- Notebook de Jakob Nissen
- Arrow, DataFrames,
- SuiteSparseGraphBLAS example
- Julia y Arduino
- Como hacer un issue en github a un repo de Julia
- Que es semantic versioning?
- Base.summarysize vs sizeof
- Covid19 y BioSequence, covid y evitar el
count_naive
, link de Discourse decount
siendo lento - Artifacts y JLLs - solvers de JuMP y datasets: erradicar fricciones para devs y usuarios
- Leer PDFs y scrapear imagenes con PDFIO.jl
- Tim Holy and the
ProgressMeters.jl