Skip to content

Commit

Permalink
notebook is running
Browse files Browse the repository at this point in the history
  • Loading branch information
eagmon committed Sep 12, 2024
1 parent f0595f6 commit 5f91532
Show file tree
Hide file tree
Showing 7 changed files with 714 additions and 144 deletions.
745 changes: 637 additions & 108 deletions demo/particle_comets.ipynb

Large diffs are not rendered by default.

4 changes: 1 addition & 3 deletions spatio_flux/processes/comets.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,11 @@
"""
COMETS composite made of dFBAs and diffusion-advection processes.
"""

from process_bigraph import Composite
from bigraph_viz import plot_bigraph
from spatio_flux import core
from spatio_flux.viz.plot import plot_time_series, plot_species_distributions_to_gif


# TODO -- need to do this to register???
from spatio_flux.processes.dfba import get_spatial_dfba_state
from spatio_flux.processes.diffusion_advection import get_diffusion_advection_spec
Expand Down Expand Up @@ -78,7 +76,7 @@ def run_comets(
# plot timeseries
plot_time_series(
comets_results,
coordinates=[(0, 0), (5, 5)],
coordinates=[(0, 0), (n_bins[0]-1, n_bins[1]-1)],
out_dir='out',
filename='comets_timeseries.png',
)
Expand Down
42 changes: 40 additions & 2 deletions spatio_flux/processes/dfba.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
Process for a pluggable dFBA simulation.
"""

import numpy as np
import warnings
import cobra
Expand Down Expand Up @@ -234,6 +233,45 @@ def get_spatial_dfba_state(
}


def run_dfba_single(
total_time=60,
mol_ids=None,
):
single_dfba_config = {
'dfba': get_single_dfba_spec(path=['fields']),
'fields': {
'glucose': 10,
'acetate': 0,
'biomass': 0.1
}
}

# make the simulation
sim = Composite({
'state': single_dfba_config,
'emitter': {'mode': 'all'}
}, core=core)

# save the document
sim.save(filename='single_dfba.json', outdir='out')

# simulate
print('Simulating...')
sim.update({}, total_time)

# gather results
dfba_results = sim.gather_results()

print('Plotting results...')
# plot timeseries
plot_time_series(
dfba_results,
# coordinates=[(0, 0), (1, 1), (2, 2)],
out_dir='out',
filename='dfba_single_timeseries.png',
)


def run_dfba_spatial(
total_time=60,
n_bins=(3, 3), # TODO -- why can't do (5, 10)??
Expand Down Expand Up @@ -271,7 +309,6 @@ def run_dfba_spatial(

# gather results
dfba_results = sim.gather_results()
# print(dfba_results)

print('Plotting results...')
# plot timeseries
Expand All @@ -293,4 +330,5 @@ def run_dfba_spatial(


if __name__ == '__main__':
run_dfba_single()
run_dfba_spatial(n_bins=(8,8))
8 changes: 0 additions & 8 deletions spatio_flux/processes/diffusion_advection.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
This process is a simple 2D diffusion-advection process. It takes a 2D field as input and returns a 2D field as output.
"""

import numpy as np
from scipy.ndimage import convolve
from process_bigraph import Process, Composite
Expand Down Expand Up @@ -122,9 +121,6 @@ def diffusion_delta(self, state, interval, diffusion_coeff, advection_coeff):
# Update the current state
updated_state += (laplacian + advective_flux_x + advective_flux_y) * dt

# # Ensure non-negativity
# current_states[species] = np.maximum(updated_state, 0)

# Update time
t += dt

Expand Down Expand Up @@ -199,9 +195,6 @@ def get_diffusion_advection_state(
'acetate': 0,
'biomass': 0.1
}
# initial_fields = {
# mol_id: np.random.uniform(low=0, high=initial_max[mol_id], size=n_bins[1], n_bins[0]))
# for mol_id in mol_ids}
initial_fields = {
mol_id: np.random.uniform(low=0, high=initial_max[mol_id], size=n_bins)
for mol_id in mol_ids}
Expand Down Expand Up @@ -267,7 +260,6 @@ def run_diffusion_process(

# gather results
diffadv_results = sim.gather_results()
# print(diffadv_results)

print('Plotting results...')
# plot 2d video
Expand Down
23 changes: 17 additions & 6 deletions spatio_flux/processes/particle_comets.py
Original file line number Diff line number Diff line change
@@ -1,38 +1,49 @@
"""
Particle-COMETS composite made of dFBAs, diffusion-advection, and particle processes.
"""

from process_bigraph import Composite
from bigraph_viz import plot_bigraph
from spatio_flux import core
from spatio_flux.viz.plot import plot_time_series, plot_species_distributions_with_particles_to_gif


# TODO -- need to do this to register???
from spatio_flux.processes.dfba import DynamicFBA, get_spatial_dfba_state
from spatio_flux.processes.diffusion_advection import DiffusionAdvection, get_diffusion_advection_spec
from spatio_flux.processes.particles import Particles, get_particles_spec, get_particles_state


default_config = {
'total_time': 10.0,
'total_time': 100.0,
# environment size
'bounds': (10.0, 20.0),
'n_bins': (8, 16),
# set fields
'mol_ids': ['glucose', 'acetate', 'biomass', 'detritus'],
'field_diffusion_rate': 1e-1,
'field_advection_rate': (0, 0),
'initial_min_max': {'glucose': (10, 10), 'acetate': (0, 0), 'biomass': (0, 0.1), 'detritus': (0, 0)},
'initial_min_max': {
'glucose': (10, 10),
'acetate': (0, 0),
'biomass': (0, 0.1),
'detritus': (0, 0)
},
# set particles
'n_particles': 10,
'particle_diffusion_rate': 1e-1,
'particle_advection_rate': (0, -0.1),
'particle_add_probability': 0.3,
'particle_boundary_to_add': ['top'],
'field_interactions': {
'biomass': {'vmax': 0.1, 'Km': 1.0, 'interaction_type': 'uptake'},
'detritus': {'vmax': -0.1, 'Km': 1.0, 'interaction_type': 'secretion'},
'biomass': {
'vmax': 0.1,
'Km': 1.0,
'interaction_type': 'uptake'
},
'detritus': {
'vmax': -0.1,
'Km': 1.0,
'interaction_type': 'secretion'
},
},
}

Expand Down
16 changes: 7 additions & 9 deletions spatio_flux/processes/particles.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
A process for simulating the motion of particles in a 2D environment.
"""

import uuid
import numpy as np
from process_bigraph import Process, Composite
Expand Down Expand Up @@ -239,13 +238,13 @@ def check_boundary_hit(self, new_x_position, new_y_position):

def get_boundary_position(self, boundary):
if boundary == 'left':
return (self.env_size[0][0], np.random.uniform(*self.env_size[1]))
return self.env_size[0][0], np.random.uniform(*self.env_size[1])
elif boundary == 'right':
return (self.env_size[0][1], np.random.uniform(*self.env_size[1]))
return self.env_size[0][1], np.random.uniform(*self.env_size[1])
elif boundary == 'top':
return (np.random.uniform(*self.env_size[0]), self.env_size[1][1])
return np.random.uniform(*self.env_size[0]), self.env_size[1][1]
elif boundary == 'bottom':
return (np.random.uniform(*self.env_size[0]), self.env_size[1][0])
return np.random.uniform(*self.env_size[0]), self.env_size[1][0]


core.register_process('Particles', Particles)
Expand Down Expand Up @@ -300,7 +299,7 @@ def get_particles_state(
}
if initial_min_max is None:
initial_min_max = {
'biomass': (0, 0.1),
'biomass': (0.1, 0.2),
'detritus': (0, 0),
}

Expand Down Expand Up @@ -329,7 +328,7 @@ def get_particles_state(


def run_particles(
total_time=20, # Total frames
total_time=100, # Total frames
bounds=(10.0, 20.0), # Bounds of the environment
n_bins=(20, 40), # Number of bins in the x and y directions
n_particles=20,
Expand Down Expand Up @@ -372,9 +371,8 @@ def run_particles(
# gather results
particles_results = sim.gather_results()
emitter_results = particles_results[('emitter',)]

# resort results
particles_history = [p['particles'] for p in emitter_results]
# print(particles_history)

print('Plotting...')
# plot particles
Expand Down
20 changes: 12 additions & 8 deletions spatio_flux/viz/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@

def sort_results(results):
emitter_results = results[('emitter',)]
if emitter_results[0] is None:
return
sorted_results = {'fields': {
key: [] for key in emitter_results[0]['fields'].keys()
}, 'time': []}
Expand Down Expand Up @@ -45,7 +47,7 @@ def plot_time_series(
:param coordinates: List of coordinates (indices) to plot.
Example: [(0, 0), (1, 2)]
"""
coordinates = coordinates or [(0, 0)]
# coordinates = coordinates or [(0, 0)]
field_names = field_names or ['glucose', 'acetate', 'biomass']
sorted_results = sort_results(results)
time = sorted_results['time']
Expand All @@ -56,13 +58,15 @@ def plot_time_series(
for field_name in field_names:
if field_name in sorted_results['fields']:
field_data = sorted_results['fields'][field_name]

for coord in coordinates:
x, y = coord
time_series = [field_data[t][x, y] for t in range(len(time))]
ax.plot(time, time_series, label=f'{field_name} at {coord}')
# plot log scale on y axis
# ax.set_yscale('log')
if coordinates is None:
ax.plot(time, field_data, label=field_name)
else:
for coord in coordinates:
x, y = coord
time_series = [field_data[t][x, y] for t in range(len(time))]
ax.plot(time, time_series, label=f'{field_name} at {coord}')
# plot log scale on y axis
# ax.set_yscale('log')
else:
print(f"Field '{field_name}' not found in results['fields']")

Expand Down

0 comments on commit 5f91532

Please sign in to comment.