src/fhn_solver.jl

@@ -9,81 +9,41 @@ using Observables
 using ..Constants
 using .Laplacian
 
-"""
+function step_fhn!(U, V, params_obs::Observable; dx=1, dt=0.01)
-    fhn(du, u, p:FHNParams, t:)
-
-    Implements the spatial dynamics of FitzHugh-Nagumo (fhn). Designed to be 
-    within a larger numerical solver of partial differential equations.
-
-    # Arguments
-    - `du`: output argument which stores the calculated derivatives 
-    - `u`: input vector containing the current state of the system at time t
-    - `p`: holds all the fixed parameters of the FHN model
-    - `t`: current time 
-
-    # Returns
-    - `du`: calculated derivatives put back into the du array
-"""
-function fhn!(du, u, p::FHNParams, t)
-    N = p.N
-    u_mat = reshape(u[1:N^2], N, N)
-    v_mat = reshape(u[N^2+1:end], N, N)
-
-    Δu = laplacian(u_mat, p.dx)
-    Δv = laplacian(v_mat, p.dx)
-
-    u_in = u_mat[2:end-1, 2:end-1]
-    v_in = v_mat[2:end-1, 2:end-1]
-
-    fu = p.Du * Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
-    fv = p.Dv * Δv .+ p.ϵ * (u_in .+ p.a .- p.b .* v_in)
-
-    # Construct output with zero boundary (padding)
-    fu_full = zeros(N, N)
-    fv_full = zeros(N, N)
-    fu_full[2:end-1, 2:end-1] .= fu
-    fv_full[2:end-1, 2:end-1] .= fv
-
-    du .= vcat(vec(fu_full), vec(fv_full))
-end
-
-function step_fhn!(U, V, params_obs::Observable; dx=1)
-    """
     p = params_obs[]
 
-    # Flatten initial condition (activation u, recovery v)
+    N = p.N
-    u0 = vec(U)
-    v0 = vec(V)
-    u_init = vcat(u0, v0)
 
-    # Define one integration step using your fhn! function
+    # Compute Laplacians on the interior
-    prob = ODEProblem((du, u, p, t) -> fhn!(du, u, p, t), u_init, (0.0, 0.1), p)
+    Δu = laplacian(U, dx)
-    sol = solve(prob, Tsit5(); dt=0.1, save_start=false, saveat=10.0)
+    Δv = laplacian(V, dx)
 
-    # Extract solution and reshape
+    # Extract interior
-    u_new = reshape(sol[end][1:p.N^2], p.N, p.N)
+    u_in = U[2:end-1, 2:end-1]
-    v_new = reshape(sol[end][p.N^2+1:end], p.N, p.N)
+    v_in = V[2:end-1, 2:end-1]
 
-    # Update matrices in-place
+    # Compute update using FHN reaction-diffusion terms
-    U .= u_new
+    fu = p.Du .* Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
-    V .= v_new
+    fv = p.Dv .* Δv .+ p.ϵ .* (u_in .+ p.a .- p.b .* v_in)
 
-    return U
+    # Euler update
-    """
+    u_new = u_in .+ dt .* fu
+    v_new = v_in .+ dt .* fv
 
-    params = params_obs[]
+    # Write back interior updates
+    U[2:end-1, 2:end-1] .= u_new
+    V[2:end-1, 2:end-1] .= v_new
 
-    Δu = laplacian(U, params.dx)
+    # Apply periodic boundary conditions
-    Δv = laplacian(V, params.dx)
+    U[1, :] .= U[end-1, :]
+    U[end, :] .= U[2, :]
+    U[:, 1] .= U[:, end-1]
+    U[:, end] .= U[:, 2]
 
-    u = U[2:end-1, 2:end-1]
+    V[1, :] .= V[end-1, :]
-    v = V[2:end-1, 2:end-1]
+    V[end, :] .= V[2, :]
-
+    V[:, 1] .= V[:, end-1]
-    fu = params.Du .* Δu .+ u .- (u .^ 3) ./ 3 .- v
+    V[:, end] .= V[:, 2]
-    fv = params.Dv .* Δv .+ params.ϵ .* (u .+ params.a .- params.b .* v)
-
-    U[2:end-1, 2:end-1] .+= 0.1 .* fu
-    V[2:end-1, 2:end-1] .+= 0.1 .* fv
 
     return U
 end

src/utils/laplacian.jl

@@ -11,6 +11,6 @@ function laplacian(U::AbstractMatrix{<:Real}, dx::Real)
     return (up .+ down .+ left .+ right .- 4 .* center) ./ h2
 end
 
-export laplacian, laplacian5
+export laplacian
 
 end # module Laplacian

src/utils/templates.jl

@@ -1,36 +1,8 @@
-
-# initial conditions for different patterns
-function zebra_conditions(N)
-    # Turing-spot parameters
-    params = FHNParams(N=N, dx=1.0, Du=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9)
-
-    # Or use this
-    u0, v0 = two_rows_edge_distance_ic(N)
-
-    return params, vcat(u0, v0)
-end
-
-function cheetah_conditions(N)
-    # Turing-spot parameters
-    #params = FHNParams(N=N, dx=1.0, Du=5e-6, Dv=2e-3, ϵ=0.025, a=0.6, b=0.15)
-    params = FHNParams(N=N, dx=1.0, Du=0.182, Dv=0.5, ϵ=0.05, a=0.2, b=2.0)
-
-    u0 = 0.05 .* (2 .* rand(N, N) .- 1)  # noise in [-0.05, 0.05]
-    v0 = 0.05 .* (2 .* rand(N, N) .- 1)
-
-    return params, vcat(vec(u0), vec(v0))
-end
-
-function coral_conditions(N)
-    # Turing-spot parameters
-    params = FHNParams(N=N, dx=1.0, Du=0.001, Dv=0.06, ϵ=0.05, a=0.0, b=1.2)
-
-    u0, v0 = coral_ic(N)
-
-    return params, vcat(u0, v0)
-end
-
 # helper functions for filling cells in specific places of the matrix
+# Not all functions are used. Those are for experimenting and debugging purposes
+
+module Templates
+
 function blocks_ic(N)
     u = fill(1.0, N, N)
     v = fill(0.0, N, N)
@@ -38,15 +10,15 @@ function blocks_ic(N)
 
     function safe_block!(u, row_center, col_center)
         row_start = max(row_center - 8, 1)
-        row_end   = min(row_center + 7, N)
+        row_end = min(row_center + 7, N)
         col_start = max(col_center - 8, 1)
-        col_end   = min(col_center + 7, N)
+        col_end = min(col_center + 7, N)
         u[row_start:row_end, col_start:col_end] .= -0.01
     end
 
     safe_block!(u, p, p)
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function column_ic(N)
@@ -61,7 +33,18 @@ function column_ic(N)
 
     u[col_start:col_end, :] .= -0.01
 
-    return vec(u), vec(v)
+    return u, v
+end
+
+function stripe_ic(N)
+    u = zeros(N, N)
+    v = zeros(N, N)
+    for i in 1:N
+        for j in 1:N
+            u[i, j] = 0.1 + 0.05 * sin(2π * j / 10) + 0.01 * randn()
+        end
+    end
+    return u, v
 end
 
 function two_rows_edge_distance_ic(N)
@@ -97,7 +80,7 @@ function two_rows_edge_distance_ic(N)
     # Apply the second column
     u[:, col2_start:col2_end] .= -0.01
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function center_band_ic(N)
@@ -111,7 +94,7 @@ function center_band_ic(N)
     u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
     v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function circle_ic(N)
@@ -126,7 +109,7 @@ function circle_ic(N)
         end
     end
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function three_circles_random_ic(N)
@@ -155,7 +138,7 @@ function three_circles_random_ic(N)
         end
     end
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function squiggle_ic(N, Lx=400.0, Ly=400.0)
@@ -180,19 +163,23 @@ function squiggle_ic(N, Lx=400.0, Ly=400.0)
     v = fill(vplus, N, N)
 
     # Apply squiggle
-    u[Z .> 0.8] .= uminus
+    u[Z.>0.8] .= uminus
 
-    return vec(u), vec(v)
+    return u, v
 end
 
 function coral_ic(N)
-    u = fill(0.534522, N, N)
+    u = ones(N, N)
-    v = fill(0.381802, N, N)
+    v = zeros(N, N)
 
-    for _ in 1:40  # place 15 noisy seeds
+    for _ in 1:150  # place 15 noisy seeds
-        i, j = rand(10:N-10), rand(10:N-10)
+        i, j = rand(5:N-5), rand(5:N-5)
         u[i-2:i+2, j-2:j+2] .= -0.534522 .+ 0.2 * rand(5, 5)
     end
 
-    return vec(u), vec(v)
+    return u, v
 end
+
+export blocks_ic, column_ic, squiggle_ic, three_circles_random_ic, circle_ic, center_band_ic, two_rows_edge_distance_ic, coral_ic, stripe_ic
+
+end

src/visualization.jl

@@ -1,11 +1,13 @@
 module Visualization
 
-include("../src/utils/constants.jl")
+include("utils/constants.jl")
+include("utils/templates.jl")
 include("gray_scott_solver.jl")
 include("fhn_solver.jl")
 
 using Observables, Makie, GLMakie
 using .Constants
+using .Templates
 using .GrayScottSolver: step_gray_scott!
 using .FHNSolver: step_fhn!
 
@@ -41,6 +43,7 @@ function reset!(U, V, heat_obs)
     radius = 10
     U[center-radius:center+radius, center-radius:center+radius] .= 0.50
     V[center-radius:center+radius, center-radius:center+radius] .= 0.25
+
     heat_obs[] = copy(U)
 end
 
@@ -71,7 +74,7 @@ end
 function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs)
 
     reset!(U, V, heat_obs)
-    fig = Figure(size=(800, 800))
+    fig = Figure(size=(1300, 950))
 
     gh = GridLayout(fig[1, 1])
     dropdown = Menu(fig, options=collect(zip(["Gray-Scott", "FHN"], [:gray_scott, :fhn])))
@@ -92,32 +95,43 @@ function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs)
     btn_step = Button(buttongrid[1, 1], width=50, label="Step")
     btn_start = Button(buttongrid[1, 2], width=50, label=run_label)
     btn_reset = Button(buttongrid[1, 3], width=50, label="Reset")
-    slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:100, format="{}x", width=350, startvalue=stepsize[]), tellwidth=false)
+    slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:50, format="{}x", width=750, startvalue=stepsize[], tellwidth=false))
 
     speed_slider = slidergrid.sliders[1].value
 
+    gh[1, 2] = templategrid = GridLayout(ax.scene, tellwidth=false)
+    templategrid[1, 1] = Label(fig, "Templates:")
+    btn_waves = Button(templategrid[1, 2], width=100, label="Wave Pattern")
+    btn_cow = Button(templategrid[1, 3], width=100, label="Cow Spots")
+    btn_cheetah = Button(templategrid[1, 4], width=100, label="Cheetah Spots")
+
     # place all the parameter boxes
     gh[2, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false)
 
-    param_box!(textboxgrid, 1, "Du", param_obs_map.Du, col=1)
+    # Create and assign column header labels
-    param_box!(textboxgrid, 2, "Dv", param_obs_map.Dv, col=1)
+    textboxgrid[1, 1] = Label(fig, "GrayScott:", halign=:center)
-    param_box!(textboxgrid, 3, "Feed", param_obs_map.F, col=1)
+    textboxgrid[1, 3] = Label(fig, "FHN:", halign=:center)
-    param_box!(textboxgrid, 4, "Kill", param_obs_map.k, col=1)
+
+    # GrayScott column (col 1)
+    param_box!(textboxgrid, 2, "Du", param_obs_map.Du, col=1)
+    param_box!(textboxgrid, 3, "Dv", param_obs_map.Dv, col=1)
+    param_box!(textboxgrid, 4, "Feed", param_obs_map.F, col=1)
+    param_box!(textboxgrid, 5, "Kill", param_obs_map.k, col=1)
 
     # FHN column (col 2)
-    param_box!(textboxgrid, 1, "Du", param_obs_map.Du, col=2)
+    param_box!(textboxgrid, 2, "Du", param_obs_map.Du, col=2)
-    param_box!(textboxgrid, 2, "Dv", param_obs_map.Dv, col=2)
+    param_box!(textboxgrid, 3, "Dv", param_obs_map.Dv, col=2)
-    param_box!(textboxgrid, 3, "ϵ", param_obs_map.ϵ, col=2)
+    param_box!(textboxgrid, 4, "ϵ", param_obs_map.ϵ, col=2)
-    param_box!(textboxgrid, 4, "a", param_obs_map.a, col=2)
+    param_box!(textboxgrid, 5, "a", param_obs_map.a, col=2)
-    param_box!(textboxgrid, 5, "b", param_obs_map.b, col=2)
+    param_box!(textboxgrid, 6, "b", param_obs_map.b, col=2)
 
-    rowsize!(gh, 1, Fixed(50))  # small row for the menu
+    rowsize!(gh, 1, Relative(0.2))  # small row for the menu
-    rowsize!(gh, 2, Relative(1.0))
+    rowsize!(gh, 2, Relative(0.8))
     for c in 1:4
-        colsize!(textboxgrid, c, Relative(1.0 / 4.0))
+        colsize!(textboxgrid, c, Relative(0.1))
     end
 
-    # Events ##############################################################
+    # ======== Events ========
     # Timer and state for animation
     running = Observable(false)
 
@@ -132,9 +146,9 @@ function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs)
             @warn "Invalid input for $s"
         end
     end
-    # Button Listeners
+    # Control Listeners
     on(btn_step.clicks) do _
-        multi_step!((U, V), stepsize[], heat_obs, params_obs; step_method=step_method[])
+        multi_step!((U, V), stepsize[] * 5, heat_obs, params_obs; step_method=step_method[])
     end
 
     on(btn_start.clicks) do _
@@ -143,16 +157,67 @@ function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs)
 
     on(btn_start.clicks) do _
         @async while running[]
-            multi_step!((U, V), stepsize[] * 5, heat_obs, params_obs; step_method=step_method[])
+            multi_step!((U, V), stepsize[], heat_obs, params_obs; step_method=step_method[])
             sleep(0.0015)  # ~20 FPS
         end
     end
 
     on(btn_reset.clicks) do _
         running[] = false
+        hm.colormap[] = :seismic
         reset!(U, V, heat_obs)
 
     end
+
+    # Template Control
+    on(btn_waves.clicks) do _
+        # add column to center of matrix
+        U, V = Templates.blocks_ic(params_obs[].N)
+
+        hm.colormap[] = :seismic
+
+        # change params
+        param_obs_map.Du[] = 0.0008
+        param_obs_map.Dv[] = 0.1
+        param_obs_map.ϵ[] = 0.01
+        param_obs_map.a[] = 0.5
+        param_obs_map.b[] = 0.15
+        heat_obs[] = copy(U)
+    end
+
+    # Template Control
+    on(btn_cow.clicks) do _
+        # fill matrix with random noise
+        U = 0.05 .* (2 .* rand(params_obs[].N, params_obs[].N) .- 1)  # noise in [-0.05, 0.05]
+        V = 0.05 .* (2 .* rand(params_obs[].N, params_obs[].N) .- 1)
+
+        hm.colormap[] = :grayC
+
+        # change params
+        param_obs_map.Du[] = 1.0
+        param_obs_map.Dv[] = 10.0
+        param_obs_map.ϵ[] = 0.01
+        param_obs_map.a[] = -0.1
+        param_obs_map.b[] = 1.2
+        heat_obs[] = copy(U)
+    end
+
+    on(btn_cheetah.clicks) do _
+        # fill matrix with random noise
+        U = 0.05 .* (2 .* rand(params_obs[].N, params_obs[].N) .- 1)  # noise in [-0.05, 0.05]
+        V = 0.05 .* (2 .* rand(params_obs[].N, params_obs[].N) .- 1)
+
+        hm.colormap[] = cgrad(:lajolla, rev=true)
+
+        # change params
+        param_obs_map.Du[] = 0.182
+        param_obs_map.Dv[] = 0.5
+        param_obs_map.ϵ[] = 0.05
+        param_obs_map.a[] = 0.2
+        param_obs_map.b[] = 2.0
+        heat_obs[] = copy(U)
+    end
+
     on(spoint) do pt
         r = 5
         if pt === nothing