For a fairer comparison:
- Python has mean, median, and geometric_mean functions that could be employed like the Julia library functions.
- There’s no arbitrary recursion depth here.
6 Likes
Thanks.
It was a snip from the 1st comment I saw on Twitter.
Language syntax comparison suffer from the same type of problem as language performance comparison:
it depends on the ability of the person writing the code & on the actual language
The world would be better off if there was a “language elegance game” just like there is a Benchmarks Game.
8 Likes
Wooaah! remarkable! Could you explain how that magic trick works?
Fⁿ(x,n)= ∘(fill(F,n)...)(x)
Thanks!
1 Like
∘ is the function composition function. ∘(f,g,h)(x) turns into f(g(h(x))).
fill(F,n) populates an array with n times the function F. They are composed by ∘ to form the function
x -> F(F(F(F...F(x)))...) .
That’s is, you just call it with any argument you like.
21 Likes
Here’s my seven lines (excluding imports). It calculates and plots pursuit curves using DifferentialEquations.jl. I’ve renamed the pursuer and pursued to fox and rabbit for clarity. It’s kind of fun to play around with different paths for the rabbit to take. I have a parameter, k, which alters the speed of the fox.
The fourth line beginning Rx, Ry,... is a little goofy to fit within the space constraint, but it works and makes what I’m plotting clear as day.
using Plots, DifferentialEquations, LinearAlgebra
rabbit(t) = [cos(t), sin(t)] # Rabbit Running in a Circle
fox(u, k, t) = k * (rabbit(t) - u) / norm(rabbit(t)-u) # Fox chases rabbit at speed k
prob = ODEProblem(fox, [2.0, 0.0], (0.0, 10), 0.8)
sol = solve(prob, saveat=0.1);
Rx, Ry, Fx, Fy = [getindex.(vcat.(rabbit.(sol.t), sol.u), i) for i=1:4]
plot(Rx, Ry, label="Rabbit", lw=sol.t)
plot!(Fx, Fy, label="Fox", lw=sol.t, aspect_ratio=:equal)
Here’s a bonus animation:
29 Likes
Don’t forget one of the most evocative tools in Julia-land, the splat, which is needed to take the array-typed output of fill and use the elements of that array to populate the individual arguments of ∘.
3 Likes
I had to submit another one. I’m still leaning about the ApproxFun.jl package, so I’m happy to hear your feedback. This my first real trial of this package other than running a few of the examples.
The following 7 lines (excluding imports) calculates deflection, angle, moment, and shear of a uniformly loaded cantilevered beam and plots the results. This uses ApproxFun.jl find to solution of an ODE. The solution to this ODE is a function which defines the vertical displacement of the beam (under certain assumptions). Figure of the setup given below.
The line beginning B = ... defines the boundary conditions. It’s a fourth order ODE so there are four boundary conditions. The vertical displacement and first derivative (ie, angle) are both zero at the fixed end of the beam (leftendpoint). The moment and shear in the beam (2nd and 3rd) derivatives are zero at the free end of the beam (rightendpoint).
The solution is generated on the next line, v = .... The boundary conditions, B, are all set to 0 using zeros(4).... The differential equation is defined E*I*Derivative()^4 and gets set to a uniform load one(z) (shown as q in the image above).
The rest is just plotting. It’s ugly again to fit into seven lines.
using ApproxFun, Plots
L, E, I = 12.0, 1.0, 1.0
d = 0..L
z = Fun(identity, d)
B = [[Evaluation(d,leftendpoint,k) for k=0:1]... ; [Evaluation(d,rightendpoint,k) for k=2:3]... ;]
v = [B; E*I*Derivative()^4] \ [ zeros(4)..., one(z)]
func_name = zip([v, v', v'', v'''], ["Deflection", "Angle", "Momement", "Shear"])
plot([plot(z, f, title=n, label="") for (f,n) in func_name]..., lw=3)
Result:
It’s easy enough to change the load function:
v = [B; E*I*Derivative()^4] \ [ zeros(4)..., sin(z)]
I haven’t really messed with the units to ensure the solution is scaled correctly, but it does seem to match the analytical solution of the uniformly loaded beam quite well. The comparison:
plot(z, v, lw=3, label="ApproxFun", title="Deflection | Comparison", legend=:left)
plot!(z-> (z^2)*(z^2 - 4*L*z + 6L^2)/(24*E*I), 0, L, lw=3, ls=:dash, label="Analytical")
19 Likes
- eight or nine easily read lines are better than seven munched up lines
9 Likes
Hey!
I love your code, but you should use
Fn(n) = reduce(\circ, fill(F, n))(x)
...
Fn(100)[1,1,2,3,5]
1 Like
This is not my code, I saw it on Twitter. But yes it’s true that reducing is faster than splatting. But I think the point here was the clarity of the code not really the performance.
2 Likes
Your code still wastes memory to create the useless vector fill(F, n). Also, it doesn’t support the case n=0. Just do this instead
Fn(x, n) = foldl((y,_) -> F(y), 1:n; init=x)
Edit: of course one can use a StaticArrays.SVector{3} to sqeeze out the last bit of performance.
3 Likes
Just slightly over, but here’s a user-defined linear regression in Soss.jl:
julia> using Soss, MeasureTheory, SampleChainsDynamicHMC
julia> a = -0.2; b=2.0;
julia> x = randn(100); y = a .+ b .* x;
julia> m = @model x begin
a ~ Normal()
b ~ Normal()
y ~ For(x) do xj Normal(a + b * xj, 1) end
end;
julia> sample(DynamicHMCChain, m(x=x) | (;y=y))
4000-element MultiChain with 4 chains and schema (b = Float64, a = Float64)
(b = 1.985±0.097, a = -0.207±0.11)
9 Likes
It’s also standard notation for function iteration in the Dynamical Systems literature.
I agree it’s super neat!
Can you help me understand, what does “unsafe” & “type piracy” mean (in general & in this context)?
Btw, for composition I define a Nest(f,x0,n) function like Mathematica:
Nest(f,x0,n) = (n == 0) ? x0 : (f∘Nest)(f,x0,n-1)
import Base.^;
h^(n::Integer) = x0 -> foldl( (x, _) -> h(x), 1:n; init = x0)
2 Likes
Type piracy is when a module defines a method but neither the function nor any of the method’s parameter types are created by that module.
1 Like
Thanks, but what does that mean in this context?
Base.:^(f::Function, n) = n==1 ? f : f ∘ f^(n-1)
That snippet uses ^ and Function, but both of those come from another module, Base. It would not be type piracy to define
struct CustomCounter end
Base.:^(f::Function, n::CustomCounter) = n==1 ? f : f ∘ f^(n-1)
because CustomCounter is defined in this module.
9 Likes
It’s also super slow and wasteful
julia> @time @btime (F^(10^3))([1,1,2,3,5])
660.172 μs (2001 allocations: 218.89 KiB)
10.475864 seconds (29.82 M allocations: 3.087 GiB, 4.50% gc time, 0.24% compilation time)
3-element Vector{Float64}:
2.089057949736859
2.089057949736859
2.089057949736859
and doomed to fail miserably
julia> (F^(10^5))([1,1,2,3,5])
...
Bonus: if one is really concerned about the type piracy, then
(f::Function)↑(n::Integer) = x0 -> foldl((x, _) -> f(x), 1:n; init = x0)
is a suggestive alternative (↑ is \uparrow).
Edit: in the first code snippet, @time @btime is a typo. The intended code was just @time, to take into consideration also the compilation time, instead of @btime.
4 Likes
Try my version w/ ‘foldl’
Yes your version with foldl is perfectly fine!
1 Like
I like Pink noise. It’s is Gaussian noise with power spectral density 1/fᵅ. One can sample it with the inverse Fourier transform
# Make pink noise
using FFTW, Statistics
function pinknoise(m, n)
A = real(ifft((Complex{Float64}[i+j==2 ? 0. : (randn()+randn()*im)/sqrt((sin((i-1)*pi/m)^2 + sin((j-1)* pi/n)^2)) for i in 1:m, j in 1:n])))
real(A)/std(A)
end
A = pinknoise(500, 500)
It most spectacular property is that it is invariant under conformal mappings, functions that locally preserves angles, but not necessarily lengths, for example the map x + y*im -> ((x + y*im)/33)^2.
# Make conformal map
getindexfc(A, i, j) = A[mod1(ceil(Int, i), size(A,1)), mod1(ceil(Int, j), size(A,2))]
finv(x) = reim(sqrt(x[1] + x[2]*im)*33) # use inverse for look-up
cfmap(A) = map(CartesianIndices(size(A))) do I
getindexfc(A, finv(I)...)
end
This is how it acts on a grid matrix G:
[G cfmap(G)]
And here with pink noise
[A cfmap(A)]
PS: I took a bit extra effort to make it tiling, try [A A A; A A A; A A A]
27 Likes