upload solutions to JuliaIntro

Author: FHoltorf

day_9/JuliaIntro/solutions/Ex1/cos_approx.jl

@@ -0,0 +1,19 @@
+using BenchmarkTools
+
+function cos_approx(x, N)
+    # approximation of cosine via power series expansion
+    # inputs:
+    #       - x : argument of cosine 
+    #       - N : truncation order of the power series approximation
+    # outputs:
+    #       - cos_val : approximation of cos(x)
+    cos_val = 0
+    for n in 0:N
+        cos_val += (-1)^n * x^(2n)/(2factorial(n))
+    end
+    return cos_val
+end
+
+@btime cos_approx($(π/3),$(10)) 
+@btime cos($(π/3))
+@btime cos(π/3)

day_9/JuliaIntro/solutions/Ex1/cos_approx.py

@@ -0,0 +1,24 @@
+import timeit
+import numpy as np
+from math import factorial
+
+
+def cos_approx(x,order):
+    return sum((-1)**n * x**(2*n)/(2*factorial(n)) for n in range(order+1))
+
+a = 3.14/3
+order = 10
+def benchmark_cos_approx():
+    cos_approx(a,order)
+    pass
+
+def benchmark_cos():
+    np.cos(a)
+    pass
+
+execution_time = timeit.Timer(benchmark_cos_approx).timeit(number = 100000)/100000
+print("Execution of cos_approx took " + str(execution_time*1e9) + " ns")
+
+execution_time = timeit.Timer(benchmark_cos).timeit(number = 100000)/100000
+print("Execution of np.cos took " + str(execution_time*1e9) + " ns")
+

day_9/JuliaIntro/solutions/Ex2/structs.jl

@@ -0,0 +1,36 @@
+import Base.+
+
+struct MyComplex
+    Re
+    Im
+end
+
+function my_times(a::MyComplex, b::MyComplex)
+    real_part = a.Re * b.Re - a.Im * b.Im
+    imaginary_part = a.Re * b.Im + a.Im * b.Re
+    return MyComplex(real_part, imaginary_part)
+end
+
+# test-case
+x = MyComplex(1.0, 2.0)
+y = MyComplex(3.0, 4.0)
+z = my_times(x,y)
+
+# compare to Julia's built in complex arithmetic!
+x_julia = 1.0 + 1.0*im  # defines a complex number: 1 + 1i
+y_julia = 0.5 + 3.0*im  # defines a complex number: 0.5 + 3i
+z_julia = x_julia * y_julia
+
+# Implement function my_add that adds
+# two objects of type MyComplex and returns the result as such
+function my_add(a::MyComplex,b::MyComplex)
+    return MyComplex(a.Re+b.Re, a.Im+b.Im)
+end
+
+# overload Julia's + function
+function +(a::MyComplex, b::MyComplex)
+    return my_add(a, b)
+end
+
+z = x + y
+sum([x,y])

day_9/JuliaIntro/solutions/Ex3/Motzkin.jl

@@ -0,0 +1,15 @@
+cd(@__DIR__)
+include("dual_number_arithmetic.jl")
+
+# The way we set things up, we can also evaluate gradients of vector-valued functions!
+# As an example, we consider Motzkin's polynomial:
+Motzkin(x,y) = x^4*y^4 + x^2*y^4 - 3*x^2*y^2 + 1
+Motzkin_gradient(x,y) = [4*x^3*y^4 + 2*x*y^4 - 6*x*y^2;
+                         4*x^4*y^3 + 4*x^2*y^3 - 6*x^2*y]
+
+x = DualNumber(rand(), [1.0, 0])
+y = DualNumber(rand(), [0, 1.0])
+
+# check if the result is correct:
+@show Motzkin(x,y).val == Motzkin(x.val, y.val)
+@show Motzkin(x,y).grad == Motzkin_gradient(x.val,y.val) 

day_9/JuliaIntro/solutions/Ex3/dual_number_arithmetic.jl

@@ -0,0 +1,88 @@
+import Base: +, -, *, ^, /, sin, cos, show
+
+#################################
+##### Defining our DualNumber type
+#################################
+struct DualNumber{Tval,Tgrad}
+    val::Tval
+    grad::Tgrad
+end
+
+#################################
+##### product rule
+#################################
+function *(a::DualNumber, b::DualNumber)
+    return DualNumber(a.val*b.val, b.val*a.grad + a.val*b.grad) # product rule
+end
+
+function *(a::Number, b::DualNumber)
+    return DualNumber(a*b.val, a*b.grad) # complete
+end
+
+function *(a::DualNumber, b::Number)
+    return b*a                          # complete
+end
+
+#################################
+##### chain rule for powers
+#################################
+function ^(a::DualNumber, b::Number)
+    return DualNumber(a.val^b, b*a.val^(b-1)*a.grad)
+end
+
+#################################
+##### quotient rule 
+#################################
+function /(a::DualNumber, b::DualNumber)
+    return DualNumber(a.val/b.val, a.grad/b.val - a.val*b.grad/b.val^2 )
+end
+
+function /(a::DualNumber, b::Number)
+    return DualNumber(a.val/b, a.grad/b) 
+end
+
+#################################
+##### addition
+#################################
+function +(a::DualNumber, b::DualNumber) 
+    return DualNumber(a.val+b.val, a.grad + b.grad) # complete
+end
+
+function +(a::DualNumber, b::Number) 
+    return DualNumber(a.val + b, a.grad) # complete
+end
+
+function +(a::Number, b::DualNumber)
+    return b+a # complete
+end
+
+#################################
+##### Chain Rule for sin/cos
+#################################
+function sin(a::DualNumber)
+    return DualNumber(sin(a.val), cos(a.val)*a.grad)
+end
+
+function cos(a::DualNumber)
+    return DualNumber(cos(a.val), -sin(a.val)*a.grad)
+end
+
+#################################
+##### subtraction of dual numbers
+#################################
+function -(a::DualNumber)
+    return DualNumber(-a.val, -a.grad) 
+end
+
+function -(a::DualNumber, b::T) where T <: Union{DualNumber, Number}
+    return a + (-b)
+end
+
+function -(a::Number, b::DualNumber)
+    return a + (-b)
+end
+
+#################################
+##### We are suckers for neat printouts in the REPL!
+#################################
+show(io::IO, a::DualNumber) = println(io, "$(a.val) + $(a.grad)ϵ")

day_9/JuliaIntro/solutions/Ex3/lotka_volterra.jl

@@ -0,0 +1,120 @@
+cd(@__DIR__)
+include("dual_number_arithmetic.jl")
+
+using CairoMakie
+function lotka_volterra(x,α)
+    # right-hand-side of lotka-volterra system
+    # inputs:
+    #       - x : current state x[1] = population of prey, 
+    #                           x[2] = populaiton of predator 
+    #       - α : hunting efficiency of predator species
+    # returns:
+    #       - dx/dt : rate of change of population sizes
+    
+    dx = [x[1] - x[1]*x[2];
+          -2*x[2] + α*x[1]*x[2]]
+    return dx
+end
+
+# or when messing around with unicode symbols ....
+function lotka_volterra(x,α)
+    # right-hand-side of lotka-volterra system
+    # inputs:
+    #       - x : current state x[1] = population of prey, 
+    #                           x[2] = populaiton of predator 
+    #       - α : hunting efficiency of predator species
+    # returns:
+    #       - dx/dt : rate of change of population sizes
+
+    🐇, 🐺 = x
+    d🐇 = 🐇 - α*🐇*🐺
+    d🐺 =  -🐺 + α*🐇*🐺
+    return [d🐇, d🐺]
+end
+
+function simulate_lotka_volterra(x0,α,h,N)
+    # simulator for lotka volterra system (explicit Euler)
+    # inputs:
+    #       - x0 : initial condition 
+    #       - α  : hunting efficiency of predator
+    #       - h  : step size 
+    #       - N  : number of steps 
+    # outputs:
+    #       - trajectory : Vector of population sizes, entry for each time point visisted 
+    #                      during integration
+    #       - time_points: Vector of time points visisted during integration
+    trajectory = [x0]
+    time_points = range(0, N*h, length=N+1)
+    for i in 1:N
+        x_current = trajectory[end] # current state
+        x_hat = x_current + h*lotka_volterra(x_current,α) # explicit euler prediction
+        push!(trajectory, x_hat)
+    end
+    return trajectory, time_points
+end
+
+# initial condition for the simulation
+x0 = [1.0, 0.25]
+# step size
+h = 0.01
+# time horizon (0.0, N*h)
+N = 1500
+# parameter range for the predator efficiency
+α_range = 1.2:0.01:2.8
+
+# gradient prediction by automatic differentiation
+# 
+# please complete the
+# to that end, recall that 
+#       x0_dual = [x0[1] + ∂x0[1]/∂α * ϵ, x0[2] + ∂x0[2]/∂α * ϵ]
+#       α_dual = [α + ∂α/∂α * ϵ]
+# since we want to compute the sensitivity of the final state wrt. parameter α
+#α = DualNumber(1.5,?)
+#x0_dual = [DualNumber(?,?), DualNumber(?,?)]
+α = DualNumber(1.5, 1.0)
+x0_dual = [DualNumber(x0[1], 0.0), DualNumber(x0[2],0.0)]
+
+
+# propagate dual numbers through simulation:
+dual_trajectory, time_points = simulate_lotka_volterra(x0_dual, α, h, N)
+
+##############
+### No more modifications needed beyond this point 
+##############
+
+# General Lotka-Volterra system oscilaltions
+trajectory, time_points = simulate_lotka_volterra(x0, 0.25, h, N)
+
+fig = Figure(fontsize=18); 
+ax_prey = Axis(fig[1,1], ylabel = "prey population", xlabel = "time")
+ax_predator = Axis(fig[2,1], ylabel = "predator population", xlabel = "time")
+lines!(ax_prey, time_points, [x[1].val for x in dual_trajectory])
+lines!(ax_predator, time_points, [x[2].val for x in dual_trajectory])
+display(fig)
+
+# parameter dependence of the final state x(N*h) on α
+x_final = []
+for α in α_range
+    trajectory, time_points = simulate_lotka_volterra(x0, α, h, N)
+    push!(x_final, trajectory[end])
+end
+
+fig = Figure(fontsize=18); 
+ax_prey = Axis(fig[1,1], ylabel = "prey population at $(N*h) months", xlabel = "α")
+ax_predator = Axis(fig[2,1], ylabel = "predator population at $(N*h) months", xlabel = "α")
+lines!(ax_prey, α_range, [x[1] for x in x_final])
+lines!(ax_predator, α_range, [x[2] for x in x_final])
+scatter!(ax_prey, [α.val], [dual_trajectory[end][1].val], color = :red)
+lines!(ax_prey, [α.val-0.2, α.val+0.2], 
+                [dual_trajectory[end][1].val - 0.2*dual_trajectory[end][1].grad,
+                dual_trajectory[end][1].val + 0.2*dual_trajectory[end][1].grad],
+                color = :red)
+
+scatter!(ax_predator, [α.val], [dual_trajectory[end][2].val], color = :red)
+lines!(ax_predator, [α.val-0.2, α.val+0.2], 
+                    [dual_trajectory[end][2].val - 0.2*dual_trajectory[end][2].grad,
+                    dual_trajectory[end][2].val + 0.2*dual_trajectory[end][2].grad],
+                    color = :red)
+display(fig)
+
+ 

day_9/JuliaIntro/solutions/Ex3/simones_test_function.jl

@@ -0,0 +1,37 @@
+cd(@__DIR__)
+include("dual_number_arithmetic.jl")
+using CairoMakie
+
+# Simone's test function
+func(x) = cos(x)+x*sin(x)
+func_dot(x) = x*cos(x)
+
+# points at which we want to evaluate Simone's test function and its derivative 
+x_range = range(-6, 6, length = 100) 
+
+# evaluate Simone's test function for each x in x_range
+vals = func.(x_range)
+
+# evaluate Simeone's test function derivative for each x in x_range
+grads = func_dot.(x_range)
+
+# Now we initialize the dual numbers for the evaluation of the derivative 
+# of Simone's test function as
+x_dual = [DualNumber(x, 1.0) for x in x_range[1:10:end]]
+
+# evaluate Simone's test function for each x in x_dual and save the result in dual_result
+dual_result = func.(x_dual) 
+
+#############################################################
+# plotting 
+#############################################################
+fig = Figure(fontsize=18);
+axs = [Axis(fig[1,1], xlabel = L"x", ylabel = L"f(x)"),
+       Axis(fig[1,2], xlabel = L"x", ylabel = L"\frac{df}{dx}(x)")]
+
+lines!(axs[1], x_range, func.(x_range), color = :black, linestyle=:dash, linewidth=4,label="Analytical")
+lines!(axs[2], x_range, func_dot.(x_range), color = :black, linestyle=:dash, linewidth=4)
+scatter!(axs[1], [x.val for x in x_dual], [x.val for x in dual_result], color = :red, linewidth=2, label = "AD")
+scatter!(axs[2], [x.val for x in x_dual], [x.grad for x in dual_result], color = :red, linewidth=2)
+axislegend(axs[1])
+display(fig)

day_9/JuliaIntro/solutions/Ex4/__init__.py

@@ -0,0 +1,9 @@
+a = [[1/5],
+     [3/40, 9/40], 
+     [44/45, -56/15, 32/9], 
+     [19372/6561, -25360/2187, 64448/6561, -212/729],  
+     [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656], 
+     [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]]
+b = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]
+c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]
+b_control = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]