BerkeleyLab · rouson · Jun 28, 2022 · Jun 22, 2022 · Jun 22, 2022 · Jun 22, 2022
diff --git a/app/main.f90 b/app/main.f90
diff --git a/app/matcha.f90 b/app/matcha.f90
@@ -0,0 +1,42 @@
+program matcha
+  !! Motility Analysis of T-Cell Histories in Activation (Matcha)
+  use matcha_m, only : t_cell_collection_t, distribution_t
+  implicit none
+
+  integer, parameter :: ncells = 100, npositions = 25, ndim = 3, nintervals = 10
+  double precision, parameter :: scale = 100.D0, dt = 0.1D0
+  double precision, allocatable :: random_positions(:,:)
+  type(t_cell_collection_t), allocatable :: history(:)
+
+  allocate(random_positions(ncells,ndim))
+  call random_number(random_positions)  
+  history = [t_cell_collection_t(scale*random_positions, time=0.D0)]
+
+  block
+    integer, parameter :: nveldim = 4, nsteps = npositions - 1
+    double precision, allocatable :: random_4vectors(:,:,:), sample_distribution(:), velocities(:,:,:)
+    type(distribution_t) distribution
+    integer step
+
+    allocate(random_4vectors(ncells,nsteps,nveldim))
+    call random_number(random_4vectors)  
+    allocate(sample_distribution(nintervals))
+    call random_number(sample_distribution)
+    sample_distribution = sample_distribution/sum(sample_distribution)
+    distribution = distribution_t(sample_distribution)
+
+    associate(random_speeds => random_4vectors(:,:,1), random_directions => random_4vectors(:,:,2:4))
+      associate(v => distribution%velocities(random_speeds, random_directions))
+        do step = 1, nsteps
+          associate(x => history(step)%positions(), t => history(step)%time())
+            history = [history, t_cell_collection_t(x + v(:,step,:)*dt, t + dt)]
+          end associate
+        end do
+      end associate
+    end associate
+  end block
+
+  print *
+  print *,"----> Matcha done. <----"
+
+end program
diff --git a/doc/class-diagram.puml b/doc/class-diagram.puml
@@ -0,0 +1,21 @@
+@startuml
+
+Title "Matcha Classes"
+
+class t_cell_collection_t{
+  - positions_ : double precision[:,:]
+  - time : double precision
+  + positions() : real[:,:]
+  + t_cell_collection_t(positions : double precision[:,:], time : double precision) : t_cell_collection_t
+}
+
+class distribution_t{
+  - vel_ : double precision[:]
+  - cumulative_distribution_ : double precision[:]
+  + distribution_t(sample_distribution : double precision[:]) : distribution_t
+}
+note right of distribution_t::distribution_t
+  context distribution_t::distribution_t() post: all([(cumulative_distribution_(i+1) >= cumulative_distribution_(i), i=1,size(cumulative_distribution_,1)-1)]
+end note
+
+@enduml
diff --git a/src/matcha/distribution_m.f90 b/src/matcha/distribution_m.f90
@@ -8,13 +8,13 @@ module distribution_m
     private
     double precision, allocatable, dimension(:) :: vel_, cumulative_distribution_
   contains
-    procedure :: vel
     procedure :: cumulative_distribution
+    procedure :: velocities
   end type  
 
   interface distribution_t
 
-   module function construct(sample_distribution) result(distribution)
+ pure module function construct(sample_distribution) result(distribution)
       implicit none
       double precision, intent(in) :: sample_distribution(:)
       type(distribution_t) distribution
@@ -24,17 +24,20 @@ module function construct(sample_distribution) result(distribution)
 
   interface
 
-    pure module function vel(self) result(my_vel)
+    pure module function cumulative_distribution(self) result(my_cumulative_distribution)
       implicit none
       class(distribution_t), intent(in) :: self
-      double precision, allocatable :: my_vel(:)
+      double precision, allocatable :: my_cumulative_distribution(:)
     end function
 
-    pure module function cumulative_distribution(self) result(my_cumulative_distribution)
+    pure module function velocities(self, speeds, directions) result(my_velocities)
+      !! Return the t_cell_collection_t object's velocity vectors
       implicit none
       class(distribution_t), intent(in) :: self
-      double precision, allocatable :: my_cumulative_distribution(:)
+      double precision, intent(in) :: speeds(:,:), directions(:,:,:)
+      double precision, allocatable :: my_velocities(:,:,:)
     end function
+
   end interface
 
 end module distribution_m
diff --git a/src/matcha/distribution_s.f90 b/src/matcha/distribution_s.f90
@@ -7,21 +7,60 @@
 
   module procedure construct
     integer i
-
-    call assert(all(sample_distribution>=0.D0), "distribution_t: sample_distribution>=0.", intrinsic_array_t(sample_distribution))
-
+
+    call assert(all(sample_distribution>=0.D0), "distribution_t%construct: sample_distribution>=0.", &
+      intrinsic_array_t(sample_distribution))
+
     associate(nintervals => size(sample_distribution,1))      
       distribution%vel_ = [(dble(i), i =1, nintervals)]  ! Assign speeds to each distribution bin         
       distribution%cumulative_distribution_ = [0.D0, [(sum(sample_distribution(1:i)), i=1, nintervals)]]
+
+      call assert(all([(distribution%cumulative_distribution_(i+1) >= distribution%cumulative_distribution_(i), i=1,nintervals)]),&
+        "distribution_t: cumulative_distribution increases monotonically", intrinsic_array_t(sample_distribution))
     end associate
+
   end procedure construct
-
-  module procedure vel
-    my_vel = self%vel_
-  end procedure
-
+
   module procedure cumulative_distribution
+    call assert(allocated(self%cumulative_distribution_), &
+      "distribution_t%cumulative_distribution: allocated(cumulative_distribution_)")
     my_cumulative_distribution = self%cumulative_distribution_
   end procedure 
-
+
+  module procedure velocities
+
+    double precision, allocatable :: sampled_speeds(:,:),  dir(:,:,:)
+    integer cell, step
+
+    ! Sample from the distribution
+    associate(ncells => size(speeds,1), nsteps => size(speeds,2))
+      allocate(sampled_speeds(ncells,nsteps))
+      do concurrent(cell = 1:ncells, step = 1:nsteps)
+        associate(k => findloc(speeds(cell,step) >= self%cumulative_distribution(), value=.true., dim=1))
+          sampled_speeds(cell,step) = self%vel_(k)
+        end associate
+      end do
+
+      ! Create unit vectors
+      dir = directions(:,1:nsteps,:)
+
+      associate(dir_mag => sqrt(dir(:,:,1)**2 +dir(:,:,2)**2 + dir(:,:,3)**2))
+        associate(dir_mag_ => merge(dir_mag, epsilon(dir_mag), dir_mag/=0.))
+          dir(:,:,1) = dir(:,:,1)/dir_mag_
+          dir(:,:,2) = dir(:,:,2)/dir_mag_
+          dir(:,:,3) = dir(:,:,3)/dir_mag_
+        end associate
+      end associate
+
+      allocate(my_velocities, mold=dir)
+
+      do concurrent(step=1:nsteps)
+        my_velocities(:,step,1) = sampled_speeds(:,step)*dir(:,step,1)
+        my_velocities(:,step,2) = sampled_speeds(:,step)*dir(:,step,2)
+        my_velocities(:,step,3) = sampled_speeds(:,step)*dir(:,step,3)
+      end do
+    end associate
+
+  end procedure velocities
+
 end submodule distribution_s
diff --git a/src/matcha/move_m.f90 b/src/matcha/move_m.f90
diff --git a/src/matcha/move_s.f90 b/src/matcha/move_s.f90
diff --git a/src/matcha/t_cell_collection_m.f90 b/src/matcha/t_cell_collection_m.f90
@@ -1,5 +1,6 @@
 module t_cell_collection_m
   !! Define a T-cell abstraction for motility simulations
+  use distribution_m, only : distribution_t
   implicit none
 
   private
@@ -9,17 +10,18 @@ module t_cell_collection_m
     !! Encapsulate the state of a collection of T cells
     private
     double precision, allocatable :: positions_(:,:) !! position vectors
-    double precision time_ !! time stample
+    double precision time_ !! time stamp
   contains
     procedure :: positions
+    procedure :: time
   end type
 
   interface t_cell_collection_t
 
-    pure module function construct(positions, scale, time) result(t_cell_collection)
+    pure module function construct(positions, time) result(t_cell_collection)
       !! Return a t_cell_collection_t object rescaled position vectors and the provided time stamp
       implicit none
-      double precision, intent(in) :: positions(:,:), scale, time
+      double precision, intent(in) :: positions(:,:), time
       type(t_cell_collection_t) t_cell_collection
     end function 
 
@@ -34,6 +36,14 @@ pure module function positions(self) result(my_positions)
       double precision, allocatable :: my_positions(:,:)
     end function
 
+
+    pure module function time(self) result(my_time)
+      !! Return the t_cell_collection_t object's time stamp
+      implicit none
+      class(t_cell_collection_t), intent(in) :: self
+      double precision my_time
+    end function
+
   end interface 
 
 end module t_cell_collection_m
diff --git a/src/matcha/t_cell_collection_s.f90 b/src/matcha/t_cell_collection_s.f90
@@ -4,12 +4,16 @@
 contains
 
   module procedure construct
-    t_cell_collection%positions_ = positions*scale
+    t_cell_collection%positions_ = positions
     t_cell_collection%time_ = time
   end procedure 
 
   module procedure positions
     my_positions = self%positions_
   end procedure
+
+  module procedure time
+    my_time = self%time_
+  end procedure
 
 end submodule t_cell_collection_s
diff --git a/src/matcha/t_cell_m.f90 b/src/matcha/t_cell_m.f90
diff --git a/src/matcha/t_cell_s.f90 b/src/matcha/t_cell_s.f90
diff --git a/src/matcha_m.f90 b/src/matcha_m.f90
@@ -1,7 +1,5 @@
 module matcha_m 
-  use t_cell_m, only : initialize_positions
+  use t_cell_collection_m, only : t_cell_collection_t
   use distribution_m, only : distribution_t
-  use move_m, only : move_tcells
   implicit none
-
 end module