Browse Source

some fixes, allow negative k in recomputation

Gaspard Jankowiak 6 months ago
parent
commit
c1974c4160
6 changed files with 70 additions and 23 deletions
  1. 4
    1
      configs_metrics/nucleus/test_variamu.yaml
  2. 18
    2
      src/CellSim.jl
  3. 5
    1
      src/CellSimCommon.jl
  4. 2
    2
      src/Cortex.jl
  5. 30
    3
      src/Nucleus.jl
  6. 11
    14
      src/Plotting.jl

+ 4
- 1
configs_metrics/nucleus/test_variamu.yaml View File

@@ -73,12 +73,15 @@ flags:
73 73
     debug:                false
74 74
     force_cortex_area:    true
75 75
     write_metrics:        true
76
+    explicit_mt_force:    false
76 77
 
77 78
 config:
78 79
     output_prefix: test_variamu60
79 80
 
80 81
     plot_period: 20
81 82
 
83
+    recompute_nucleus_each: 5
84
+
82 85
     metrics:
83 86
         start_iteration: 1
84 87
         periods: 2
@@ -94,4 +97,4 @@ config:
94 97
 
95 98
         # shift_nucleus: 0.25 # 3/4 of a period = 3/4 * 2π / f.ω0
96 99
 
97
-   
100
+   

+ 18
- 2
src/CellSim.jl View File

@@ -25,6 +25,10 @@ import NLsolve
25 25
 import YAML
26 26
 import Dates
27 27
 
28
+# we need python's YAML to dump config files
29
+import PyCall
30
+const pyyaml = PyCall.pyimport("yaml")
31
+
28 32
 import DelimitedFiles: readdlm, writedlm
29 33
 
30 34
 import SparseArrays
@@ -223,7 +227,10 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
223 227
     mkpath("runs")
224 228
     mkpath("dumps")
225 229
     mkpath("$(config["output_prefix"])")
226
-    run(`cp $(config["config_filename"]) $(config["output_prefix"])Run_$(config["date_string"]).yaml`)
230
+    config_dict = Dict([("params", CSC.to_dict(P)), ("flags", CSC.to_dict(F)), ("config", config)])
231
+    output_config_file = open(config["output_prefix"] * "Run_" * config["date_string"] * ".yaml", "w")
232
+    pyyaml.dump(config_dict, output_config_file, allow_unicode=true)
233
+    close(output_config_file)
227 234
     println("Parameters:")
228 235
     for s in fieldnames(typeof(P))
229 236
         println(string("    ", s, ": ", getfield(P, s)))
@@ -367,6 +374,7 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
367 374
 
368 375
 
369 376
     stepping = F.DEBUG
377
+    breakpoint = -1
370 378
     if !stepping
371 379
         input_task = @async read(stdin, Char)
372 380
     end
@@ -383,12 +391,20 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
383 391
             break
384 392
         end
385 393
 
394
+        if k == breakpoint
395
+            stepping = true
396
+        end
397
+
386 398
         if stepping
387
-            println("debug: 'q' to quit, 'c' to run continuously, any other key to step, 'b' to run step by step")
399
+            println("debug: 'q' to quit, 'c' to run continuously, 'C' to continue to ITER, any other key to step, 'b' to run step by step")
388 400
             key = read(stdin, 1)[1]
389 401
             if key == 0x63 # c
390 402
                 stepping = false
391 403
                 input_task = @async read(stdin, Char)
404
+            elseif key == 0x43 # C
405
+                stepping = false
406
+                breakpoint = parse(Int, readline())
407
+                input_task = @async read(stdin, Char)
392 408
             elseif key == 0x71 # q
393 409
                 break
394 410
             end

+ 5
- 1
src/CellSimCommon.jl View File

@@ -124,6 +124,10 @@ mutable struct Flags
124 124
     DEBUG::Bool
125 125
 end
126 126
 
127
+function to_dict(s::Union{Params,Flags})
128
+    return Dict([(String(field), getfield(s, field)) for field in fieldnames(typeof(s))])
129
+end
130
+
127 131
 struct InteractionPotentials
128 132
     N_W::Vector{Float64}
129 133
     N_∇W::Matrix{Float64}
@@ -291,7 +295,7 @@ end
291 295
 """
292 296
 function pointwise_dot_prod(v)
293 297
     N = size(v, 1)
294
-    return [v[:, 1].*SA.sparse(SA.I, N, N) v[:,2].*SA.sparse(SA.I, N, N)]
298
+    return [v[:, 1].*SA.sparse(1.0*SA.I, N, N) v[:,2].*SA.sparse(1.0*SA.I, N, N)]
295 299
 end
296 300
 
297 301
 """

+ 2
- 2
src/Cortex.jl View File

@@ -204,7 +204,7 @@ function init_FD_matrices(P::Params)
204 204
     N = P.N
205 205
     Δσ = P.Δσ
206 206
 
207
-    global M_perp = SA.blockdiag(-SA.sparse(SA.I, N, N), SA.sparse(SA.I, N, N))
207
+    global M_perp = SA.blockdiag(-SA.sparse(1.0*SA.I, N, N), SA.sparse(1.0*SA.I, N, N))
208 208
     global M_cs_plus = SA.spdiagm(-N+1 => [1; zeros(N-1); 1],
209 209
                                    1 => [ones(N-1); 0; ones(N-1)])
210 210
     global M_cs_minus = SA.spdiagm(N-1 => [1; zeros(N-1); 1],
@@ -716,7 +716,7 @@ function compute_residuals_J(x::Vector{Float64},
716 716
                                 differentials, dst_Df, true)
717 717
     end
718 718
 
719
-    dst_Df[:] = SA.sparse(SA.I, 2P.N, 2P.N) .- P.δt*dst_Df
719
+    dst_Df[:] = SA.sparse(1.0*SA.I, 2P.N, 2P.N) .- P.δt*dst_Df
720 720
 end
721 721
 
722 722
 function wrap_residuals(coords::PointCoords, coords_s::PointCoordsShifted,

+ 30
- 3
src/Nucleus.jl View File

@@ -72,7 +72,18 @@ function g_p(x::Vector{Float64}, α::Float64)
72 72
     return -(2α*min.(α*x.-1, 0.0).*log.(α*x) .+ min.(α*x.-1, 0.0).^2 ./x)
73 73
 end
74 74
 
75
-function compute_curvature(c::NucleusCoords)
75
+function compute_curvature_3(c::NucleusCoords, sign_k::Vector{Float64})
76
+    # http://cral-labo.univ-lyon1.fr/labo/fc/Ateliers_archives/ateliers_2005-06/cercle_3pts.pdf
77
+    xc = (0.5*(c.Y[c.circ_idx.p1,1].^2 - c.Y[:,1].^2 + c.Y[c.circ_idx.p1,2].^2 - c.Y[:,2].^2) ./ (c.Y[c.circ_idx.p1,2]-c.Y[:,2])
78
+          - 0.5*(c.Y[:,1].^2 - c.Y[c.circ_idx.m1,1].^2 + c.Y[:,2].^2 - c.Y[c.circ_idx.m1,2].^2) ./ (c.Y[:,2]-c.Y[c.circ_idx.m1,2])) ./ ((c.Y[:,1] - c.Y[c.circ_idx.m1,1])./(c.Y[:,2] - c.Y[c.circ_idx.m1,2]) - (c.Y[c.circ_idx.p1,1] - c.Y[:,1])./(c.Y[c.circ_idx.p1,2] - c.Y[:,2]))
79
+
80
+    yc = -(c.Y[:,1] - c.Y[c.circ_idx.m1,1])./(c.Y[:,2] - c.Y[c.circ_idx.m1,2]) .* xc + 0.5*(c.Y[:,1].^2 - c.Y[c.circ_idx.m1,1].^2 + c.Y[:,2].^2 - c.Y[c.circ_idx.m1,2].^2) ./ (c.Y[:,2]-c.Y[c.circ_idx.m1,2])
81
+    radius = sqrt.((xc - c.Y[:,1]).^2 + (yc - c.Y[:,2]).^2)
82
+    c.k[:] = sign_k ./radius
83
+    c.k[:] = 0.5*(c.k + c.k[c.circ_idx.m1])
84
+end
85
+
86
+function compute_curvature_4(c::NucleusCoords, sign_k::Vector{Float64})
76 87
     A = zeros(3,3)
77 88
     b = zeros(3)
78 89
 
@@ -123,7 +134,7 @@ function compute_curvature(c::NucleusCoords)
123 134
         b[3] = -0.25*(A[1,1] + A[2,2])
124 135
 
125 136
         x = A\b
126
-        c.k[i] = 1/sqrt(x[1]^2 + x[2]^2 - x[3])
137
+        c.k[i] = sign_k[i]./sqrt(x[1]^2 + x[2]^2 - x[3])
127 138
     end
128 139
 end
129 140
 
@@ -136,7 +147,9 @@ function recompute_nucleus_coords(c::NucleusCoords)
136 147
     delta_θ = rem.(c.θ - c.θ[c.circ_idx.m1], 2π, RoundNearest)
137 148
     c.θ[:] = cumsum([c.θ[1]; delta_θ[1:end-1]])
138 149
 
139
-    compute_curvature(c)
150
+    sign_k = [sign(c.θ[1] + 2π - c.θ[end]); sign.(delta_θ)]
151
+
152
+    compute_curvature_4(c, sign_k)
140 153
     c.η[:] = log.(c.r)
141 154
     c.L = sum(c.r)
142 155
 end
@@ -494,10 +507,24 @@ function update_coords(c::NucleusCoords, new_c::NucleusCoords,
494 507
     update_θ(c, new_c, N_W, N_∇W, P, F, temparrays)
495 508
     update_Y(c, new_c, N_W, N_∇W, P, F, temparrays)
496 509
 
510
+    if F.DEBUG & recompute
511
+        println("PRE recompute alpha, beta, r, q, K, θ")
512
+        # display([c.α c.β c.r c.q c.k c.θ])
513
+        display([sum(new_c.α)/P.Nnuc sum(new_c.β)/P.Nnuc sum(new_c.r)/P.Nnuc sum(new_c.q)/P.Nnuc sum(new_c.k)/P.Nnuc sum(new_c.θ)/P.Nnuc])
514
+        println()
515
+        print("kpre = ")
516
+        println(new_c.k)
517
+    end
518
+
497 519
     if recompute
498 520
         recompute_nucleus_coords(new_c)
499 521
     end
500 522
 
523
+    if F.DEBUG & recompute
524
+        print("kpost = ")
525
+        println(new_c.k)
526
+    end
527
+
501 528
     if F.DEBUG
502 529
         println("POST alpha, beta, r, q, K, θ")
503 530
         # display([new_c.α new_c.β new_c.r new_c.q new_c.k new_c.θ])

+ 11
- 14
src/Plotting.jl View File

@@ -22,7 +22,7 @@ function init_plot(coords::Cortex.PointCoords, P::CellSimCommon.Params, F::CellS
22 22
 
23 23
     x = coords.x
24 24
 
25
-    fig = PyPlot.figure(figsize=(12.8, 10))
25
+    fig = PyPlot.figure(figsize=(12.8, 10), dpi=50)
26 26
 
27 27
     PyPlot.show()
28 28
     # figManager = PyPlot.get_current_fig_manager()
@@ -36,7 +36,7 @@ function init_plot(coords::Cortex.PointCoords, P::CellSimCommon.Params, F::CellS
36 36
     ax.plot(x[:,1], x[:,2], ".-", zorder=20)[1]
37 37
 
38 38
     # Cortex fillin
39
-    ax.fill(x[:,1], x[:,2], color="#f713e0", zorder=10)[1]
39
+    ax.fill(x[:,1], x[:,2], color="#f713e033", zorder=10)[1]
40 40
 
41 41
     if F.nucleus
42 42
         # Nucleus
@@ -54,7 +54,9 @@ function init_plot(coords::Cortex.PointCoords, P::CellSimCommon.Params, F::CellS
54 54
     end
55 55
 
56 56
     # Drag force
57
-    # ax.scatter(x[:,1], x[:,2], color="black", zorder=30)
57
+    if F.plot_drag
58
+        ax.scatter(x[:,1], x[:,2], color="black", zorder=30)
59
+    end
58 60
 
59 61
     # Initial condition
60 62
     # ax.plot(x[:,1], x[:,2], color="black", lw=0.5, zorder=1)[1]
@@ -182,8 +184,11 @@ function update_plot(coords::Cortex.PointCoords, nucleus_coords::Union{Nucleus.N
182 184
     end
183 185
 
184 186
     # Drag force
185
-    # scatters[idx_s].set_offsets(x)
186
-    # idx_s += 1
187
+    if F.plot_drag
188
+        scatters[idx_s].set_sizes(10CellSimCommon.@entry_norm(plotables.drag_force))
189
+        scatters[idx_s].set_offsets(x)
190
+        idx_s += 1
191
+    end
187 192
 
188 193
     # Initial condition
189 194
     # idx_l += 1
@@ -230,7 +235,7 @@ function update_plot(coords::Cortex.PointCoords, nucleus_coords::Union{Nucleus.N
230 235
 
231 236
 
232 237
     if F.follow_cam
233
-        if F.follow_nucleus
238
+        if F.nucleus && F.follow_nucleus
234 239
             follow_x = nucleus_coords.Y
235 240
         else
236 241
             follow_x = coords.x
@@ -246,14 +251,6 @@ function update_plot(coords::Cortex.PointCoords, nucleus_coords::Union{Nucleus.N
246 251
         ax.set_ylim((y_mid-0.7*y_span, y_mid+0.7*y_span))
247 252
     end
248 253
 
249
-    # if F.plot_drag
250
-        # scatters[1].set_sizes(10CellSimCommon.@entry_norm(plotables.drag_force))
251
-    # else
252
-        # scatters[1].set_sizes(0CellSimCommon.@entry_norm(plotables.drag_force))
253
-    # end
254
-    # scatters[1].set_facecolor(JankoUtils.scale_cm(plotables.mass_source, PyPlot.get_cmap("RdYlGn");
255
-                                               # range_min=-P.c, range_max=P.c))
256
-
257 254
     PyPlot.draw()
258 255
     sleep(0.001)
259 256
 end

Loading…
Cancel
Save