From 6157bf3ecdbe2b0927716bde1287c9400811b31b Mon Sep 17 00:00:00 2001 From: ilija dukovski Date: Mon, 24 Oct 2016 10:32:28 -0400 Subject: [PATCH] Version 2.3.0. Mike Hasson's work on material dependent transport. Also .mat file support is here. --- .../src/edu/bu/segrelab/comets/Comets.java | 2 +- .../bu/segrelab/comets/CometsSimRunner.java | 2 +- .../edu/bu/segrelab/comets/fba/FBACell.java | 3 +- .../segrelab/comets/fba/FBACometsLoader.java | 264 ++++++++- .../edu/bu/segrelab/comets/fba/FBAWorld.java | 251 +++++++- .../bu/segrelab/comets/fba/FBAWorld3D.java | 128 +++- .../edu/bu/segrelab/comets/util/Utility.java | 551 +++++++++++++++++- 7 files changed, 1170 insertions(+), 31 deletions(-) diff --git a/comets_simplified/src/edu/bu/segrelab/comets/Comets.java b/comets_simplified/src/edu/bu/segrelab/comets/Comets.java index ae430e8..dba0cf1 100644 --- a/comets_simplified/src/edu/bu/segrelab/comets/Comets.java +++ b/comets_simplified/src/edu/bu/segrelab/comets/Comets.java @@ -94,7 +94,7 @@ public class Comets implements CometsConstants, * by each cell just runs through a diffusion routine. */ public static boolean DIFFUSION_TEST_MODE = false; - private String versionString = "2.2.3, 2 August 2016"; + private String versionString = "2.3.0, 24 October 2016"; // The setup pane private CometsSimRunner runner; diff --git a/comets_simplified/src/edu/bu/segrelab/comets/CometsSimRunner.java b/comets_simplified/src/edu/bu/segrelab/comets/CometsSimRunner.java index 15ca787..413058c 100644 --- a/comets_simplified/src/edu/bu/segrelab/comets/CometsSimRunner.java +++ b/comets_simplified/src/edu/bu/segrelab/comets/CometsSimRunner.java @@ -114,7 +114,7 @@ else if (c.getParameters().getNumLayers()>1) c.fireSimulationStateChangeEvent(new SimulationStateChangeEvent(SimulationStateChangeEvent.State.PAUSE)); } if (c.getParameters().showCycleTime()) - // System.out.println("Cycle complete in " + (System.currentTimeMillis() - start)/1000.0 + "s"); + System.out.println("Cycle complete in " + (System.currentTimeMillis() - start)/1000.0 + "s"); totalTime += System.currentTimeMillis() - start; c.updateOnCycle(); } diff --git a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACell.java b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACell.java index 733a566..b97e00f 100644 --- a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACell.java +++ b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACell.java @@ -590,7 +590,6 @@ public synchronized int run(FBAModel[] models) media = world.getModelMediaAt(x, y, i); else if (cParams.getNumLayers() > 1) media = world3D.getModelMediaAt(x, y, z, i); - double[] lb = ((FBAModel)models[i]).getBaseExchLowerBounds(); double[] ub = ((FBAModel)models[i]).getBaseExchUpperBounds(); // String[] exchNames = ((FBAModel)models[i]).getExchangeReactionNames(); @@ -757,7 +756,7 @@ else if (cParams.getNumLayers() > 1) if(deltaBiomass[i]<0.0)deltaBiomass[i]=0.0; // deltaBiomass[i] = (double)(((FBAModel)models[i]).getObjectiveFluxSolution()) * cParams.getTimeStep(); // deltaBiomass[i] = (double)(((FBAModel)models[i]).getObjectiveFluxSolution()); - //System.out.println("solution: " + ((FBAModel)models[i]).getObjectiveSolution()); +// System.out.println("solution: " + ((FBAModel)models[i]).getObjectiveSolution()); if (cParams.showGraphics()) cellColor = calculateColor(); diff --git a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACometsLoader.java b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACometsLoader.java index 7813798..a71da28 100644 --- a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACometsLoader.java +++ b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBACometsLoader.java @@ -70,7 +70,15 @@ private enum LoaderState private boolean useGui; private Comets c; private int lineCount; - + private boolean substrate = false; + private boolean modelDiffusion = false; + private boolean specific = false; + private boolean friction = false; + private double[][] substrateDiffConsts; + private double[][] modelDiffConsts; + private int[][] substrateLayout; + private double[][] specificMedia; + private double[][] substrateFrictionConsts; private static final String MODEL_FILE = "model_file", MODEL_WORLD = "model_world", GRID_SIZE = "grid_size", @@ -87,7 +95,12 @@ private enum LoaderState BARRIER = "barrier", MEDIA = "media", PARAMETERS = "parameters", - DIFFUSION_CONST = "diffusion_constants"; + DIFFUSION_CONST = "diffusion_constants", + SUBSTRATE_DIFFUSIVITY = "substrate_diffusivity", + SUBSTRATE_FRICTION = "substrate_friction", + MODEL_DIFFUSIVITY = "model_diffusivity", + SUBSTRATE_LAYOUT = "substrate_layout", + SPECIFIC_MEDIA = "specific_media"; /** * Returns the recently loaded World2D. */ @@ -202,7 +215,7 @@ public int loadLayoutFile(String filename, Comets c, boolean useGui) throws IOEx /* * makes and stores internally: an FBAWorld, an array of FBAModels, an * ArrayList of FBACells - */ + */ try { BufferedReader reader = new BufferedReader(new FileReader(filename)); @@ -270,7 +283,6 @@ public int loadLayoutFile(String filename, Comets c, boolean useGui) throws IOEx if (state == LoaderState.CANCELED) throw new LayoutFileException("Layout file loading canceled", lineCount); } - else if (parsed[0].equalsIgnoreCase(MODEL_WORLD)) { Map media = new HashMap(); @@ -363,10 +375,10 @@ else if (worldParsed[0].equalsIgnoreCase(GRID_SIZE)) else if(worldParsed.length == 3) c.getParameters().setNumLayers(Integer.valueOf(1)); - if (Integer.valueOf(worldParsed[3])==1) - { - throw new IOException("In a 3D simulation the number of layers along the 3rd coordinate must be larger than 1."); - } + //if (Integer.valueOf(worldParsed[3])==1) + //{ + // throw new IOException("In a 3D simulation the number of layers along the 3rd coordinate must be larger than 1."); + //} } /****************** BARRIER *************************/ @@ -396,6 +408,47 @@ else if (worldParsed[0].equalsIgnoreCase(DIFFUSION_CONST)) List lines = collectLayoutFileBlock(reader); state = parseMediaDiffusionConstantsBlock(lines, diffConsts); } + /****************** DIFFUSION CONSTANTS BY SUBSTRATE ********************/ + + else if (worldParsed[0].equalsIgnoreCase(SUBSTRATE_DIFFUSIVITY)) + { + substrate = true; + List lines = collectLayoutFileBlock(reader); + state = parseSubstrateDiffusionConstantsBlock(lines,numMedia); + } + /****************** FRICTION CONSTANTS BY SUBSTRATE ********************/ + + else if (worldParsed[0].equalsIgnoreCase(SUBSTRATE_FRICTION)) + { + friction = true; + List lines = collectLayoutFileBlock(reader); + state = parseSubstrateFrictionConstantsBlock(lines); + } + /****************** DIFFUSION CONSTANTS BY MODEL ********************/ + + else if (worldParsed[0].equalsIgnoreCase(MODEL_DIFFUSIVITY)) + { + modelDiffusion = true; + List lines = collectLayoutFileBlock(reader); + state = parseModelDiffusionConstantsBlock(lines,numMedia); + } + + /****************** SUBSTRATE LAYOUT ********************/ + + else if (worldParsed[0].equalsIgnoreCase(SUBSTRATE_LAYOUT)) + { + List lines = collectLayoutFileBlock(reader); + state = parseSubstrateLayoutBlock(lines,c.getParameters().getNumCols(),c.getParameters().getNumRows()); + } + + /****************** SPECIFIC MEDIA ********************/ + + else if (worldParsed[0].equalsIgnoreCase(SPECIFIC_MEDIA)) + { + specific = true; + List lines = collectLayoutFileBlock(reader); + state = parseSpecificMediaBlock(lines); + } /****************** MEDIA ***********************/ @@ -634,7 +687,6 @@ else if (worldParsed[0].equalsIgnoreCase("prevent_media_in")) } } } - /* Set diffusion constants */ double[] diffusionConsts = new double[numMedia]; double defaultDiffConst = pParams.getDefaultDiffusionConstant(); @@ -649,6 +701,22 @@ else if (worldParsed[0].equalsIgnoreCase("prevent_media_in")) } } world.setDiffusionConstants(diffusionConsts); + if(substrate) + { + world.setSubstrateDiffusion(substrateDiffConsts); + world.setSubstrateLayout(substrateLayout); + } + if(modelDiffusion){ + world.setModelDiffusivity(modelDiffConsts); + } + if(specific) + { + world.setSpecificMedia(specificMedia); + } + if(friction) + { + world.setSubstrateFriction(substrateFrictionConsts); + } System.out.println("Done!"); } else if(c.getParameters().getNumLayers()>1) @@ -1109,6 +1177,184 @@ private LoaderState parseMediaDiffusionConstantsBlock(List lines, Map lines, int numMedia) throws LayoutFileException, + NumberFormatException + { + /* the 'diffusion_constants' block is taken from the way of doing it in the Model files, so it looks like this: + * + * diffusion_constants + * + * + * // + */ + Integer i = 0; + substrateDiffConsts = new double[lines.size()][numMedia]; + for (String line : lines) + { + lineCount++; + if (line.length() == 0) + continue; + + String[] diffConstParsed = line.split("\\s+"); + if (diffConstParsed.length != numMedia) + throw new LayoutFileException("Each line of the 'diffusion_constants' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + + else + { + for(int j = 0;j lines) throws LayoutFileException, + NumberFormatException + { + /* the 'diffusion_constants' block is taken from the way of doing it in the Model files, so it looks like this: + * + * diffusion_constants + * + * + * // + */ + Integer i = 0; + substrateFrictionConsts = new double[lines.size()][models.length]; + for (String line : lines) + { + lineCount++; + if (line.length() == 0) + continue; + + String[] diffConstParsed = line.split("\\s+"); + if (diffConstParsed.length != models.length) + throw new LayoutFileException("Each line of the 'diffusion_constants' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + + else + { + for(int j = 0;j lines, int numMedia) throws LayoutFileException, + NumberFormatException + { + /* the 'diffusion_constants' block is taken from the way of doing it in the Model files, so it looks like this: + * + * diffusion_constants + * + * + * // + */ + Integer i = 0; + modelDiffConsts = new double[lines.size()][numMedia]; + for (String line : lines) + { + lineCount++; + if (line.length() == 0) + continue; + + String[] diffConstParsed = line.split("\\s+"); + if (diffConstParsed.length != numMedia) + throw new LayoutFileException("Each line of the 'diffusion_constants' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + + else + { + for(int j = 0;j lines,int cols,int rows) throws LayoutFileException, + NumberFormatException + { + /* the 'diffusion_constants' block is taken from the way of doing it in the Model files, so it looks like this: + * + * diffusion_constants + * + * + * // + */ + if (lines.size() != rows) + throw new LayoutFileException("Each line of the 'diffusion_constants' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + else + { + Integer i = 0; + substrateLayout = new int[cols][rows]; + for (String line : lines) + { + lineCount++; + if (line.length() == 0) + continue; + + String[] layoutParsed = line.split("\\s+"); + if (layoutParsed.length != cols) + throw new LayoutFileException("Each line of the 'diffusion_constants' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + + else + { + for(int j = 0;j lines) throws LayoutFileException, + NumberFormatException + { + /* the 'diffusion_constants' block is taken from the way of doing it in the Model files, so it looks like this: + * + * diffusion_constants + * + * + * // + */ + Integer i = 0; + specificMedia = new double[4][lines.size()]; + for (String line : lines) + { + lineCount++; + if (line.length() == 0) + continue; + + String[] sMediaParsed = line.split("\\s+"); + if (sMediaParsed.length != 4) + throw new LayoutFileException("Each line of the 'specific_media' block should have two tokens.\n The first should be the index of the medium component, followed by its (non-default) diffusion constant.", lineCount); + + else + { + for(int j = 0;j<4;j++) + { + specificMedia[j][i] = Double.parseDouble(sMediaParsed[j]); + } + i++; + } + } + return LoaderState.OK; + } + private LoaderState parseInitialPopBlock(String[] header, List lines) throws LayoutFileException, NumberFormatException { diff --git a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBAWorld.java b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBAWorld.java index 61b94a4..17edcfb 100644 --- a/comets_simplified/src/edu/bu/segrelab/comets/fba/FBAWorld.java +++ b/comets_simplified/src/edu/bu/segrelab/comets/fba/FBAWorld.java @@ -67,6 +67,8 @@ public class FBAWorld extends World2D private boolean[][][] diffuseBiomassOut; // as above, but diffusing out private boolean[][][] diffuseMediaIn; // as above, but for media private boolean[][][] diffuseMediaOut; + private double[][][] diffusionRHS1; + private double[][][] diffusionRHS2; private int threadLock; // while threaded FBA is running, this = number of cells remaining private FBAParameters pParams; @@ -98,6 +100,14 @@ public class FBAWorld extends World2D private double defaultDiffConst; + private int[][] substrateLayout; + private double[][] modelDiffusivity; + private boolean diffuseContext = false; + private boolean frictionContext = false; + private boolean diffuseModel = false; + private FBASubstrate[] substrates; + + private int numSubstrates; /** * Initialize a new, empty world, tied the current Comets with a given * number of medium components @@ -134,6 +144,8 @@ public FBAWorld(Comets c, String[] mediaNames, double[] startingMedia, dirichlet = new boolean[numCols][numRows]; diffuseMediaIn = new boolean[numCols][numRows][numMedia]; diffuseMediaOut = new boolean[numCols][numRows][numMedia]; + diffusionRHS1 = new double[numMedia][numCols][numRows]; + diffusionRHS2 = new double[numMedia][numCols][numRows]; diffuseBiomassIn = new boolean[numCols][numRows][numModels]; diffuseBiomassOut = new boolean[numCols][numRows][numModels]; nutrientDiffConsts = new double[numMedia]; @@ -149,6 +161,8 @@ public FBAWorld(Comets c, String[] mediaNames, double[] startingMedia, for (int k = 0; k < numMedia; k++) { media[i][j][k] = startingMedia[k]; + diffusionRHS1[k][i][j] = 0; + diffusionRHS2[k][i][j] = 0; diffuseMediaIn[i][j][k] = true; diffuseMediaOut[i][j][k] = true; } @@ -1449,6 +1463,129 @@ private void diffuseMediaFick() } } } + /** + * Diffuses media according to RHS and predictor-corrector 2D system. Each media layer + * is diffused separately, by calling Utility.getDiffusionRHS() on each one in turn. + */ + private void diffuseMediaContext() + { +// long time = System.currentTimeMillis(); + /** + * Foreach of the media types, copy each media layer into a + * new double[][], pipe over to Utility.diffusionFick(), then + * copy the solutions back into that layer. + * + * The added io time for copying, etc, might make it take longer, + * but it will more than make up for it in readability/ + * maintainability. + */ + + double dT = cParams.getTimeStep() * 3600/pParams.getNumDiffusionsPerStep(); // time units = hours ( as in fba ), convert to seconds + double dX = cParams.getSpaceWidth(); + double[][][] diffusionRHS = new double[numMedia][numCols][numRows]; + + if (DEBUG) + { + System.out.println("dT: " + dT + "\tdX: " + dX); + System.out.println("media diff consts:"); + for (int i=0; i=0){ + media[(int)specificMedia[0][j]][(int)specificMedia[1][j]][(int)specificMedia[2][j]] = specificMedia[3][j]; + }else{ + setSubstrateMedia((int)specificMedia[1][j],(int)specificMedia[2][j],specificMedia[3][j]); + } + } + } + public void setSubstrateMedia(int substrate, int mediaNum, double amount) + { + for (int i=0; ipackDensity*dX*dX) + { + pressure[i][j]=elasticModulusConst*Math.pow(1.0-packDensity*dX*dX/biomass[i][j],1.5); + } + else + { + pressure[i][j]=0.0; + } + + //System.out.println("pressure"+" "+i+","+j+" "+pressure[i][j]); + } + } + return pressure; + } + + /** + * Calculates the pressure according to the Farrel et. al. model: PRL 111, 168101(2013). + * @param biomass + * @return + */ + public static double[][] pressure2Dc(double[][] biomass,double elasticModulusConst, double[][] frictionConst, double packDensity) { double[][] pressure=new double[biomass.length][biomass[0].length]; for(int i=0;ipackDensity) { - pressure[i][j]=elasticModulusConst*Math.pow(1.0-packDensity/biomass[i][j],1.5); + if (frictionConst[i][j] > 0){ + pressure[i][j]=(elasticModulusConst/frictionConst[i][j])*Math.pow(1.0-packDensity/biomass[i][j],1.5); + } else{ + pressure[i][j] = 0; + } } else { pressure[i][j]=0.0; } + //System.out.println(i + " " + j + " " + pressure[i][j]); + //System.out.println("pressure"+" "+i+","+j+" "+pressure[i][j]); + } + } + return pressure; + } + + /** + * Calculates the pressure according to the Farrel et. al. model: PRL 111, 168101(2013). + * @param biomass + * @return + */ + public static double[][] pressure2Da(double[][] biomass,double elasticModulusConst, double[][] frictionConst, double packDensity, double dX) + { + double[][] pressure=new double[biomass.length][biomass[0].length]; + for(int i=0;ipackDensity*dX*dX) + { + pressure[i][j]=Math.pow(1.0-packDensity*dX*dX/biomass[i][j],1.5); + } + else + { + pressure[i][j]=0.0; + } + //System.out.println(i + " " + j + " " + pressure[i][j]); //System.out.println("pressure"+" "+i+","+j+" "+pressure[i][j]); } } return pressure; } + /** + * Calculates the pressure according to the Farrel et. al. model: PRL 111, 168101(2013). + * @param biomass + * @return + */ + public static double[][] frictElast2D(double[][] biomass,double elasticModulusConst, double[][] frictionConst, double packDensity) + { + double[][] feConst=new double[biomass.length][biomass[0].length]; + for(int i=0;i 0){ + feConst[i][j]=(elasticModulusConst/frictionConst[i][j]); + } else{ + feConst[i][j] = 0; + } + + } + //System.out.println(i + " " + j + " " + pressure[i][j]); + //System.out.println("pressure"+" "+i+","+j+" "+pressure[i][j]); + } + return feConst; + } + /** * 3D case: Calculates the pressure according to the Farrel et. al. model: PRL 111, 168101(2013). * @param biomass