Merge pull request #99 from drbergman/ic-ecm-xml

ecm variations

Project.toml

@@ -19,6 +19,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PhysiCellCellCreator = "eac91067-0c6e-4a29-b848-3b4e7cf67017"
+PhysiCellECMCreator = "515ed340-d2f8-4e13-b70e-c9ebab3e992f"
 PhysiCellXMLRules = "bad69253-e2d7-49f9-a0cf-2fd158c7c32b"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
@@ -44,6 +45,7 @@ LinearAlgebra = "1.10"
 MAT = "0.10"
 Parameters = "0.12"
 PhysiCellCellCreator = "0.0.3"
+PhysiCellECMCreator = "0.0.3"
 PhysiCellXMLRules = "0.0.4"
 Pkg = "1.10.0"
 QuasiMonteCarlo = "0.3"

docs/make.jl

@@ -20,13 +20,14 @@ makedocs(;
             "Data directory" => "man/data_directory.md",
             "Known limitations" => "man/known_limitations.md",
             "PhysiCell Studio" => "man/physicell_studio.md",
-            "Sensitivity Analysis" => "man/sensitivity_analysis.md",
+            "Sensitivity analysis" => "man/sensitivity_analysis.md",
         ],
         "Documentation" => map(
             s -> "lib/$(s)",
             sort(readdir(joinpath(@__DIR__, "src/lib")))
         ),
         "Miscellaneous" => Any[
+            "Database upgrades" => "misc/database_upgrades.md",
             "Renaming" => "misc/renaming.md",
         ],
     ],

docs/src/lib/VCTAnalysis.md

@@ -8,5 +8,5 @@ Analyze output from a pcvct project.
 
 ```@autodocs
 Modules = [pcvct]
-Pages = ["population.jl", "substrate.jl"]
+Pages = ["population.jl", "substrate.jl", "motility.jl"]
 ```

docs/src/lib/VCTICCell.md

@@ -4,7 +4,7 @@ CollapsedDocStrings = true
 
 # VCTICCell
 
-Functionality for the pcvct-specific XML-based cell initial conditions.
+Functionality for the using [PhysiCellCellCreator.jl](https://github.com/drbergman/PhysiCellECMCreator.jl) to create and use PhysiCell IC cell XML files.
 
 ```@autodocs
 Modules = [pcvct]

docs/src/lib/VCTICECM.md

@@ -0,0 +1,12 @@
+```@meta
+CollapsedDocStrings = true
+```
+
+# VCTICECM
+
+Functionality for the using [PhysiCellECMCreator.jl](https://github.com/drbergman/PhysiCellECMCreator.jl) to create and use PhysiCell ECM XML files.
+
+```@autodocs
+Modules = [pcvct]
+Pages = ["VCTICECM.jl"]
+```

docs/src/man/database_upgrades.md → docs/src/misc/database_upgrades.md

@@ -4,6 +4,11 @@ Not every release will change the database structure, but when one does in a way
 The warning will link to this page and the function will wait for user input to proceed.
 Changes are listed in reverse chronological order.
 
+## to v0.0.15
+Introduce XML-based ECM initial conditions. This introduces `ic_ecm_variations`.
+Also, introduce Dirichlet initial conditions from file, which introduces the `ic_dc_id` in the database.
+For any simulations in the database before upgrading, both of these will be set to `-1` (i.e., no initial conditions) except if  `ic_ecm_id` is not `-1`, in which case `ic_ecm_variation_id` will be set to `0` (i.e., the default ECM initial conditions which is all the original CSV version can handle).
+
 ## to v0.0.10
 Start tracking the PhysiCell version used in the simulation.
 This introduces the `physicell_versions` table which tracks the PhysiCell versions used in simulations.
@@ -15,7 +20,8 @@ Key changes include:
   - If `PhysiCell` is not a git-tracked repo, it will read the `VERSION.txt` file and store that as the `commit_hash` with `-download` appended to the version.
 
 ## to v0.0.3
-Introduce XML-based cell initial conditions. This introduces `ic_cell_variations`. Also, standardized the use of `config_variation` in place of `variation`. Key changes include:
+Introduce XML-based cell initial conditions. This introduces `ic_cell_variations`.
+Also, standardized the use of `config_variation` in place of `variation`. Key changes include:
 - Renaming the `variation_id` column in the `simulations` and `monads` tables to `config_variation_id`.
 - Adding the `ic_cell_variation_id` column to the `simulations` and `monads` tables.
 - In `data/inputs/configs`, renaming all instances of "variation" to "config_variation" in filenames and databases.

src/VCTAnalysis.jl

@@ -1,2 +1,3 @@
 include("VCTAnalysis/population.jl")
-include("VCTAnalysis/substrate.jl")
+include("VCTAnalysis/substrate.jl")
+include("VCTAnalysis/motility.jl")

src/VCTAnalysis/motility.jl

@@ -0,0 +1,76 @@
+export getCellPositionSequence, computeMeanSpeed
+
+"""
+    getCellPositionSequence(sequence::PhysiCellSequence; include_dead::Bool=false, include_cell_type::Bool=false)
+
+Return a dictionary where the keys are cell IDs from the PhysiCell simulation and the values are NamedTuples containing the time and the position of the cell.
+
+This is a convenience function for `getCellDataSequence(sequence, "position"; include_dead=include_dead, include_cell_type=include_cell_type)`.
+"""
+function getCellPositionSequence(sequence::PhysiCellSequence; include_dead::Bool=false, include_cell_type::Bool=false)
+    return getCellDataSequence(sequence, "position"; include_dead=include_dead, include_cell_type=include_cell_type)
+end
+
+function meanSpeed(p; direction=:any)::NTuple{3,Dict{String,Float64}}
+    x, y, z = [col for col in eachcol(p.position)]
+    cell_type_name = p.cell_type_name
+    dx = x[2:end] .- x[1:end-1]
+    dy = y[2:end] .- y[1:end-1]
+    dz = z[2:end] .- z[1:end-1]
+    if direction == :x
+        dist_fn = (dx, dy, dz) -> abs(dx)
+    elseif direction == :y
+        dist_fn = (dx, dy, dz) -> abs(dy)
+    elseif direction == :z
+        dist_fn = (dx, dy, dz) -> abs(dz)
+    elseif direction == :any
+        dist_fn = (dx, dy, dz) -> sqrt(dx ^ 2 + dy ^ 2 + dz ^ 2)
+    else
+        error("Invalid direction: $direction")
+    end
+    type_change = cell_type_name[2:end] .!= cell_type_name[1:end-1]
+    start_ind = 1
+    cell_type_names = unique(cell_type_name)
+    distance_dict = Dict{String, Float64}(zip(cell_type_names, zeros(Float64, length(cell_type_names))))
+    time_dict = Dict{String, Float64}(zip(cell_type_names, zeros(Float64, length(cell_type_names))))
+    while start_ind <= length(type_change) 
+        I = findfirst(type_change[start_ind:end]) # from s to I, cell_type_name is constant. at I+1 it changes
+        I = isnothing(I) ? length(type_change)+2-start_ind : I # if the cell_type_name is constant till the end, set I to be at the end
+        # If start_ind = 1 (at start of sim) and I = 2 (so cell_type_name[3] != cell_type_name[2], meaning that for steps [1,2] cell_type_name is constnat), only use dx in stepping from 1->2 since somewhere in 2->3 the type changes. That is, use dx[1]
+        distance_dict[cell_type_name[start_ind]] += sum(dist_fn.(dx[start_ind:I-1], dy[start_ind:I-1], dz[start_ind:I-1]))  # only count distance travelled while remaining in the initial cell_type_name
+        time_dict[cell_type_name[start_ind]] += p.time[start_ind+I-1] - p.time[start_ind] # record time spent in this cell_type_name (note p.time is not diffs like dx and dy are, hence the difference in indices)
+        start_ind += I # advance the start to the first instance of a new cell_type_name
+    end
+    speed_dict = [k => distance_dict[k] / time_dict[k] for k in cell_type_names] |> Dict{String,Float64} # convert to speed
+    return speed_dict, distance_dict, time_dict
+end
+
+function computeMeanSpeed(folder::String; direction=:any)::NTuple{3,Vector{Dict{String,Float64}}}
+    sequence = PhysiCellSequence(folder; include_cells=true)
+    pos = getCellPositionSequence(sequence; include_dead=false, include_cell_type=true)
+    dicts = [meanSpeed(p; direction=direction) for p in values(pos) if length(p.time) > 1]
+    return [dict[1] for dict in dicts], [dict[2] for dict in dicts], [dict[3] for dict in dicts]
+end
+
+"""
+    computeMeanSpeed(simulation_id::Integer[; direction=:any])
+
+Return dictionaries containing the mean speed, total distance traveled, and total time spent for each cell type in the PhysiCell simulation.
+
+The time is counted from when the cell first appears in simulation output until it dies or the simulation ends, whichever comes first.
+
+To account for cells that may change cell type during the simulation, the dictionaries returned are keyed by cell type.
+So, a dictionary with key "A" and value 2.0 indicates that the mean speed of this cell while it was of type "A" is 2.0.
+
+# Arguments
+- `simulation_id::Integer`: The ID of the PhysiCell simulation.
+- `direction::Symbol`: The direction to compute the mean speed. Can be `:x`, `:y`, `:z`, or `:any` (default). If `:x`, for example, the mean speed is calculated using only the x component of the cell's movement.
+
+# Returns
+- `mean_speed_dicts::Vector{Dict{String,Float64}}`: A vector of dictionaries where each dictionary is specific to a single cell. The key is the cell type and the value is the mean speed of that cell.
+- `distance_dicts::Vector{Dict{String,Float64}}`: A vector of dictionaries where each dictionary is specific to a single cell. The key is the cell type and the value is the total distance traveled by that cell.
+- `time_dicts::Vector{Dict{String,Float64}}`: A vector of dictionaries where each dictionary is specific to a single cell. The key is the cell type and the value is the total time in the simulation for that cell.
+"""
+function computeMeanSpeed(simulation_id::Integer; direction=:any)::NTuple{3,Vector{Dict{String,Float64}}}
+    return joinpath(outputFolder("simulation", simulation_id), "output") |> x -> computeMeanSpeed(x; direction=direction)
+end

src/VCTClasses.jl

@@ -91,6 +91,7 @@ struct VariationIDs
     config::Int # integer identifying which variation on the base config file to use (config_variations.db)
     rulesets_collection::Int # integer identifying which variation on the ruleset file to use (rulesets_collection_variations.db)
     ic_cell::Int # integer identifying which variation on the ic cell file to use (ic_cell_variations.db) (only used if cells.xml, not used for cells.csv)
+    ic_ecm::Int # integer identifying which variation on the ic ecm file to use (ic_ecm_variations.db) (only used if ecm.xml, not used for ecm.csv)
 end
 
 function VariationIDs(inputs::InputFolders)
@@ -139,6 +140,7 @@ struct Simulation <: AbstractMonad
         @assert variation_ids.config >= 0 "config variation id must be non-negative"
         @assert variation_ids.rulesets_collection >= -1 "rulesets_collection variation id must be non-negative or -1 (indicating no rules)"
         @assert variation_ids.ic_cell >= -1 "ic_cell variation id must be non-negative or -1 (indicating no ic cells)"
+        @assert variation_ids.ic_ecm >= -1 "ic_ecm variation id must be non-negative or -1 (indicating no ic ecm)"
         if variation_ids.rulesets_collection != -1
             @assert inputs.rulesets_collection.folder != "" "rulesets_collection folder must be provided if rulesets_collection variation id is not -1 (indicating that the rules are in use)"
         end
@@ -148,6 +150,12 @@ struct Simulation <: AbstractMonad
             @assert inputs.ic_cell.folder != "" "ic_cell folder must be provided if ic_cell variation_id is not -1 (indicating that the cells are in use)"
             @assert variation_ids.ic_cell == 0 || isfile(joinpath(data_dir, "inputs", "ics", "cells", inputs.ic_cell.folder, "cells.xml")) "cells.xml must be provided if ic_cell variation_id is >1 (indicating that the cell ic parameters are being varied)"
         end
+        if variation_ids.ic_ecm == -1
+            @assert inputs.ic_ecm.folder == "" "ic_ecm variation_id must be >=0 if ic_ecm folder is provided"
+        else
+            @assert inputs.ic_ecm.folder != "" "ic_ecm folder must be provided if ic_ecm variation_id is not -1 (indicating that the ecm is in use)"
+            @assert variation_ids.ic_ecm == 0 || isfile(joinpath(data_dir, "inputs", "ics", "ecms", inputs.ic_ecm.folder, "ecm.xml")) "ecm.xml must be provided if ic_ecm variation_id is >1 (indicating that the ecm ic parameters are being varied)"
+        end
         return new(id, inputs, variation_ids)
     end
 end
@@ -181,7 +189,7 @@ function getSimulation(simulation_id::Int)
         error("Simulation $(simulation_id) not in the database.")
     end
     inputs = InputFolders(df.config_id[1], df.custom_code_id[1], df.rulesets_collection_id[1], df.ic_cell_id[1], df.ic_substrate_id[1], df.ic_ecm_id[1], df.ic_dc_id[1])
-    variation_ids = VariationIDs(df.config_variation_id[1], df.rulesets_collection_variation_id[1], df.ic_cell_variation_id[1])
+    variation_ids = VariationIDs(df.config_variation_id[1], df.rulesets_collection_variation_id[1], df.ic_cell_variation_id[1], df.ic_ecm_variation_id[1])
     return Simulation(simulation_id, inputs, variation_ids)
 end
 
@@ -313,7 +321,7 @@ function getMonad(monad_id::Int, n_replicates::Int)
         error("Monad $(monad_id) not in the database.")
     end
     inputs = InputFolders(df.config_id[1], df.custom_code_id[1], df.rulesets_collection_id[1], df.ic_cell_id[1], df.ic_substrate_id[1], df.ic_ecm_id[1], df.ic_dc_id[1])
-    variation_ids = VariationIDs(df.config_variation_id[1], df.rulesets_collection_variation_id[1], df.ic_cell_variation_id[1])
+    variation_ids = VariationIDs(df.config_variation_id[1], df.rulesets_collection_variation_id[1], df.ic_cell_variation_id[1], df.ic_ecm_variation_id[1])
     use_previous = true
     return Monad(monad_id, n_replicates, inputs, variation_ids, use_previous)
 end
@@ -455,21 +463,24 @@ function Sampling(inputs::InputFolders;
                 config_variation_ids::Union{Int,AbstractArray{<:Integer}}=Int[], 
                 rulesets_collection_variation_ids::Union{Int,AbstractArray{<:Integer}}=fill(inputs.rulesets_collection.folder=="" ? -1 : 0, size(config_variation_ids)),
                 ic_cell_variation_ids::Union{Int,AbstractArray{<:Integer}}=fill(inputs.ic_cell.folder=="" ? -1 : 0, size(config_variation_ids)),
+                ic_ecm_variation_ids::Union{Int,AbstractArray{<:Integer}}=fill(inputs.ic_ecm.folder=="" ? -1 : 0, size(config_variation_ids)),
                 use_previous::Bool=true) 
     # allow for passing in a single config_variation_id and/or rulesets_collection_variation_id
     # later, can support passing in (for example) a 3x6 config_variation_ids and a 3x1 rulesets_collection_variation_ids and expanding the rulesets_collection_variation_ids to 3x6, but that can get tricky fast
-    if all(x->x isa Integer, [config_variation_ids, rulesets_collection_variation_ids, ic_cell_variation_ids])
+    if all(x->x isa Integer, [config_variation_ids, rulesets_collection_variation_ids, ic_cell_variation_ids, ic_ecm_variation_ids])
         config_variation_ids = [config_variation_ids]
         rulesets_collection_variation_ids = [rulesets_collection_variation_ids]
         ic_cell_variation_ids = [ic_cell_variation_ids]
+        ic_ecm_variation_ids = [ic_ecm_variation_ids]
     else
-        ns = [length(x) for x in [config_variation_ids, rulesets_collection_variation_ids, ic_cell_variation_ids] if !(x isa Integer)]
-        @assert all(x->x==ns[1], ns) "config_variation_ids, rulesets_collection_variation_ids, and ic_cell_variation_ids must have the same length if they are not integers"
+        ns = [length(x) for x in [config_variation_ids, rulesets_collection_variation_ids, ic_cell_variation_ids, ic_ecm_variation_ids] if !(x isa Integer)]
+        @assert all(x->x==ns[1], ns) "config_variation_ids, rulesets_collection_variation_ids, ic_cell_variation_ids, ic_ecm_variation_ids must have the same length if they are not integers"
         config_variation_ids = config_variation_ids isa Integer ? fill(config_variation_ids, ns[1]) : config_variation_ids
         rulesets_collection_variation_ids = rulesets_collection_variation_ids isa Integer ? fill(rulesets_collection_variation_ids, ns[1]) : rulesets_collection_variation_ids
         ic_cell_variation_ids = ic_cell_variation_ids isa Integer ? fill(ic_cell_variation_ids, ns[1]) : ic_cell_variation_ids
+        ic_ecm_variation_ids = ic_ecm_variation_ids isa Integer ? fill(ic_ecm_variation_ids, ns[1]) : ic_ecm_variation_ids
     end
-    variation_ids = [VariationIDs(config_variation_ids[i], rulesets_collection_variation_ids[i], ic_cell_variation_ids[i]) for i in 1:length(config_variation_ids)]
+    variation_ids = [VariationIDs(config_variation_ids[i], rulesets_collection_variation_ids[i], ic_cell_variation_ids[i], ic_ecm_variation_ids[i]) for i in 1:length(config_variation_ids)]
     return Sampling(n_replicates, inputs, variation_ids; use_previous=use_previous) 
 end
 
@@ -505,7 +516,7 @@ function getSampling(sampling_id::Int, n_replicates::Int)
     monad_ids = readSamplingMonadIDs(sampling_id)
     inputs = InputFolders(df.config_id[1], df.custom_code_id[1], df.rulesets_collection_id[1], df.ic_cell_id[1], df.ic_substrate_id[1], df.ic_ecm_id[1], df.ic_dc_id[1])
     monad_df = constructSelectQuery("monads", "WHERE monad_id IN ($(join(monad_ids,",")))") |> queryToDataFrame
-    variation_ids = [VariationIDs(monad_df.config_variation_id[i], monad_df.rulesets_collection_variation_id[i], monad_df.ic_cell_variation_id[i]) for i in 1:length(monad_ids)]
+    variation_ids = [VariationIDs(monad_df.config_variation_id[i], monad_df.rulesets_collection_variation_id[i], monad_df.ic_cell_variation_id[i], monad_df.ic_ecm_variation_id[i]) for i in 1:length(monad_ids)]
     return Sampling(sampling_id, n_replicates, monad_ids, inputs, variation_ids)
 end
 

src/VCTCompilation.jl

@@ -40,9 +40,6 @@ function loadCustomCode(S::AbstractSampling; force_recompile::Bool=false)
 
     executable_name = baseToExecutable("project_ccid_$(S.inputs.custom_code.id)")
     cmd = `make -j 8 CC=$(PHYSICELL_CPP) PROGRAM_NAME=$(executable_name) CFLAGS=$(cflags)`
-    # if Sys.isapple() # hacky way to say the -j flag works on my machine but not on the HPC
-    #     cmd = `$cmd -j 20`
-    # end
 
     println("Compiling custom code for $(S.inputs.custom_code.folder) with flags: $cflags")