BioFVM_basic_agent.cpp

/*
#############################################################################
# If you use BioFVM in your project, please cite BioFVM and the version     #
# number, such as below:                                                    #
#                                                                           #
# We solved the diffusion equations using BioFVM (Version 1.1.6) [1]        #
#                                                                           #
# [1] A. Ghaffarizadeh, S.H. Friedman, and P. Macklin, BioFVM: an efficient #
#    parallelized diffusive transport solver for 3-D biological simulations,#
#    Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 #
#                                                                           #
#############################################################################
#                                                                           #
# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause)   #
#                                                                           #
# Copyright (c) 2015-2017, Paul Macklin and the BioFVM Project              #
# All rights reserved.                                                      #
#                                                                           #
# Redistribution and use in source and binary forms, with or without        #
# modification, are permitted provided that the following conditions are    #
# met:                                                                      #
#                                                                           #
# 1. Redistributions of source code must retain the above copyright notice, #
# this list of conditions and the following disclaimer.                     #
#                                                                           #
# 2. Redistributions in binary form must reproduce the above copyright      #
# notice, this list of conditions and the following disclaimer in the       #
# documentation and/or other materials provided with the distribution.      #
#                                                                           #
# 3. Neither the name of the copyright holder nor the names of its          #
# contributors may be used to endorse or promote products derived from this #
# software without specific prior written permission.                       #
#                                                                           #
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       #
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED #
# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A           #
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER #
# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,  #
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       #
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        #
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    #
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      #
# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        #
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              #
#                                                                           #
#############################################################################
*/

#include "BioFVM_basic_agent.h"
#include "BioFVM_agent_container.h"
#include "BioFVM_vector.h" 

namespace BioFVM{

std::vector<Basic_Agent*> all_basic_agents(0); 

Basic_Agent::Basic_Agent()
{
	/*-----------------------------------------------------------------------------------------------------*/
    /* give the agent a unique ID - but now there is no need for this as                                   */
    /* the ID is being generated in create_sources() or create_sinks() functions                           */
    /* but if an agent is created directly i.e. Basic_Agent x = new Basic_Agent;                           */
    /* then there will be a problem !                                                                      */
    /* eepecially at non-root processes !                                                                  */
    /*-----------------------------------------------------------------------------------------------------*/
    
	static int max_basic_agent_ID = 0; 
	ID = max_basic_agent_ID; // 
	max_basic_agent_ID++; 
	
	is_active=true;
	
	volume = 1.0; 
	
	position.assign( 3 , 0.0 );                        // initialize position and velocity
	velocity.assign( 3 , 0.0 );
	previous_velocity.assign( 3 , 0.0 ); 
                                                        // link into the microenvironment, if one is defined 
	secretion_rates= new std::vector<double>(0);
	uptake_rates= new std::vector<double>(0);
	saturation_densities= new std::vector<double>(0);
	extern Microenvironment* default_microenvironment;
	
	register_microenvironment( default_microenvironment ); 

	return;	
}

void Basic_Agent::update_position(double dt){ 
//make sure to update current_voxel_index if you are implementing this function
};

 /*---------------------------------------------------------------------------------------*/
 /* This is a dummy function which calls another version of assign_position(...) function */
 /* to set the x,y, and z coordinates of a Basic_Agent object 														*/
 /*---------------------------------------------------------------------------------------*/

bool Basic_Agent::assign_position(std::vector<double> new_position, int mpi_Rank, int *mpi_Dims) //-->Gaurav Saxena
{
	return assign_position(new_position[0], new_position[1], new_position[2], mpi_Rank, mpi_Dims);
}

/*-----------------------------------------------------------------------------------*/
/* Sets x,y, and z position of the Basic_Agent on a specific process and also checks */
/* if the basic_agent has gone out of the domain. Further, it updates the voxel 		 */
/* index in which the basic_agent now resides 																			 */
/*-----------------------------------------------------------------------------------*/

bool Basic_Agent::assign_position(double x, double y, double z, int mpi_Rank, int *mpi_Dims)
{
	// std::cout << __FILE__ << " " << __LINE__ << std::endl; 
	if( !get_microenvironment()->mesh.is_position_valid(x,y,z))    //No need to parallelize is_position_valid(x,y,z) - it will work for parallel implementation
	{	
		// std::cout<<"Error: the new position for agent "<< ID << " is invalid: "<<x<<","<<y<<","<<"z"<<std::endl;
		return false;
	}
	position[0]=x;
	position[1]=y;
	position[2]=z;
	update_voxel_index(mpi_Rank, mpi_Dims);
	
	// make sure the agent is not already registered
	get_microenvironment()->agent_container->register_agent(this);
	return true;
}

/*---------------------------------------------------------------*/
/* The function checks if the basic_agent is still within domain */
/* If yes then it updates the current_voxel_index of basic_agent */
/* The current_voxel_index is the LOCAL index and NOT the GLOBAL */
/* index 																												 */
/*---------------------------------------------------------------*/

void Basic_Agent::update_voxel_index(int mpi_Rank, int *mpi_Dims)
{
	
	if( !get_microenvironment()->mesh.is_position_valid(position[0],position[1],position[2]))
	{	
		std::cout << "Invalid position detected "  << position[0] << " " << position[1] << " " << position[2] << std::endl;
		std::cout << "ID " << ID << std::endl;
		current_voxel_index=-1;
		is_active=false;
		return;
	}
	int aux = get_microenvironment()->mesh.nearest_voxel_local_index( position, mpi_Rank, mpi_Dims );
	if (aux < 0) aux = 0;
	if (aux >= get_microenvironment()->mesh.voxels.size()) aux = get_microenvironment()->mesh.voxels.size() -1;
	//current_voxel_index= microenvironment->mesh.nearest_voxel_local_index( position, mpi_Rank, mpi_Dims );
	current_voxel_index = aux;
}

void Basic_Agent::set_internal_uptake_constants( double dt )
{
	// overall form: dp/dt = S*(T-p) - U*p 
	//   p(n+1) - p(n) = dt*S(n)*T(n) - dt*( S(n) + U(n))*p(n+1)
	//   p(n+1)*temp2 =  p(n) + temp1
	//   p(n+1) = (  p(n) + temp1 )/temp2
	//int nearest_voxel= current_voxel_index;
	double internal_constant_to_discretize_the_delta_approximation = dt * volume / ( (microenvironment->voxels(current_voxel_index)).volume ) ; // needs a fix 

	// temp1 = dt*(V_cell/V_voxel)*S*T 
	cell_source_sink_solver_temp1.assign( (*secretion_rates).size() , 0.0 ); 
	cell_source_sink_solver_temp1 += *secretion_rates; 
	cell_source_sink_solver_temp1 *= *saturation_densities; 
	cell_source_sink_solver_temp1 *= internal_constant_to_discretize_the_delta_approximation; 

	// temp2 = 1 + dt*(V_cell/V_voxel)*( S + U )
	cell_source_sink_solver_temp2.assign( (*secretion_rates).size() , 1.0 ); 
	axpy( &(cell_source_sink_solver_temp2) , internal_constant_to_discretize_the_delta_approximation , *secretion_rates );
	axpy( &(cell_source_sink_solver_temp2) , internal_constant_to_discretize_the_delta_approximation , *uptake_rates );	
	
	volume_is_changed = false; 
}

void Basic_Agent::register_microenvironment( Microenvironment* microenvironment_in )
{
	
	microenvironment = microenvironment_in; 	
	secretion_rates->resize( microenvironment->number_of_densities() , 0.0 );
	saturation_densities->resize( microenvironment->number_of_densities() , 0.0 );
	uptake_rates->resize( microenvironment->number_of_densities() , 0.0 );	

	// some solver temporary variables 
	cell_source_sink_solver_temp1.resize( microenvironment->number_of_densities() , 0.0 );
	cell_source_sink_solver_temp2.resize( microenvironment->number_of_densities() , 1.0 );
	
	return; 
}

Microenvironment* Basic_Agent::get_microenvironment( void )
{ return microenvironment; }

Basic_Agent::~Basic_Agent()
{
 return; 
}

Basic_Agent* create_basic_agent( void )
{
	Basic_Agent* pNew; 
	pNew = new Basic_Agent;	 
	all_basic_agents.push_back( pNew ); 
	pNew->index=all_basic_agents.size()-1;         //index is local to this process
	return pNew; 
}

void delete_basic_agent( int index )
{
	// deregister agent in microenvironment
	all_basic_agents[index]->get_microenvironment()->agent_container->remove_agent(all_basic_agents[index]);
	// de-allocate (delete) the Basic_Agent; 
	
	delete all_basic_agents[index]; 

	// next goal: remove this memory address. 

	// performance goal: don't delete in the middle -- very expensive reallocation
	// alternative: copy last element to index position, then shrink vector by 1 at the end O(constant)

	// move last item to index location  
	all_basic_agents[ all_basic_agents.size()-1 ]->index=index;
	all_basic_agents[index] = all_basic_agents[ all_basic_agents.size()-1 ];

	// shrink the vector
	all_basic_agents.pop_back();
	
	return; 
}

void delete_basic_agent( Basic_Agent* pDelete )
{
	// First, figure out the index of this agent. This is not efficient. 

	// int delete_index = 0; 
	// while( all_basic_agents[ delete_index ] != pDelete )
	// { delete_index++; }

	delete_basic_agent(pDelete->index);
	return; 
}

int Basic_Agent::get_current_voxel_index( void )
{
	return current_voxel_index;
}

double *Basic_Agent::nearest_density_vector( void ) 
{  
	return microenvironment->nearest_density_vector( current_voxel_index ); 
}


// directly access the gradient of substrate n nearest to the cell 
std::vector<double>& Basic_Agent::nearest_gradient( int substrate_index )
{
	return microenvironment->gradient_vector(current_voxel_index)[substrate_index]; 
}

	// directly access a vector of gradients, one gradient per substrate 
std::vector<gradient>& Basic_Agent::nearest_gradient_vector( void )
{
	return microenvironment->gradient_vector(current_voxel_index); 
}


void Basic_Agent::set_total_volume(double volume)
{
	this->volume = volume;
	volume_is_changed = true;
}

double Basic_Agent::get_total_volume()
{
	return volume;
}
void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt )
{
	
	if(!is_active)
	{ return; }
	
	if( volume_is_changed )
	{
		set_internal_uptake_constants(dt);
		volume_is_changed = false;
	}

	int d_size = (*pS).number_of_densities();
	for(int d = 0; d < d_size; ++d) {
		(*pS)(current_voxel_index)[d] += cell_source_sink_solver_temp1[d];
	}
	for(int d = 0; d < d_size; ++d) {
		(*pS)(current_voxel_index)[d] /= cell_source_sink_solver_temp2[d];
	}

	return; 
}

};