diff --git a/app/main.f90 b/app/main.f90 deleted file mode 100644 index 618d8d1..0000000 --- a/app/main.f90 +++ /dev/null @@ -1,39 +0,0 @@ -program matcha - use matcha_m, only : initialize_positions, move_tcells - use distribution_m, only : distribution_t - implicit none - - integer, parameter :: ncells = 100, npositions = 25, nveldim = 4, ndim = 3, nintervals = 10 - integer, parameter :: nsteps = npositions - 1 - double precision, allocatable, dimension(:,:) :: x, y, z, random_positions - - allocate(x(ncells,npositions)) - allocate(y(ncells,npositions)) - allocate(z(ncells,npositions)) - allocate(random_positions(ncells,ndim)) - - call random_number(random_positions) - call initialize_positions(random_positions, x(:,1),y(:,1),z(:,1)) - - block - double precision, allocatable :: random_velocities(:,:,:), sample_distribution(:) - type(distribution_t) distribution - - allocate(random_velocities(ncells,nsteps,nveldim)) - allocate(sample_distribution(nintervals)) - - call random_number(sample_distribution) - sample_distribution = sample_distribution/sum(sample_distribution) - - call random_number(random_velocities) - - distribution = distribution_t(sample_distribution) - associate(random_speeds => random_velocities(:,:,1), random_directions => random_velocities(:,:,2:4)) - call move_tcells(distribution, random_speeds, random_directions, x,y,z) - end associate - end block - - print * - print *,"----> Matcha done. <----" - -end program diff --git a/app/matcha.f90 b/app/matcha.f90 new file mode 100644 index 0000000..01c3c06 --- /dev/null +++ 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 new file mode 100644 index 0000000..872bb42 --- /dev/null +++ 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 index 39151bb..1835ed5 100644 --- 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 index 0ada611..d8e4aef 100644 --- 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 deleted file mode 100644 index 4d084f3..0000000 --- a/src/matcha/move_m.f90 +++ /dev/null @@ -1,16 +0,0 @@ -module move_m - use distribution_m, only : distribution_t - implicit none - - interface - - module subroutine move_tcells(distribution, random_speeds, random_directions, x,y,z) - implicit none - type(distribution_t), intent(in) :: distribution - double precision, intent(in) :: random_speeds(:,:), random_directions(:,:,:) - double precision, intent(inout), dimension(:,:) :: x, y, z - end subroutine - - end interface - -end module move_m diff --git a/src/matcha/move_s.f90 b/src/matcha/move_s.f90 deleted file mode 100644 index db89e87..0000000 --- a/src/matcha/move_s.f90 +++ /dev/null @@ -1,54 +0,0 @@ -submodule(move_m) move_s - implicit none - -contains - - module procedure move_tcells - ! Local variables - integer i, j - double precision, allocatable, dimension(:,:) :: speed - - associate(nsteps => size(random_speeds,2)) - associate( & - ncells => size(random_speeds,1), & - cumulative_distribution => distribution%cumulative_distribution(), & - vel => distribution%vel() & - ) - allocate(speed(ncells, nsteps)) - - ! Sample from the distribution - do concurrent(i = 1:ncells, j = 1:nsteps) - associate(k => findloc(random_speeds(i,j) >= cumulative_distribution, value=.true., dim=1)) - speed(i,j) = vel(k) - end associate - end do - end associate - - block - ! Time step - double precision, allocatable, dimension(:,:,:) :: dir - double precision, parameter :: dt = .1D0 - - ! Create a random unit vector - dir = random_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 - - ! Use a forward Euler to advance the cell position - do i=1,nsteps - x(:,i+1) = x(:,i) + dt*speed(:,i)*dir(:,i,1) - y(:,i+1) = y(:,i) + dt*speed(:,i)*dir(:,i,2) - z(:,i+1) = z(:,i) + dt*speed(:,i)*dir(:,i,3) - end do - end block - end associate - - end procedure move_tcells - -end submodule move_s diff --git a/src/matcha/t_cell_collection_m.f90 b/src/matcha/t_cell_collection_m.f90 index 2a59abc..69e1f3f 100644 --- 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 index cc05054..e059b67 100644 --- 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 deleted file mode 100644 index 438d96f..0000000 --- a/src/matcha/t_cell_m.f90 +++ /dev/null @@ -1,14 +0,0 @@ -module t_cell_m - implicit none - - interface - - module subroutine initialize_positions(random_positions,x,y,z) - implicit none - double precision, intent(in) :: random_positions(:,:) - double precision, intent(out) :: x(:), y(:), z(:) - end subroutine - - end interface - -end module t_cell_m diff --git a/src/matcha/t_cell_s.f90 b/src/matcha/t_cell_s.f90 deleted file mode 100644 index d3796db..0000000 --- a/src/matcha/t_cell_s.f90 +++ /dev/null @@ -1,17 +0,0 @@ -submodule(t_cell_m) t_cell_s - implicit none - -contains - - module procedure initialize_positions - - double precision, parameter :: scaling_factor = 100. - - ! Assign initial positions to T cells randomly in a cube with sides of length scaling_factor - x = random_positions(:,1)*scaling_factor - y = random_positions(:,2)*scaling_factor - z = random_positions(:,3)*scaling_factor - - end procedure initialize_positions - -end submodule t_cell_s diff --git a/src/matcha_m.f90 b/src/matcha_m.f90 index c43828b..6903748 100644 --- 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 diff --git a/test/t_cell_collection_test.f90 b/test/t_cell_collection_test.f90 index fb4ec95..7b6f9a7 100644 --- a/test/t_cell_collection_test.f90 +++ b/test/t_cell_collection_test.f90 @@ -1,7 +1,7 @@ module t_cell_collection_test !! summary: unit tests for T-cell collections use garden, only: & - result_t, test_item_t, describe, it, assert_that + result_t, test_item_t, describe, it, assert_that, succeed use t_cell_collection_m, only : t_cell_collection_t implicit none @@ -23,12 +23,17 @@ function check_constructed_domain() result(result_) type(result_t) result_ integer, parameter :: ncells = 100, ndim = 3 double precision random_positions(ncells,ndim) - type(t_cell_collection_t) t_cell_collection double precision, parameter :: scale_factor=100.D0 + type(t_cell_collection_t) t_cell_collection + + call random_number(random_positions) + t_cell_collection = t_cell_collection_t(scale_factor*random_positions,time=0.D0) - call random_number(random_positions) - t_cell_collection = t_cell_collection_t(positions=random_positions, scale=scale_factor, time=0.D0) - result_ = assert_that(all(0.D0 <= t_cell_collection%positions() .and. t_cell_collection%positions() <= scale_factor)) + associate(constructed_positions => t_cell_collection%positions()) + result_ = assert_that( & + all(0.D0 <= constructed_positions .and. constructed_positions <= scale_factor), "position(s) out of range" & + ) + end associate end function end module t_cell_collection_test