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.

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:

28 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