diff --git a/src/arcade/core/util/Graph.java b/src/arcade/core/util/Graph.java index 77e37bfc8..25c292f96 100644 --- a/src/arcade/core/util/Graph.java +++ b/src/arcade/core/util/Graph.java @@ -31,6 +31,9 @@ public enum Strategy { /** Collection of all {@code Edge} objects in a graph. */ private final Bag allEdges; + /** Map of {@code Node} objects, for lookup. */ + private final Map nodes; + /** Map of {@code Node} OUT to bag of {@code Edge} objects. */ private final Map nodeToOutBag; @@ -42,6 +45,7 @@ public Graph() { allEdges = new Bag(); nodeToOutBag = new HashMap<>(); nodeToInBag = new HashMap<>(); + nodes = new HashMap<>(); } /** @@ -53,6 +57,36 @@ public void update(Graph graph) { allEdges.addAll(graph.allEdges); nodeToOutBag.putAll(graph.nodeToOutBag); nodeToInBag.putAll(graph.nodeToInBag); + updateNodes(); + } + + /** Updates nodes from graph. */ + private void updateNodes() { + nodes.clear(); + for (Object obj : allEdges) { + Edge edge = (Edge) obj; + Node from = edge.getFrom(); + Node to = edge.getTo(); + + if (!nodes.containsKey(from.toString())) { + nodes.put(from.toString(), from); + } + if (!nodes.containsKey(to.toString())) { + nodes.put(to.toString(), to); + } + } + } + + public void combine(Graph graph1, Graph graph2) { + clear(); + for (Object obj : graph2.getAllEdges()) { + Edge edge = (Edge) obj; + addEdge(edge); + } + for (Object obj : graph1.getAllEdges()) { + Edge edge = (Edge) obj; + addEdge(edge); + } } /** Clear edges and nodes from graph. */ @@ -60,6 +94,7 @@ public void clear() { allEdges.clear(); nodeToOutBag.clear(); nodeToInBag.clear(); + nodes.clear(); } /** @@ -91,6 +126,28 @@ public boolean contains(Edge edge) { return checkEdge(edge); } + /** + * Retrieve the node object for the given coordinates. Returns null if no node exists. + * + * @param x the x coordinate + * @param y the y coordinate + * @param z the z coordinate + * @return the node object + */ + public Node lookup(int x, int y, int z) { + return nodes.get("(" + x + "," + y + "," + z + ")"); + } + + /** + * Retrieve the node object for the given coordinates. Returns null if no node exists. + * + * @param node to lookup an original node + * @return the node object + */ + public Node lookup(Node node) { + return nodes.get(node.toString()); + } + /** * Gets all edges in the graph. * @@ -206,6 +263,7 @@ public void getSubgraph(Graph g, GraphFilter f) { for (Object obj : allEdges) { Edge edge = (Edge) obj; if (f.filter(edge)) { + g.addNodes(edge); g.allEdges.add(edge); g.setOutMap(edge.getFrom(), edge); g.setInMap(edge.getTo(), edge); @@ -255,6 +313,8 @@ public void mergeNodes() { e.setTo(join); } } + + nodes.put(join.toString(), join); } } @@ -274,14 +334,44 @@ public void addEdge(Node from, Node to) { * @param edge the edge to add */ public void addEdge(Edge edge) { + addNodes(edge); allEdges.add(edge); setOutMap(edge.getFrom(), edge); setInMap(edge.getTo(), edge); setLinks(edge); } + public void addEdge(Edge edge, boolean duplicate) { + allEdges.add(edge); + setOutMap(edge.getFrom(), edge, duplicate); + setInMap(edge.getTo(), edge, duplicate); + setLinks(edge); + addNodes(edge); + } + /** - * Adds the edge to the bag for the mapping of OUT node to edge. + * Helper function to adds the nodes from an edge to a graph. + * + * @param edge the edge to add + */ + private void addNodes(Edge edge) { + Node from = edge.getFrom(); + Node to = edge.getTo(); + if (!nodes.containsKey(from.toString())) { + nodes.put(from.toString(), from); + } else { + edge.setFrom(nodes.get(from.toString())); + } + if (!nodes.containsKey(to.toString())) { + nodes.put(to.toString(), to); + } else { + edge.setTo(nodes.get(to.toString())); + } + } + + /** + * Adds the edge to the bag for the mapping of OUT node to edge. Default behavior is to + * duplicate nodes. * * @param node the node hash * @param edge the edge @@ -296,7 +386,25 @@ private void setOutMap(Node node, Edge edge) { } /** - * Adds the edge to the bag for the mapping of IN node to edge. + * Adds the edge to the bag for the mapping of OUT node to edge. Used in cases new node object + * is not needed. + * + * @param node the node hash + * @param edge the edge + */ + private void setOutMap(Node node, Edge edge, boolean duplicate) { + Bag objs = nodeToOutBag.get(node); + if (objs == null) { + objs = new Bag(10); + Node bagNode = duplicate ? node.duplicate() : node; + nodeToOutBag.put(bagNode, objs); + } + objs.add(edge); + } + + /** + * Adds the edge to the bag for the mapping of IN node to edge. Default behavior is to duplicate + * nodes. * * @param node the node hash * @param edge the edge @@ -310,6 +418,23 @@ private void setInMap(Node node, Edge edge) { objs.add(edge); } + /** + * Adds the edge to the bag for the mapping of OUT node to edge. Used in cases new node object + * is not needed. + * + * @param node the node hash + * @param edge the edge + */ + private void setInMap(Node node, Edge edge, boolean duplicate) { + Bag objs = nodeToInBag.get(node); + if (objs == null) { + objs = new Bag(10); + Node bagNode = duplicate ? node.duplicate() : node; + nodeToInBag.put(bagNode, objs); + } + objs.add(edge); + } + /** * Adds links between edges in and out of the nodes for a given edge. * @@ -353,6 +478,19 @@ public void removeEdge(Edge edge) { unsetOutMap(edge.getFrom(), edge); unsetInMap(edge.getTo(), edge); unsetLinks(edge); + removeNodeIfDetached(edge.getFrom()); + removeNodeIfDetached(edge.getTo()); + } + + /** + * Remove a node if it is detached from the graph. + * + * @param node the node to check + */ + private void removeNodeIfDetached(Node node) { + if (!nodeToInBag.containsKey(node) && !nodeToOutBag.containsKey(node)) { + nodes.remove(node.toString()); + } } /** @@ -730,7 +868,8 @@ public static class Edge { private final ArrayList edgesOut; /** - * Creates an {@code Edge} between two {@link Node} objects. + * Creates an {@code Edge} between two {@link Node} objects. Default behavior is to + * duplicate nodes. * * @param from the node the edge is from * @param to the node the edge is to @@ -742,6 +881,21 @@ public Edge(Node from, Node to) { edgesOut = new ArrayList<>(); } + /** + * Creates an {@code Edge} object for graph sites. Used in cases where a new node object is + * not needed. + * + * @param from the node the edge is from + * @param to the node the edge is to + * @param duplicate {@code true} if nodes should be duplicated, {@code false} otherwise + */ + public Edge(Node from, Node to, boolean duplicate) { + this.from = duplicate ? from.duplicate() : from; + this.to = duplicate ? to.duplicate() : to; + edgesIn = new ArrayList<>(); + edgesOut = new ArrayList<>(); + } + /** * Gets the node the edge points to based on the calculation strategy (e.g. upstream or * downstream). @@ -822,7 +976,7 @@ public ArrayList getEdgesOut() { * * @return the reversed edge */ - Edge reverse() { + public Edge reverse() { Node tempTo = to; Node tempFrom = from; to = tempFrom; diff --git a/src/arcade/core/util/MatrixArray.java b/src/arcade/core/util/MatrixArray.java new file mode 100644 index 000000000..c66e512ce --- /dev/null +++ b/src/arcade/core/util/MatrixArray.java @@ -0,0 +1,129 @@ +package arcade.core.util; + +import java.util.ArrayList; +import java.util.Collections; +import arcade.core.util.Matrix.Value; + +/** + * Container class for array-based sparse matrix representation. + * + *

Class provides a subset of matrix operations needed for solving a system of linear equations + * using the successive over-relaxation method in {@link arcade.core.util.Solver}. + */ +public class MatrixArray { + public MatrixArray(double[][] a) { + nRows = a.length; + if (nRows == 0) { + throw new IllegalArgumentException("MatrixArray ctor: input is empty"); + } + nColumns = a[0].length; + + if (nColumns == 0) { + throw new IllegalArgumentException("Matrix array ctor: empty columns"); + } + + ArrayList alValues = new ArrayList<>(); + int nNonZero = 0; + for (int row = 0; row < nRows; row++) { + if (a[row].length != nColumns) { + throw new IllegalArgumentException( + "Matrix array ctor: not all columns are the same length"); + } + + for (int column = 0; column < nColumns; column++) { + if (a[row][column] != 0.) { + alValues.add(new Value(row, column, a[row][column])); + nNonZero++; + } + } + } + + // + // No need to sort, since we built them in order. + // + + values = new Value[nNonZero]; + int i = 0; + for (Value v : alValues) { + values[i] = v; + i++; + } + } // ctor + + public MatrixArray(ArrayList alValues, int i_nRows, int i_Columns) { + nRows = i_nRows; + nColumns = i_Columns; + values = new Value[alValues.size()]; + + Collections.sort( + alValues, + (v1, v2) -> (v1.i == v2.i ? Integer.compare(v1.j, v2.j) : (v1.i > v2.i ? 1 : -1))); + + int i = 0; + for (Value v : alValues) { + values[i] = v; + i++; + } + } + + /** + * Solves the equation {@code Lx = b} using forward substitution for an array-based sparse + * matrix. + * + * @param b the right-hand side vector + * @return the left-hand side vector + */ + public double[] forwardSubstitution(double[] b) { + + int n = b.length; + double[] subbed = new double[n]; + double[] diag = new double[n]; + + // Group lower diagonal by row. + ArrayList> rowsL = new ArrayList>(); + for (int r = 0; r < n; r++) { + rowsL.add(new ArrayList<>()); + } + for (Value v : values) { + rowsL.get(v.i).add(v); + } + + // Get values along diagonal. + for (Value v : values) { + if (v.i == v.j) { + diag[v.i] = v.v; + } + } + + // Iterate only through non-zero entries in the lower diagonal matrix. + for (int i = 0; i < n; i++) { + ArrayList rowL = rowsL.get(i); + double val = 0; + for (Value v : rowL) { + val += subbed[v.j] * v.v; + } + val = b[i] - val; + subbed[i] = val / diag[i]; + } + + return subbed; + } + + public double[] multiply(double[] b) { + if (b.length != nRows) { + throw new IllegalArgumentException("MatrixArray.multiply (by a vector): conformation"); + } + double[] multiplied = new double[nRows]; + + // Iterate through all entries and multiply. + for (Value a : values) { + multiplied[a.i] += a.v * b[a.j]; + } + + return multiplied; + } + + int nRows; + int nColumns; + Value[] values; +} diff --git a/src/arcade/core/util/Solver.java b/src/arcade/core/util/Solver.java index 019aefff8..214c299bf 100644 --- a/src/arcade/core/util/Solver.java +++ b/src/arcade/core/util/Solver.java @@ -38,13 +38,13 @@ public class Solver { private static final double OMEGA = 1.4; /** Maximum number of iterations. */ - private static final int MAX_ITERS = 10000; + private static final int MAX_ITERS = 20000; /** Error tolerance for SOR. */ - private static final double TOLERANCE = 1E-8; + private static final double TOLERANCE = 1E-5; /** Convergence delta for bisection method. */ - private static final double DELTA = 1E-5; + private static final double DELTA = 1E-6; /** Matrix size threshold for dense representation. */ private static final int MATRIX_THRESHOLD = 100; @@ -370,7 +370,6 @@ private static double[] denseSOR( double[] r = subtract(vec, multiply(mat, xCurr)); error = normalize(r); } - return xCurr; } @@ -396,6 +395,7 @@ private static double[] sparseSOR( double[] c = forwardSubstitution(sparseA, vec); ArrayList t = forwardSubstitution(sparseA); t = scale(t, -1); + MatrixArray tArray = new MatrixArray(t, mat.length, mat[0].length); // Set initial guess. double[] xCurr = x0; @@ -404,7 +404,7 @@ private static double[] sparseSOR( // Iterate until convergence. while (i < maxIters && error > tolerance) { // Calculate new guess for x. - xCurr = add(scale(add(multiply(t, xPrev), c), OMEGA), scale(xPrev, 1 - OMEGA)); + xCurr = add(scale(add(tArray.multiply(xPrev), c), OMEGA), scale(xPrev, 1 - OMEGA)); // Set previous to copy of current and increment iteration count. xPrev = xCurr; @@ -414,7 +414,6 @@ private static double[] sparseSOR( double[] r = subtract(vec, multiply(sparseA, xCurr)); error = normalize(r); } - return xCurr; } @@ -428,9 +427,11 @@ private static double[] sparseSOR( * @param a the lower bound on the interval * @param b the upper bound on the interval * @param maxIters the maximum number of iterations + * @param tolerance the error * @return the root of the function */ - public static double bisection(Function func, double a, double b, int maxIters) { + public static double bisection( + Function func, double a, double b, int maxIters, double tolerance) { double c; double fc; int i = 0; @@ -452,7 +453,7 @@ public static double bisection(Function func, double a, double b, int maxIters) fc = func.f(c); // Check for exit conditions. - if (fc == 0 || (b - a) / 2 < DELTA) { + if (fc == 0 || (b - a) / 2 < tolerance) { return c; } else { if (Math.signum(fc) == Math.signum(func.f(a))) { @@ -469,7 +470,7 @@ public static double bisection(Function func, double a, double b, int maxIters) } /** - * Finds root using bisection method with default maximum iterations. + * Finds root using bisection method with default maximum iterations and tolerance. * *

Root is found by repeatedly bisecting the interval and selecting the interval in which the * function changes sign. If no root is found, the simulation will throw an ArithmeticException. @@ -480,6 +481,38 @@ public static double bisection(Function func, double a, double b, int maxIters) * @return the root of the function */ public static double bisection(Function func, double a, double b) { - return bisection(func, a, b, MAX_ITERS); + return bisection(func, a, b, MAX_ITERS, DELTA); + } + + /** + * Finds root using bisection method with default maximum iterations. + * + *

Root is found by repeatedly bisecting the interval and selecting the interval in which the + * function changes sign. If no root is found, the simulation will throw an ArithmeticException. + * + * @param func the function + * @param a the lower bound on the interval + * @param b the upper bound on the interval + * @param delta the error tolerance + * @return the root of the function + */ + public static double bisection(Function func, double a, double b, double delta) { + return bisection(func, a, b, MAX_ITERS, delta); + } + + /** + * Finds root using bisection method with default tolerance. + * + *

Root is found by repeatedly bisecting the interval and selecting the interval in which the + * function changes sign. If no root is found, the simulation will throw an ArithmeticException. + * + * @param func the function + * @param a the lower bound on the interval + * @param b the upper bound on the interval + * @param maxIters the maximum number of iterations + * @return the root of the function + */ + public static double bisection(Function func, double a, double b, int maxIters) { + return bisection(func, a, b, maxIters, DELTA); } } diff --git a/src/arcade/core/util/exceptions/IncompatibleFeatureException.java b/src/arcade/core/util/exceptions/IncompatibleFeatureException.java new file mode 100644 index 000000000..da417158d --- /dev/null +++ b/src/arcade/core/util/exceptions/IncompatibleFeatureException.java @@ -0,0 +1,22 @@ +package arcade.core.util.exceptions; + +/** Exception thrown when incompatible features are given. */ +public class IncompatibleFeatureException extends RuntimeException { + /** + * Constructs an {@code IncompatibleFeatureException} with the specified detail message. + * + * @param feature the feature name + * @param given the given associated feature + * @param expected the expected associated feature + */ + public IncompatibleFeatureException(String feature, String given, String expected) { + super( + String.format( + feature + + "is incompatible with " + + given + + ", must be associated with " + + expected + + ".")); + } +} diff --git a/src/arcade/core/util/exceptions/MissingSpecificationException.java b/src/arcade/core/util/exceptions/MissingSpecificationException.java new file mode 100644 index 000000000..4a175a798 --- /dev/null +++ b/src/arcade/core/util/exceptions/MissingSpecificationException.java @@ -0,0 +1,13 @@ +package arcade.core.util.exceptions; + +/** Exception thrown when a specification is missing. */ +public class MissingSpecificationException extends RuntimeException { + /** + * Constructs an {@code MissingSpecificationException} with the specified detail message. + * + * @param missing the expected specification + */ + public MissingSpecificationException(String missing) { + super("The input file must contain the following specification: " + missing); + } +} diff --git a/src/arcade/patch/agent/cell/PatchCell.java b/src/arcade/patch/agent/cell/PatchCell.java index 4a84bb520..0cd02cd06 100644 --- a/src/arcade/patch/agent/cell/PatchCell.java +++ b/src/arcade/patch/agent/cell/PatchCell.java @@ -454,6 +454,9 @@ public PatchLocation selectBestLocation(Simulation sim, MersenneTwisterFast rand // Calculate score by introducing error to the location check // and adding affinity to move toward center. double normConc = sim.getLattice("GLUCOSE").getAverageValue(location) / maxGlucose; + if (Double.isNaN(normConc)) { + normConc = 0; + } double gluc = (accuracy * normConc + (1 - accuracy) * random.nextDouble()); double dist = ((currR - loc.getPlanarDistance()) + 1) / 2.0; double score = affinity * dist + (1 - affinity) * gluc; diff --git a/src/arcade/patch/agent/cell/PatchCellTissue.java b/src/arcade/patch/agent/cell/PatchCellTissue.java index 352260b28..398b05147 100644 --- a/src/arcade/patch/agent/cell/PatchCellTissue.java +++ b/src/arcade/patch/agent/cell/PatchCellTissue.java @@ -134,6 +134,11 @@ public void step(SimState simstate) { if (module != null) { module.step(simstate.random, sim); } + + // Step remaining processes. + if (processes.get(Domain.SENSING) != null) { + processes.get(Domain.SENSING).step(simstate.random, sim); + } } /** diff --git a/src/arcade/patch/agent/process/PatchProcessSensingSimple.java b/src/arcade/patch/agent/process/PatchProcessSensingSimple.java index c81ba8cca..c9d2c3895 100644 --- a/src/arcade/patch/agent/process/PatchProcessSensingSimple.java +++ b/src/arcade/patch/agent/process/PatchProcessSensingSimple.java @@ -1,5 +1,6 @@ package arcade.patch.agent.process; +import java.util.logging.Logger; import ec.util.MersenneTwisterFast; import arcade.core.agent.process.Process; import arcade.core.sim.Simulation; @@ -13,6 +14,8 @@ * added to the environment at the rate specified by the input parameter VEGF_SECRETION_RATE. */ public class PatchProcessSensingSimple extends PatchProcessSensing { + private static final Logger LOGGER = + Logger.getLogger(PatchProcessSensingSimple.class.getName()); /** Rate of secretion of VEGF [VEGF/min]. */ private final double secretionRate; @@ -32,7 +35,6 @@ public class PatchProcessSensingSimple extends PatchProcessSensing { */ public PatchProcessSensingSimple(PatchCell cell) { super(cell); - secretionRate = cell.getParameters().getDouble("sensing/VEGF_SECRETION_RATE"); } diff --git a/src/arcade/patch/env/component/PatchComponentDegrade.java b/src/arcade/patch/env/component/PatchComponentDegrade.java index 35da75184..cfdc99be9 100644 --- a/src/arcade/patch/env/component/PatchComponentDegrade.java +++ b/src/arcade/patch/env/component/PatchComponentDegrade.java @@ -2,6 +2,7 @@ import java.util.ArrayList; import java.util.HashSet; +import java.util.logging.Logger; import sim.engine.Schedule; import sim.engine.SimState; import sim.util.Bag; @@ -31,6 +32,8 @@ * properties are recalculated. */ public class PatchComponentDegrade implements Component { + private static final Logger LOGGER = Logger.getLogger(PatchComponentDegrade.class.getName()); + /** Interval between degradation steps [min]. */ private final int degradationInterval; diff --git a/src/arcade/patch/env/component/PatchComponentGrowth.java b/src/arcade/patch/env/component/PatchComponentGrowth.java new file mode 100644 index 000000000..2e3507d07 --- /dev/null +++ b/src/arcade/patch/env/component/PatchComponentGrowth.java @@ -0,0 +1,1409 @@ +package arcade.patch.env.component; + +import java.util.ArrayList; +import java.util.EnumMap; +import java.util.HashMap; +import java.util.LinkedHashSet; +import java.util.Map; +import java.util.logging.Logger; +import sim.engine.Schedule; +import sim.engine.SimState; +import sim.util.Bag; +import ec.util.MersenneTwisterFast; +import arcade.core.env.component.Component; +import arcade.core.env.lattice.Lattice; +import arcade.core.sim.Series; +import arcade.core.sim.Simulation; +import arcade.core.util.Graph; +import arcade.core.util.MiniBox; +import arcade.core.util.Solver; +import arcade.core.util.Solver.Function; +import arcade.core.util.exceptions.IncompatibleFeatureException; +import arcade.core.util.exceptions.MissingSpecificationException; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteEdge; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteNode; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeDirection; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeLevel; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeType; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.Root; +import arcade.patch.env.location.CoordinateXYZ; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.CAPILLARY_RADIUS; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.MAXIMUM_CAPILLARY_RADIUS; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.MINIMUM_CAPILLARY_RADIUS; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.calculateCurrentState; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.calculateLocalFlow; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.calculatePressure; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.calculateThickness; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.getPath; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.path; +import static arcade.patch.util.PatchEnums.Ordering; + +/** + * Implementation of {@link Component} for adding graph edges and resolving mass balances. + * + *

This component can only be used with {@link PatchComponentsSitesGraph}. The component is + * stepped according to a reasonable interval based on the specified {@code MIGRATION_RATE}. + * Generally, if the average VEGF concentration in the immediate neighborhood of a node is greater + * than the threshold {@code VEGF_THRESHOLD} the node will be flagged for sprouting. The sprout + * direction is determined by the {@code WALK_TYPE} and the maximum length of migration is + * determined by the {@code MAX_LENGTH}. The walk types can bvalidateecified as {@code RANDOM}, + * {@code DETERMINISTIC}, or {@code BIASED}. + * + *

New edges are added based on the specified {@link CAPILLARY_RADIUS}. The calculation of the + * resulting hemodynamics is based on two strategies. The first is {@code COMPENSATE} which + * increases the flow of the original artery to compensate for the new edge. The second is {@code + * DIVERT} which calculates the diverted flow rate and pressure of the resulting edge by subtracting + * the flow rate of the new edge from the original edge. + */ +public class PatchComponentGrowth implements Component { + /** Logger for {@code PatchComponentGrowth}. */ + private static final Logger LOGGER = Logger.getLogger(PatchComponentGrowth.class.getName()); + + /** Default edge level to add to the graph from this component. */ + private static final EdgeLevel DEFAULT_EDGE_LEVEL = EdgeLevel.LEVEL_2; + + /** Default edge type to add to the graph from this component. */ + private static final EdgeType DEFAULT_EDGE_TYPE = EdgeType.ANGIOGENIC; + + /** + * Value used to determine how many steps a node should be delayed before being considered for + * growth after attempting growth. + */ + private static final int DELAY = 10; + + /** Calculation strategies. */ + public enum Calculation { + /** Code for upstream calculation strategy. */ + COMPENSATE, + + /** Code for downstream direction strategy. */ + DIVERT + } + + /** Strategy for direction of migration. */ + public enum MigrationDirection { + /** Code for random direction. */ + RANDOM, + + /** Code for deterministic direction. */ + DETERMINISTIC, + + /** Code for biased random direction. */ + BIASED; + } + + /** How to adjust the flow rate during radius calculations. */ + private enum Adjustment { + /** Code for increasing flow. */ + INCREASE, + + /** Code for decreasing flow. */ + DECREASE; + } + + private enum Outcome { + /** Code for successful calculation. */ + SUCCESS, + + /** Code for unsuccessful calculation. */ + FAILURE; + } + + /** Rate of migration. */ + private double migrationRate; + + /** Angiogenesis threshold for vegf concentration near SiteNode to initiate migration. */ + private double vegfThreshold; + + /** Lattice containing vegf concentration. */ + private Lattice vegfLattice; + + /** Direction of migration. */ + private MigrationDirection walkType; + + /** Strategy for calculation of boundary conditions. */ + private Calculation calculationStrategy; + + /** Maximum length of migration. */ + private double maxLength; + + /** Maximum number of edges to be added to the graph from a sprout, based on maxLength. */ + private int maxEdges; + + /** Interval at which to step the component, based on migration rate. */ + private int interval; + + /** The size of an edge, based on the grid size. */ + private double edgeSize; + + /** The associated {@link PatchComponentSitesGraph} object. */ + private PatchComponentSitesGraph sites; + + /** The {@link Graph} object representing the sites. */ + private Graph graph; + + /** Persistent map of angiographic edges, keyed by the node they originate from. */ + private HashMap> angiogenicNodeMap = new HashMap<>(); + + /** List of temporary edges to be added to the graph this step. */ + private ArrayList tempEdges = new ArrayList<>(); + + /** List of directions to be used in the migration. */ + private ArrayList offsetDirections; + + /** Map of offsets to be used in the migration. */ + private EnumMap offsets; + + /** Tick for the current step. */ + private int tick; + + /** List of nodes to be removed from the angiogenic node map this time step. */ + private ArrayList keyNodesToRemove = new ArrayList<>(); + + /** Map of nodes to a delay to prevent nodes from constantly trying to add the same edge. */ + private HashMap nodeDelays = new HashMap<>(); + + /** + * Creates a growth component for a {@link PatchComponentSitesGraph}. + * + *

Loaded parameters include: + * + *

    + *
  • MIGRATION_RATE = How quickly the tip cells migrate. + *
  • VEGF_THRESHOLD = The threshold for the average VEGF concentration surrounding a node. + *
  • WALK_TYPE = How the directional migration is to be performed. + *
  • MAX_LENGTH = The maximum length of migration. + *
+ * + * @param series {@link Series} object. + * @param parameters {@link MiniBox} object. + */ + public PatchComponentGrowth(Series series, MiniBox parameters) { + // Set loaded parameters. + migrationRate = parameters.getDouble("MIGRATION_RATE"); + vegfThreshold = parameters.getDouble("VEGF_THRESHOLD"); + walkType = MigrationDirection.valueOf(parameters.get("WALK_TYPE")); + maxLength = parameters.getDouble("MAX_LENGTH"); + calculationStrategy = Calculation.valueOf(parameters.get("CALCULATION_STRATEGY")); + } + + @Override + public void schedule(Schedule schedule) { + interval = migrationRate < edgeSize ? 60 : (int) (edgeSize / migrationRate * 60); + schedule.scheduleRepeating(this, Ordering.LAST_COMPONENT.ordinal() + 1, interval); + } + + @Override + public void register(Simulation sim, String componentID) { + Component component = sim.getComponent(componentID); + + // validate + if (!(component instanceof PatchComponentSitesGraph)) { + throw new IncompatibleFeatureException( + "Growth Component", component.getClass().getName(), "PatchComponentSitesGraph"); + } + + vegfLattice = sim.getLattice("VEGF"); + if (vegfLattice == null) { + throw new MissingSpecificationException("VEGF layer must be included."); + } + + sites = (PatchComponentSitesGraph) component; + graph = sites.graph; + + offsets = sites.graphFactory.getOffsets(); + offsetDirections = new ArrayList<>(offsets.keySet()); + offsetDirections.remove(EdgeDirection.UNDEFINED); + + edgeSize = sim.getSeries().ds; + maxEdges = (int) Math.floor(maxLength / edgeSize); + } + + /** + * Returns a set of valid nodes. A node is valid if it not associated with an ignored edge, is + * not a root, and has fewer than 3 degrees (effectively 2 degrees exactly). + * + * @return a set of valid nodes + */ + private LinkedHashSet getValidNodes() { + LinkedHashSet set = new LinkedHashSet<>(); + for (Object obj : graph.getAllEdges()) { + SiteEdge edge = (SiteEdge) obj; + if (edge.isIgnored) { + continue; + } + SiteNode from = edge.getFrom(); + SiteNode to = edge.getTo(); + + if ((graph.getDegree(from) == 2) && !from.isRoot) { + set.add(from); + } + if ((graph.getDegree(to) == 2) && !to.isRoot) { + set.add(to); + } + } + return set; + } + + /** + * Helper function to handle which version of walk to use for the direction of migration. + * + * @param random the random number generator + * @param node the node to start from + * @param spanMap the map of directions to their corresponding values + * @return the direction to sprout + */ + private EdgeDirection performWalk( + MersenneTwisterFast random, SiteNode node, EnumMap spanMap) { + switch (walkType) { + case RANDOM: + return performRandomWalk(random, node, spanMap); + case BIASED: + return performBiasedWalk(random, node, spanMap); + case DETERMINISTIC: + return performDeterministicWalk(random, node, spanMap); + default: + return performDeterministicWalk(random, node, spanMap); + } + } + + @Override + public void step(SimState simstate) { + tick = (int) simstate.schedule.getTime(); + MersenneTwisterFast random = simstate.random; + + LinkedHashSet validNodes = getValidNodes(); + for (SiteNode node : validNodes) { + if (checkNodeSkipStatus(node)) { + continue; + } + + EnumMap> vegfMap = + buildDirectionalSpanMap(vegfLattice, node); + + if (averageDirectionalMap(vegfMap) > vegfThreshold) { + Bag in = graph.getEdgesIn(node); + Bag out = graph.getEdgesOut(node); + SiteEdge inEdge = (SiteEdge) in.get(0); + vegfMap.remove(sites.graphFactory.getOppositeDirection(inEdge, DEFAULT_EDGE_LEVEL)); + + SiteEdge outEdge = (SiteEdge) out.get(0); + vegfMap.remove(sites.graphFactory.getDirection(outEdge, DEFAULT_EDGE_LEVEL)); + + if (vegfMap.isEmpty()) { + continue; + } + + angiogenicNodeMap.put(node, new ArrayList<>()); + + EnumMap vegfAverages = getDirectionalAverages(vegfMap); + + node.sproutDir = performWalk(random, node, vegfAverages); + } + } + + boolean addFlag = propogateEdges(); + if (addFlag) { + for (SiteNode sproutNode : angiogenicNodeMap.keySet()) { + if (keyNodesToRemove.contains(sproutNode)) { + continue; + } + if (sproutNode.anastomosis) { + + nodeDelays.putIfAbsent(sproutNode, 0); + + if (graph.getDegree(sproutNode) > 2) { + keyNodesToRemove.add(sproutNode); + continue; + } + + int leadingIndex = angiogenicNodeMap.get(sproutNode).size() - 1; + assert leadingIndex >= 0; + SiteNode finalNode = + angiogenicNodeMap.get(sproutNode).get(leadingIndex).getTo(); + SiteNode init; + SiteNode fin; + + if (!graph.contains(finalNode)) { + // Connecting two angiogenic nodes + SiteNode targetNode = findKeyNodeInMap(finalNode, sproutNode); + + if (targetNode == null) { + sproutNode.anastomosis = false; + keyNodesToRemove.add(sproutNode); + continue; + } + if (keyNodesToRemove.contains(targetNode)) { + sproutNode.anastomosis = false; + keyNodesToRemove.add(sproutNode); + continue; + } + keyNodesToRemove.add(sproutNode); + keyNodesToRemove.add(targetNode); + + // Connecting sprout to existing node + sproutNode = validateNodeObject(sproutNode); + targetNode = validateNodeObject(targetNode); + + if (sproutNode.pressure == 0 + || targetNode.pressure == 0 + || Double.isNaN(sproutNode.pressure) + || Double.isNaN(targetNode.pressure)) { + continue; + } + + if (sproutNode.pressure < targetNode.pressure) { + reverseAllEdges(sproutNode); + init = targetNode; + fin = sproutNode; + } else { + reverseAllEdges(targetNode); + init = sproutNode; + fin = targetNode; + } + angiogenicNodeMap.get(sproutNode).addAll(angiogenicNodeMap.get(targetNode)); + for (SiteEdge e : angiogenicNodeMap.get(sproutNode)) { + if (e.getTo().equals(finalNode)) { + e.setTo(finalNode); + } else if (e.getFrom().equals(finalNode)) { + e.setFrom(finalNode); + } + } + + } else { + keyNodesToRemove.add(sproutNode); + keyNodesToRemove.add(finalNode); + + if (finalNode.sproutDir != null + || finalNode.anastomosis + || finalNode.isRoot) { + continue; + } + + // Connecting sprout to existing node + sproutNode = validateNodeObject(sproutNode); + finalNode = validateNodeObject(finalNode); + + if (sproutNode.pressure == 0 + || finalNode.pressure == 0 + || Double.isNaN(sproutNode.pressure) + || Double.isNaN(finalNode.pressure)) { + continue; + } + + if (sproutNode.pressure < finalNode.pressure) { + reverseAllEdges(sproutNode); + init = finalNode; + fin = sproutNode; + } else { + init = sproutNode; + fin = finalNode; + } + } + + assert init.pressure != 0.0; + assert fin.pressure != 0.0; + assert init != fin; + + addAngiogenicEdges( + angiogenicNodeMap.get(sproutNode), init, fin, calculationStrategy); + + sproutNode.anastomosis = false; + } + } + calculateCurrentState(graph); + } + + for (SiteNode n : keyNodesToRemove) { + if (angiogenicNodeMap.containsKey(n)) { + angiogenicNodeMap.remove(n); + } + } + keyNodesToRemove.clear(); + } + + private SiteNode validateNodeObject(SiteNode node) { + // Connecting sprout to existing node + if (node.pressure == 0.0) { + node = (SiteNode) graph.lookup(node); + if (graph.getEdgesOut(node) != null && node.pressure == 0.0) { + node = ((SiteEdge) graph.getEdgesOut(node).get(0)).getFrom(); + return node; + } + if (graph.getEdgesIn(node) != null && node.pressure == 0.0) { + node = ((SiteEdge) graph.getEdgesIn(node).get(0)).getTo(); + return node; + } + } + return node; + } + + /** + * Propogates the edges from each of the nodes in the angiogenic node map. + * + *

A node is stepped according to the sprout direction of the node. If the added edge + * connects to another node in the temporary graph, it becomes anastomotic. If the node is + * anastomotic, the add flag is set to true. If the node is not anastomotic, before reaching the + * max length, or cannot be added to the graph for another reason, the node added to the removal + * queue. + * + * @return {@code true} if an angiogenic edge becomes perfused, {@code false} otherwise + */ + private boolean propogateEdges() { + boolean addFlag = false; + addTemporaryEdges(); + + for (SiteNode sprout : angiogenicNodeMap.keySet()) { + if (checkForIgnoredEdges(sprout)) { + keyNodesToRemove.add(sprout); + continue; + } + + ArrayList edgeList = angiogenicNodeMap.get(sprout); + SiteNode tipNode; + SiteEdge newEdge; + + if (!edgeList.isEmpty()) { + tipNode = edgeList.get(edgeList.size() - 1).getTo(); + } else { + tipNode = sprout; + } + + if (tick - tipNode.lastUpdate < migrationRate) { + continue; + } + + newEdge = createNewEdge(sprout.sproutDir, tipNode); + + if (edgeList.size() == maxEdges || newEdge == null) { + keyNodesToRemove.add(sprout); + } else { + edgeList.add(newEdge); + if (newEdge.isAnastomotic) { + sprout.anastomosis = true; + addFlag = true; + } + } + } + + removeTemporaryEdges(); + return addFlag; + } + + /** + * Checks if the node has any edges that are marked as ignored. + * + * @param node {@link SiteNode} object. + * @return {@code true} if the node has any connected ignored edges. + */ + private boolean checkForIgnoredEdges(SiteNode node) { + Bag in = graph.getEdgesIn(node); + Bag out = graph.getEdgesOut(node); + if (in != null) { + for (Object edge : in) { + SiteEdge inEdge = (SiteEdge) edge; + if (inEdge.isIgnored) { + return true; + } + } + } + if (out != null) { + for (Object edge : out) { + SiteEdge outEdge = (SiteEdge) edge; + if (outEdge.isIgnored) { + return true; + } + } + } + return false; + } + + /** + * Criteria for skipping a node during the migration checks. If a node is in the node delay map, + * it will increment the value and return true. + * + * @param node the node to check + * @return {@code true} if the node should be skipped, {@code false} otherwise + */ + private boolean checkNodeSkipStatus(SiteNode node) { + if (angiogenicNodeMap.keySet().contains(node)) { + return true; + } + // if the node has been added in the last 72 hours, it is not considered for growth. + if ((tick - node.addTime) < (72 * 60)) { + return true; + } + // If a node has been previously marked as angiogenic, it is delayed for DELAY steps. + if (nodeDelays.containsKey(node)) { + nodeDelays.put(node, nodeDelays.get(node) + 1); + if (nodeDelays.get(node) > DELAY) { + nodeDelays.remove(node); + return false; + } + return true; + } + return false; + } + + /** + * Private helper function for reversing all edges in the angiogenicNodeMap for a given node. + * + * @param node {@link SiteNode} key to reverse edges for + */ + private void reverseAllEdges(SiteNode node) { + for (SiteEdge edge : angiogenicNodeMap.get(node)) { + edge.reverse(); + } + } + + /** + * Private helper function for lazily searching the key node in the angiogenicNodeMap for a + * given target node and skip node. + * + * @param targetNode {@link SiteNode} the node to look for in the map list values + * @param skipNode {@link SiteNode} to ignore in the map + * @return {@link SiteNode} key for the targetNode object + */ + SiteNode findKeyNodeInMap(SiteNode targetNode, SiteNode skipNode) { + for (SiteNode keyNode : angiogenicNodeMap.keySet()) { + if (keyNode.equals(skipNode)) { + continue; + } + if (keyNode.equals(targetNode)) { + return keyNode; + } + if (edgeListcontains(angiogenicNodeMap.get(keyNode), targetNode)) { + return keyNode; + } + } + return null; + } + + /** + * Private helper function for checking if an edge list contains a given target node. + * + * @param edgeList {@link ArrayList} of {@link SiteEdge}s + * @param targetNode {@link SiteNode} to search for + * @return {@code true} if the edge list contains the target node + */ + private boolean edgeListcontains(ArrayList edgeList, SiteNode targetNode) { + for (SiteEdge edge : edgeList) { + if (edge.getTo().equals(targetNode)) { + return true; + } + } + return false; + } + + /** + * Creates a temporary version of the graph that contains all the proposed edges from the + * angiogenic node map. + */ + private void addTemporaryEdges() { + for (Map.Entry> entry : angiogenicNodeMap.entrySet()) { + ArrayList edgeList = entry.getValue(); + tempEdges.addAll(edgeList); + addEdgeList(edgeList); + } + } + + /** Removes the proposed edges from the temporary version of the graph. */ + private void removeTemporaryEdges() { + if (tempEdges.isEmpty()) { + return; + } + removeEdgeList(tempEdges); + tempEdges.clear(); + } + + /** + * Removes all edges in an edge list from the graph. + * + * @param edgeList {@link SiteEdge}s to remove + */ + private void removeEdgeList(ArrayList edgeList) { + for (SiteEdge edge : edgeList) { + graph.removeEdge(edge); + } + } + + /** + * Private helper function for building a map of VEGF concentration values for a given node. + * + * @param lattice the VEGF lattice object + * @param node the angiogenic node object + * @return map of {@link EdgeDirection} to {@link ArrayList} of double values from the span of + * the edge in that direction + */ + private EnumMap> buildDirectionalSpanMap( + Lattice lattice, SiteNode node) { + double[][][] field = lattice.getField(); + EnumMap> vegfMap = new EnumMap<>(EdgeDirection.class); + for (EdgeDirection dir : offsetDirections) { + SiteNode proposed = sites.graphFactory.offsetNode(node, dir, DEFAULT_EDGE_LEVEL); + if (sites.graphFactory.checkNode(proposed)) { + ArrayList span = sites.getSpan(node, proposed); + vegfMap.put(dir, new ArrayList<>()); + for (CoordinateXYZ coordinate : span) { + int i = coordinate.x; + int j = coordinate.y; + int k = coordinate.z; + vegfMap.get(dir).add(field[k][i][j]); + } + } else { + vegfMap.put(dir, new ArrayList<>(0)); + } + } + return vegfMap; + } + + /** + * Private helper function for choosing a sprout direction randomly on the from the node. + * + * @param random simulation random number generator + * @param node the angiogenic node object + * @return a random edge direction + */ + static EdgeDirection performRandomWalk( + MersenneTwisterFast random, SiteNode node, EnumMap spanMap) { + ArrayList possibleDirections = new ArrayList<>(spanMap.keySet()); + return possibleDirections.get(random.nextInt(possibleDirections.size())); + } + + /** + * Private helper function for choosing a sprout direction biased on the VEGF concentration + * around the node. + * + * @param random simulation random number generator + * @param node the angiogenic node object + * @param spanMap map of direction to their respective VEGF concentration value + * @return a biased random edge direction + */ + static EdgeDirection performBiasedWalk( + MersenneTwisterFast random, SiteNode node, EnumMap spanMap) { + EnumMap seqMap = normalizeDirectionalMap(spanMap); + double val = random.nextDouble(); + for (EdgeDirection dir : seqMap.keySet()) { + if (val <= seqMap.get(dir)) { + return dir; + } + } + assert false; + return EdgeDirection.UNDEFINED; + } + + /** + * Private helper function for choosing a sprout direction deterministically on the VEGF + * concentration around the node. + * + * @param random simulation random number generator + * @param node the angiogenic node object + * @param spanMap map of direction to their respective VEGF concentration value + * @return the edge direction with highest concentration + */ + static EdgeDirection performDeterministicWalk( + MersenneTwisterFast random, SiteNode node, EnumMap spanMap) { + return getMaxKey(spanMap); + } + + /** + * Private helper function for getting the average VEGF concentration values for a given map. + * + * @param map map of {@link EdgeDirection} to {@link ArrayList} of double values from the span + * of the edge in that direction + * @return map of {@link EdgeDirection} to the average VEGF concentration value across the span + * in that direction + */ + static EnumMap getDirectionalAverages( + EnumMap> map) { + EnumMap averageMap = new EnumMap<>(EdgeDirection.class); + for (EdgeDirection dir : map.keySet()) { + double sum = 0; + for (double value : map.get(dir)) { + sum += value; + } + if (map.get(dir).size() == 0) { + averageMap.put(dir, 0.0); + } else { + averageMap.put(dir, sum / map.get(dir).size()); + } + } + return averageMap; + } + + /** + * Private helper function for getting the maximum VEGF concentration value for a given map. + * + * @param map map of {@link EdgeDirection} to the average VEGF concentration value across the + * span in that direction + * @return direction with the highest concentration + */ + static EdgeDirection getMaxKey(EnumMap map) { + EdgeDirection maxDir = EdgeDirection.UNDEFINED; + double maxVal = 0; + for (EdgeDirection dir : map.keySet()) { + if (map.get(dir) > maxVal) { + maxDir = dir; + maxVal = map.get(dir); + } + } + assert maxDir != EdgeDirection.UNDEFINED; + return maxDir; + } + + /** + * Private helper function for normalizing a map of VEGF concentration values. + * + * @param map {@link EnumMap} of {@link EdgeDirection} to real values. + * @return {@link EnumMap} of {@link EdgeDirection} to normalized values (range between 0 and + * 1). + */ + static EnumMap normalizeDirectionalMap( + EnumMap map) { + EnumMap normalizedMap = new EnumMap<>(EdgeDirection.class); + double norm = sumMap(map); + double prev = 0; + for (EdgeDirection dir : map.keySet()) { + prev = prev + map.get(dir) / norm; + normalizedMap.put(dir, prev); + } + return normalizedMap; + } + + /** + * Private helper function for summing a map of VEGF concentration values. + * + * @param map {@link EnumMap} of {@link EdgeDirection} to double values. + * @return sum of the values in the map + */ + static double sumMap(EnumMap map) { + double sum = 0; + for (EdgeDirection dir : map.keySet()) { + sum += map.get(dir); + } + return sum; + } + + /** + * Private helper function for averaging a map of VEGF concentration values. + * + * @param map {@link EnumMap} of {@link EdgeDirection} to {@link ArrayList} of double values. + * @return {@code double} object. + */ + static double averageDirectionalMap(EnumMap> map) { + double sum = 0; + int count = 0; + for (EdgeDirection dir : map.keySet()) { + for (double value : map.get(dir)) { + sum += value; + count++; + } + } + return sum / count; + } + + /** + * Private helper function for adding an edge list to the graph. + * + * @param list list of angiogenic edges to add to the graph + * @param start the starting site node object + * @param end the ending site node object + * @param calc code for the type of calculation to perform + */ + private void addAngiogenicEdges( + ArrayList list, SiteNode start, SiteNode end, Calculation calc) { + + path(graph, end, start); + if (end.prev != null) { + return; + } + + double otherRadius; + Bag outEdges = graph.getEdgesOut(start); + if (outEdges != null) { + otherRadius = ((SiteEdge) outEdges.get(0)).radius; + } else { + return; + } + + Graph tempGraph = new Graph(); + for (SiteEdge e : list) { + tempGraph.addEdge(e); + } + + // update edges in the minimal path between start and end + ArrayList angioPath = getPath(tempGraph, start, end); + + for (SiteEdge edge : angioPath) { + if (graph.contains(edge)) { + return; + } + edge.getTo().addTime = tick; + edge.radius = + (otherRadius > CAPILLARY_RADIUS) + ? CAPILLARY_RADIUS + : calculateEvenSplitRadius((SiteEdge) outEdges.get(0)); + if (Double.isNaN(edge.radius)) { + return; + } + edge.wall = calculateThickness(edge); + edge.span = sites.getSpan(edge.getFrom(), edge.getTo()); + edge.transport.putIfAbsent("GLUCOSE", 0.); + edge.transport.putIfAbsent("OXYGEN", 0.); + edge.fraction.putIfAbsent("GLUCOSE", 1.); + edge.length = sites.graphFactory.getLength(edge, DEFAULT_EDGE_LEVEL); + edge.isPerfused = true; + } + + addEdgeList(angioPath); + + switch (calc) { + case COMPENSATE: + updateRootsAndRadii(angioPath, start, end); + break; + case DIVERT: + default: + SiteNode intersection = + (SiteNode) + graph.findDownstreamIntersection( + (SiteEdge) outEdges.get(0), (SiteEdge) angioPath.get(0)); + Outcome result = recalculateRadii(angioPath, start, end, intersection); + if (result == Outcome.FAILURE) { + removeEdgeList(angioPath); + } + break; + } + } + + // private void addNodesInEdgeList(ArrayList edgeList) { + // if (edgeList == null || edgeList.isEmpty()) { + // return; + // } + // for (SiteEdge edge : edgeList) { + // SiteNode to = edge.getTo(); + // addedNodes.add(to); + // } + // } + + /** + * Private helper function for calculating the new radius of two edges after splitting flow + * evenly. + * + * @param edge the original edge with specified radius + * @return new radius for the edges + */ + private static double calculateEvenSplitRadius(SiteEdge edge) { + double radius = edge.radius; + double length = edge.length; + double deltaP = edge.getFrom().pressure - edge.getTo().pressure; + double flow = calculateLocalFlow(radius, length, deltaP); + double newRadius; + + try { + newRadius = + Solver.bisection( + (double r) -> flow - 2 * calculateLocalFlow(r, length, deltaP), + 1E-6, + 5 * MAXIMUM_CAPILLARY_RADIUS, + 1E-6); + } catch (Exception e) { + return Double.NaN; + } + return newRadius; + } + + /** + * Private helper function for updating the roots and radii of an edge list. + * + * @param addedEdges list of edges that were added + * @param start the starting {@link SiteNode} + * @param end the ending {@link SiteNode} + */ + private void updateRootsAndRadii(ArrayList addedEdges, SiteNode start, SiteNode end) { + ArrayList oldRadii = new ArrayList<>(); + ArrayList updatedEdges = new ArrayList<>(); + Boolean failed = false; + + Bag edges = graph.getEdgesOut(start); + + if (edges == null) { + return; + } + if (edges.size() < 2) { + return; + } + + Double deltaP = start.pressure - end.pressure; + Double newFlow = calculateLocalFlow(CAPILLARY_RADIUS, addedEdges, deltaP); + + ArrayList arteries = sites.graphFactory.arteries; + ArrayList veins = sites.graphFactory.veins; + + ArrayList> pathsArteries = new ArrayList<>(); + for (Root artery : arteries) { + ArrayList path = getPath(graph, artery.node, start); + if (path == null) { + continue; + } + pathsArteries.add(path); + } + + Double arteryFlow = newFlow / arteries.size(); + + for (ArrayList path : pathsArteries) { + assert !path.get(0).getFrom().isRoot; + + updatedEdges.addAll(path); + for (SiteEdge e : path) { + oldRadii.add(e.radius); + } + + SiteEdge rootEdge = path.remove(0); + + if (calculateArteryRootRadius(rootEdge, arteryFlow, Adjustment.INCREASE) + == Outcome.FAILURE) { + failed = true; + break; + } + if (updateRadiiOfEdgeList(path, arteryFlow, Adjustment.INCREASE) == Outcome.FAILURE) { + failed = true; + break; + } + } + + if (failed) { + resetRadii(updatedEdges, oldRadii); + return; + } + + ArrayList> pathsVeins = new ArrayList<>(); + for (Root vein : veins) { + ArrayList path = getPath(graph, end, vein.node); + if (path == null) { + continue; + } + pathsVeins.add(path); + } + + Double veinFlow = newFlow / veins.size(); + + for (ArrayList path : pathsVeins) { + assert !path.get(0).getFrom().isRoot; + + SiteEdge rootEdge = path.remove(path.size() - 1); + oldRadii.add(rootEdge.radius); + updatedEdges.add(rootEdge); + if (calculateArteryRootRadius(rootEdge, veinFlow, Adjustment.INCREASE) + == Outcome.FAILURE) { + failed = true; + break; + } + if (updateRadiiOfEdgeList(path, veinFlow, Adjustment.DECREASE) == Outcome.FAILURE) { + failed = true; + break; + } + } + + if (failed) { + resetRadii(updatedEdges, oldRadii); + return; + } + } + + /** + * Private helper function for recalculating the radii of an edge list. + * + * @param ignoredEdges list of edges that should not be changed + * @param start starting {@link SiteNode} + * @param end ending {@link SiteNode} + * @param intersection the intersecting node object from the edges out of start + */ + private Outcome recalculateRadii( + ArrayList ignoredEdges, SiteNode start, SiteNode end, SiteNode intersection) { + + Bag edges = graph.getEdgesOut(start); + + assert edges != null; + assert edges.size() >= 2; + + edges = graph.getEdgesOut(start); + if (edges.size() != 2) { + return Outcome.FAILURE; + } + + Integer angioIndex = ignoredEdges.contains(edges.get(0)) ? 0 : 1; + Integer nonAngioIndex = angioIndex ^ 1; + double deltaP = start.pressure - end.pressure; + + if (start.pressure * end.pressure == 0 || Double.isNaN(deltaP)) { + return Outcome.FAILURE; + } + + Double divertedFlow = calculateLocalFlow(CAPILLARY_RADIUS, ignoredEdges, deltaP); + + for (SiteEdge e : ignoredEdges) { + e.flow = divertedFlow; + } + + Double originalFlow = ((SiteEdge) edges.get(nonAngioIndex)).flow; + if (divertedFlow > originalFlow) { + return Outcome.FAILURE; + } + + if (intersection != null) { + if (intersection.isRoot) { + return updateRadiusToRoot( + (SiteEdge) edges.get(angioIndex), + sites.graphFactory.veins.get(0).node, + divertedFlow, + Adjustment.INCREASE, + ignoredEdges); + } + + if (updateRadius( + (SiteEdge) edges.get(nonAngioIndex), + intersection, + divertedFlow, + Adjustment.DECREASE, + ignoredEdges) + == Outcome.FAILURE) { + return Outcome.FAILURE; + } + + if (updateRadius( + (SiteEdge) edges.get(angioIndex), + intersection, + divertedFlow, + Adjustment.INCREASE, + ignoredEdges) + == Outcome.FAILURE) { + return Outcome.FAILURE; + } + + } else { + for (Root vein : sites.graphFactory.veins) { + SiteNode boundary = vein.node; + path(graph, start, boundary); + if (boundary.prev != null + && ((SiteEdge) edges.get(angioIndex)).radius > MINIMUM_CAPILLARY_RADIUS) { + return updateRadiusToRoot( + (SiteEdge) edges.get(angioIndex), + sites.graphFactory.veins.get(0).node, + divertedFlow, + Adjustment.INCREASE, + ignoredEdges); + } + } + return Outcome.FAILURE; + } + return Outcome.SUCCESS; + } + + /** + * Private helper function for updating the radius of an edge. + * + * @param edge the {@link SiteEdge} + * @param intersection downstream intersection {@link SiteNode} + * @param flow flow adjustment through the edge + * @param adjustment code for flow change + * @param ignored list of {@link SiteEdge} to not update + * @return Outcome.FAILURE for unsuccessful calculation, 0 for successful calculation + */ + private Outcome updateRadius( + SiteEdge edge, + SiteNode intersection, + double flow, + Adjustment adjustment, + ArrayList ignored) { + ArrayList edgesToUpdate = getPath(graph, edge.getTo(), intersection); + + if (edgesToUpdate == null) { + return Outcome.FAILURE; + } + + edgesToUpdate.add(0, edge); + + return updateRadiiOfEdgeList(edgesToUpdate, flow, adjustment, ignored); + } + + /** + * Private helper function for updating the radii of an edge list without ignoring edges. + * + * @param edges list of edges to update + * @param flow new flow + * @param adjustment code for flow change + * @return Outcome.FAILURE for unsuccessful calculation, 0 for successful calculation + */ + private Outcome updateRadiiOfEdgeList( + ArrayList edges, double flow, Adjustment adjustment) { + return updateRadiiOfEdgeList(edges, flow, adjustment, new ArrayList<>()); + } + + /** + * Private helper function for updating the radii of an edge list. + * + * @param edges list of edges to update + * @param flow new flow + * @param adjustment code for flow change + * @param ignored list of {@link SiteEdge} to not update + * @return Outcome.FAILURE for unsuccessful calculation, 0 for successful calculation + */ + private Outcome updateRadiiOfEdgeList( + ArrayList edges, + double flow, + Adjustment adjustment, + ArrayList ignored) { + ArrayList oldRadii = new ArrayList<>(); + for (SiteEdge e : edges) { + oldRadii.add(e.radius); + if (ignored.contains(e)) { + continue; + } + if (calculateRadius(e, flow, adjustment) == Outcome.FAILURE) { + resetRadii(edges, oldRadii); + return Outcome.FAILURE; + } + } + return Outcome.SUCCESS; + } + + /** + * Private helper function for calculating the radius of an edge based on a change in flow. + * + * @param edge edge to update + * @param flow new flow + * @param adjustment code for flow change + * @return Outcome.SUCCESS for successful update, Outcome.FAILURE for failure + */ + private Outcome calculateRadius(SiteEdge edge, double flow, Adjustment adjustment) { + if (flow <= 0) { + return Outcome.FAILURE; + } + + double adjustedFlow = (adjustment == Adjustment.DECREASE) ? -1 * flow : flow; + double originalRadius = edge.radius; + double deltaP = edge.getFrom().pressure - edge.getTo().pressure; + + assert edge.getFrom().pressure > 0.0; + assert edge.getTo().pressure > 0.0; + + double originalFlow = calculateLocalFlow(originalRadius, edge.length, deltaP); + + // double originalFlow = edge.flow; + Function f = + (double r) -> + originalFlow + adjustedFlow - calculateLocalFlow(r, edge.length, deltaP); + double newRadius; + + try { + newRadius = Solver.bisection(f, 1E-6, 5 * MAXIMUM_CAPILLARY_RADIUS, 1E-6); + } catch (Exception e) { + return Outcome.FAILURE; + } + + if (newRadius == 1E-6) { + return Outcome.FAILURE; + } + edge.radius = newRadius; + return Outcome.SUCCESS; + } + + /** + * Private helper function for calculating the radius of an edge associated with a vein root. + * + * @param edge {@link SiteEdge} object. + * @param flow new flow + * @param adjustment code for flow change + * @return Outcome.SUCCESS for successful update, Outcome.FAILURE for failure + */ + private Outcome calculateVeinRootRadius(SiteEdge edge, double flow, Adjustment adjustment) { + assert flow >= 0; + + double adjustedFlow = (adjustment == Adjustment.DECREASE) ? -1 * flow : flow; + double originalRadius = edge.radius; + double deltaP = edge.getFrom().pressure - edge.getTo().pressure; + double originalFlow = calculateLocalFlow(originalRadius, edge.length, deltaP); + + Function f = + (double r) -> + originalFlow + + adjustedFlow + - calculateLocalFlow( + r, + edge.length, + edge.getFrom().pressure + - calculatePressure(r, edge.type.category)); + + double newRadius; + try { + newRadius = Solver.bisection(f, .5 * originalRadius, 1.5 * originalRadius); + } catch (Exception e) { + return Outcome.FAILURE; + } + + if (newRadius == .5 * originalRadius || Double.isNaN(newRadius)) { + return Outcome.FAILURE; + } + + edge.radius = newRadius; + edge.getTo().pressure = calculatePressure(newRadius, edge.type.category); + return Outcome.SUCCESS; + } + + /** + * Private helper function for calculating the radius of an edge when updating the source root. + * + * @param edge {@link SiteEdge} object. + * @param flow new flow + * @param adjustment code for flow change + * @return Outcome.SUCCESS for successful update, Outcome.FAILURE for failure + */ + private Outcome calculateArteryRootRadius(SiteEdge edge, double flow, Adjustment adjustment) { + assert flow >= 0; + + double adjustedFlow = (adjustment == Adjustment.DECREASE) ? -1 * flow : flow; + double originalRadius = edge.radius; + double deltaP = edge.getFrom().pressure - edge.getTo().pressure; + double originalFlow = calculateLocalFlow(originalRadius, edge.length, deltaP); + + Function f = + (double r) -> + originalFlow + + adjustedFlow + - calculateLocalFlow( + r, + edge.length, + calculatePressure(r, edge.type.category) + - edge.getTo().pressure); + + double newRadius; + try { + newRadius = Solver.bisection(f, .5 * originalRadius, 1.5 * originalRadius); + } catch (Exception e) { + return Outcome.FAILURE; + } + + if (newRadius == .5 * originalRadius + || Double.isNaN(newRadius) + || newRadius == 1.5 * originalRadius) { + return Outcome.FAILURE; + } + + edge.radius = newRadius; + edge.getFrom().pressure = calculatePressure(newRadius, edge.type.category); + return Outcome.SUCCESS; + } + + /** + * Private helper function for updating the radius of an edge to the sink root node. + * + * @param edge edge to update + * @param intersection downstream {@link SiteNode} intersection + * @param flow new flow + * @param adjustment code for flow change + * @param ignored list of {@link SiteEdge} to not update + * @return Outcome.SUCCESS for successful update, Outcome.FAILURE for failure + */ + private Outcome updateRadiusToRoot( + SiteEdge edge, + SiteNode intersection, + double flow, + Adjustment adjustment, + ArrayList ignored) { + + assert flow >= 0; + + ArrayList veins = sites.graphFactory.veins; + ArrayList oldRadii = new ArrayList<>(); + for (Root vein : veins) { + ArrayList path = getPath(graph, edge.getTo(), vein.node); + if (path == null) { + continue; + } + path.add(0, edge); + for (SiteEdge e : path) { + oldRadii.add(e.radius); + if (ignored.contains(e)) { + continue; + } + if (e.getTo().isRoot) { + if (calculateVeinRootRadius(e, flow, adjustment) == Outcome.FAILURE) { + resetRadii(path, oldRadii); + return Outcome.FAILURE; + } + } else { + if (calculateRadius(e, flow, adjustment) == Outcome.FAILURE) { + resetRadii(path, oldRadii); + return Outcome.FAILURE; + } + } + } + break; + } + return Outcome.SUCCESS; + } + + /** + * Private helper function for resetting the radii of an edge list. + * + * @param edges {@link ArrayList} of {@link SiteEdge} objects to be reset + * @param oldRadii {@link ArrayList} of radii + */ + private void resetRadii(ArrayList edges, ArrayList oldRadii) { + for (int i = 0; i < oldRadii.size(); i++) { + edges.get(i).radius = oldRadii.get(i); + } + } + + /** + * Private helper function for adding an edge list to the graph. + * + * @param list {@link ArrayList} of {@link SiteEdge}s to add to the graph + */ + private void addEdgeList(ArrayList list) { + for (SiteEdge edge : list) { + graph.addEdge(edge); + } + } + + /** + * Private helper function for creating a valid new edge. + * + * @param direction {@link EdgeDirection} describing the offset direction + * @param node the original {@link SiteNode} object + * @return the created {@link SiteEdge} object, or null if the edge would not be valid + */ + private SiteEdge createNewEdge(EdgeDirection direction, SiteNode node) { + SiteNode proposed = sites.graphFactory.offsetNode(node, direction, DEFAULT_EDGE_LEVEL); + proposed.lastUpdate = tick; + if (sites.graphFactory.checkNode(proposed) && graph.getDegree(node) < 3) { + SiteEdge edge; + if (graph.contains(proposed)) { + if (graph.getDegree(proposed) > 2 + || graph.getEdgesOut(proposed) == null + || graph.getEdgesIn(proposed) == null) { + return null; + } + // SiteNode existing = validateNodeObject(proposed); + SiteNode existing = (SiteNode) graph.lookup(proposed); + + assert existing != null; + + if (Double.isNaN(existing.pressure)) { + return null; + } + + edge = new SiteEdge(node, existing, DEFAULT_EDGE_TYPE, DEFAULT_EDGE_LEVEL); + edge.isAnastomotic = true; + return edge; + } + edge = new SiteEdge(node, proposed, DEFAULT_EDGE_TYPE, DEFAULT_EDGE_LEVEL); + + if (graph.contains(edge)) { + return null; + } + return edge; + } + return null; + } +} diff --git a/src/arcade/patch/env/component/PatchComponentRemodel.java b/src/arcade/patch/env/component/PatchComponentRemodel.java index b2e1b4f6d..0c6b172cb 100644 --- a/src/arcade/patch/env/component/PatchComponentRemodel.java +++ b/src/arcade/patch/env/component/PatchComponentRemodel.java @@ -1,5 +1,6 @@ package arcade.patch.env.component; +import java.util.logging.Logger; import sim.engine.Schedule; import sim.engine.SimState; import sim.util.Bag; @@ -13,6 +14,8 @@ import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.MAXIMUM_WALL_RADIUS_FRACTION; import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.MINIMUM_CAPILLARY_RADIUS; import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.MINIMUM_WALL_THICKNESS; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.calculateCurrentState; +import static arcade.patch.env.component.PatchComponentSitesGraphUtilities.updateGraph; import static arcade.patch.util.PatchEnums.Ordering; /** @@ -28,6 +31,8 @@ * All hemodynamic properties are recalculated at the end of the step. */ public class PatchComponentRemodel implements Component { + private static final Logger LOGGER = Logger.getLogger(PatchComponentRemodel.class.getName()); + /** Interval between remodeling steps [min]. */ private final int remodelingInterval; @@ -175,15 +180,9 @@ public void step(SimState state) { // If any edges are removed, update the graph edges that are ignored. // Otherwise, recalculate pressure, flow, and stresses. if (removed) { - PatchComponentSitesGraphUtilities.updateGraph(graph); + updateGraph(graph); } else { - PatchComponentSitesGraphUtilities.calculatePressures(graph); - boolean reversed = PatchComponentSitesGraphUtilities.reversePressures(graph); - if (reversed) { - PatchComponentSitesGraphUtilities.calculatePressures(graph); - } - PatchComponentSitesGraphUtilities.calculateFlows(graph); - PatchComponentSitesGraphUtilities.calculateStresses(graph); + calculateCurrentState(graph); } } diff --git a/src/arcade/patch/env/component/PatchComponentSitesGraph.java b/src/arcade/patch/env/component/PatchComponentSitesGraph.java index 771147dd6..83c1f457b 100644 --- a/src/arcade/patch/env/component/PatchComponentSitesGraph.java +++ b/src/arcade/patch/env/component/PatchComponentSitesGraph.java @@ -3,6 +3,7 @@ import java.util.ArrayList; import java.util.HashMap; import java.util.LinkedHashSet; +import java.util.logging.Logger; import sim.util.Bag; import ec.util.MersenneTwisterFast; import arcade.core.env.location.Location; @@ -13,6 +14,7 @@ import arcade.core.util.MiniBox; import arcade.core.util.Solver; import arcade.core.util.Solver.Function; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeDirection; import arcade.patch.env.location.CoordinateXYZ; import arcade.patch.sim.PatchSeries; import static arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeLevel; @@ -44,6 +46,9 @@ * {@code A} / {@code a} for an artery or {@code V} / {@code v} for a vein. */ public abstract class PatchComponentSitesGraph extends PatchComponentSites { + /** Logger for {@code PatchComponentSitesGraph}. */ + private static final Logger LOGGER = Logger.getLogger(PatchComponentSitesGraph.class.getName()); + /** Tolerance for difference in internal and external concentrations. */ private static final double DELTA_TOLERANCE = 1E-8; @@ -239,29 +244,11 @@ void simpleStep() { * @param random the random number generator */ void complexStep(MersenneTwisterFast random) { - Bag allEdges = new Bag(graph.getAllEdges()); - // Check if graph has become unconnected. - boolean isConnected = false; - for (Object obj : allEdges) { - SiteEdge edge = (SiteEdge) obj; - if (edge.getFrom().isRoot && !edge.isIgnored) { - isConnected = true; - break; - } - } - if (!isConnected) { - for (SiteLayer layer : layers) { - for (int k = 0; k < latticeHeight; k++) { - for (int i = 0; i < latticeLength; i++) { - for (int j = 0; j < latticeWidth; j++) { - layer.delta[k][i][j] = 0; - } - } - } - } + if (checkDisconnected()) { return; } + ; // Iterate through each molecule. for (SiteLayer layer : layers) { @@ -281,10 +268,11 @@ void complexStep(MersenneTwisterFast random) { } } - allEdges.shuffle(random); + Bag currentEdges = new Bag(graph.getAllEdges()); + currentEdges.shuffle(random); // Iterate through each edge in graph. - for (Object obj : allEdges) { + for (Object obj : currentEdges) { SiteEdge edge = (SiteEdge) obj; if (edge.isIgnored) { continue; @@ -329,6 +317,10 @@ void complexStep(MersenneTwisterFast random) { intConc = edge.fraction.get(layer.name) * concentration; // fmol/um^3 intConcNew = intConc; // fmol/um^3 extConcNew = extConc; // fmol/um^3 + + if (Double.isNaN(intConc) || Double.isNaN(extConc)) { + LOGGER.info(layer.name + ": NaN (1)"); + } } if (Math.abs(intConc - extConc) > DELTA_TOLERANCE) { @@ -346,6 +338,9 @@ void complexStep(MersenneTwisterFast random) { dmdt = pa * (intConcNew - extConcNew); extConcNew += dmdt / latticePatchVolume; } + if (Double.isNaN(intConcNew) || Double.isNaN(extConcNew)) { + LOGGER.info(layer.name + ": NaN (2)"); + } } // Update external concentrations. @@ -370,6 +365,9 @@ void complexStep(MersenneTwisterFast random) { if (layer.name.equalsIgnoreCase("OXYGEN")) { edge.transport.put(layer.name, (intConc - intConcNew) * edge.flow); } else { + if (Double.isNaN((intConc - intConcNew) / concentration)) { + LOGGER.info(layer.name + ": NaN (3)"); + } edge.transport.put(layer.name, (intConc - intConcNew) / concentration); } } @@ -377,6 +375,34 @@ void complexStep(MersenneTwisterFast random) { } } + private boolean checkDisconnected() { + Bag allEdges = new Bag(graph.getAllEdges()); + + // Check if graph has become unconnected. + boolean isConnected = false; + for (Object obj : allEdges) { + SiteEdge edge = (SiteEdge) obj; + if (edge.getFrom().isRoot && !edge.isIgnored) { + isConnected = true; + break; + } + } + if (!isConnected) { + for (SiteLayer layer : layers) { + for (int k = 0; k < latticeHeight; k++) { + for (int i = 0; i < latticeLength; i++) { + for (int j = 0; j < latticeWidth; j++) { + layer.delta[k][i][j] = 0; + } + } + } + } + return true; + } + + return false; + } + /** * Extension of {@link arcade.core.util.Graph.Node} for site nodes. * @@ -398,6 +424,20 @@ public static class SiteNode extends Node { /** Distance for Dijkstra's algorithm. */ int distance; + /** Tick for the last update during growth. */ + int lastUpdate; + + /** Tick for when the node was added to the graph. */ + int addTime; + + /** Direction of the angiogenic sprout. */ + EdgeDirection sproutDir; + + /** + * {@code true} if the angiogenic sprout is anastomotic/perfused, {@code false} otherwise. + */ + boolean anastomosis; + /** Parent node. */ SiteNode prev; @@ -415,7 +455,7 @@ public static class SiteNode extends Node { } @Override - public Node duplicate() { + public SiteNode duplicate() { return new SiteNode(x, y, z); } @@ -466,6 +506,9 @@ public static class SiteEdge extends Edge { /** {@code true} if edge is ignored, {@code false} otherwise. */ boolean isIgnored; + /** {@code true} if edge is anastomotic, {@code false} otherwise. */ + boolean isAnastomotic; + /** Edge type. */ final EdgeType type; @@ -513,8 +556,8 @@ public static class SiteEdge extends Edge { * @param type the edge type * @param level the graph resolution level */ - SiteEdge(Node from, Node to, EdgeType type, EdgeLevel level) { - super(from, to); + SiteEdge(SiteNode from, SiteNode to, EdgeType type, EdgeLevel level) { + super(from, to, type != EdgeType.ANGIOGENIC); this.type = type; this.level = level; isVisited = false; @@ -605,6 +648,22 @@ public double getCircum() { public double getFlow() { return flow; } + + public String getFraction() { + StringBuilder sb = new StringBuilder(); + for (String key : fraction.keySet()) { + sb.append(key + ":" + fraction.get(key) + ","); + } + return sb.toString(); + } + + public String getTransport() { + StringBuilder sb = new StringBuilder(); + for (String key : transport.keySet()) { + sb.append(key + ":" + transport.get(key) + ","); + } + return sb.toString(); + } } /** @@ -759,6 +818,10 @@ private ArrayList traverseEdge(SiteNode node, String code) { } } + // If the flow out is zero exit and return no children. + if (flowOut == 0) { + return children; + } // Assign new fractions. for (Object obj : out) { SiteEdge edge = (SiteEdge) obj; @@ -841,6 +904,9 @@ private ArrayList traverseNode(SiteNode node, String code) { node.oxygen = Solver.bisection(func, 0, MAX_OXYGEN_PARTIAL_PRESSURE); } + assert node.oxygen >= 0; + assert !Double.isNaN(node.oxygen); + // Recurse through output edges. for (Object obj : out) { SiteEdge edge = (SiteEdge) obj; diff --git a/src/arcade/patch/env/component/PatchComponentSitesGraphFactory.java b/src/arcade/patch/env/component/PatchComponentSitesGraphFactory.java index e9ccfb69c..87b42fe73 100644 --- a/src/arcade/patch/env/component/PatchComponentSitesGraphFactory.java +++ b/src/arcade/patch/env/component/PatchComponentSitesGraphFactory.java @@ -1,6 +1,7 @@ package arcade.patch.env.component; import java.util.ArrayList; +import java.util.EnumMap; import java.util.HashMap; import java.util.LinkedHashSet; import java.util.regex.Matcher; @@ -56,7 +57,10 @@ public enum EdgeType { VEIN(EdgeCategory.VEIN), /** Code for venule edge type. */ - VENULE(EdgeCategory.VEIN); + VENULE(EdgeCategory.VEIN), + + /** Code for angiogenic edge type. */ + ANGIOGENIC(EdgeCategory.CAPILLARY); /** Edge category corresponding to the edge type. */ final EdgeCategory category; @@ -188,6 +192,12 @@ enum EdgeMotif { /** Width of the array (y direction). */ final int latticeWidth; + /** List of pointers to artery node and edge objects. */ + ArrayList arteries; + + /** List of pointers to vein node and edge objects. */ + ArrayList veins; + /** * Creates a factory for making {@link Graph} sites. * @@ -245,6 +255,17 @@ public PatchComponentSitesGraphFactory(Series series) { */ abstract EdgeDirection getDirection(int fromX, int fromY, int toX, int toY); + /** + * Gets the opposite direction of the given edge. + * + * @param edge the edge object + * @param level the graph resolution level + * @return the code for the opposite edge direction + */ + public EdgeDirection getOppositeDirection(SiteEdge edge, EdgeLevel level) { + return getDirection(edge.getTo(), edge.getFrom(), level); + } + /** * Adds a root motif to the graph. * @@ -327,6 +348,13 @@ abstract void addConnection( */ abstract int[] getOffset(EdgeDirection offset); + /** + * Get a map of possible offset directions to their corresponding coordinate changes. + * + * @return the map of offset directions to their corresponding coordinate changes + */ + abstract EnumMap getOffsets(); + /** * Gets the length of the given edge. * @@ -604,8 +632,8 @@ public Graph initializeRootGraph(MersenneTwisterFast random, String graphLayout) leaves.addAll(bag); } - ArrayList arteries = new ArrayList<>(); - ArrayList veins = new ArrayList<>(); + arteries = new ArrayList<>(); + veins = new ArrayList<>(); boolean hasArtery = false; boolean hasVein = false; @@ -643,19 +671,19 @@ public Graph initializeRootGraph(MersenneTwisterFast random, String graphLayout) addMotifs(graph, leaves2, EdgeLevel.LEVEL_1, EdgeMotif.SINGLE, random); // Calculate radii, pressure, and shears. - updateRootGraph(graph, arteries, veins, EdgeLevel.LEVEL_1, random); + updateRootGraph(graph, EdgeLevel.LEVEL_1, random); // Iterative remodeling. int iter = 0; double frac = 1.0; while (frac > REMODELING_FRACTION && iter < MAX_ITERATIONS) { frac = remodelRootGraph(graph, EdgeLevel.LEVEL_1, random); - updateRootGraph(graph, arteries, veins, EdgeLevel.LEVEL_1, random); + updateRootGraph(graph, EdgeLevel.LEVEL_1, random); iter++; } // Prune network for perfused segments and recalculate properties. - refineRootGraph(graph, arteries, veins); + refineRootGraph(graph); // Subdivide growth sites and add new motifs. Bag midpoints = subdivideRootGraph(graph, EdgeLevel.LEVEL_1); @@ -664,10 +692,10 @@ public Graph initializeRootGraph(MersenneTwisterFast random, String graphLayout) addMotifs(graph, midpoints2, EdgeLevel.LEVEL_2, EdgeMotif.SINGLE, random); // Calculate radii, pressure, and shears. - updateRootGraph(graph, arteries, veins, EdgeLevel.LEVEL_2, random); + updateRootGraph(graph, EdgeLevel.LEVEL_2, random); // Prune network for perfused segments and recalculate properties. - refineRootGraph(graph, arteries, veins); + refineRootGraph(graph); return graph; } @@ -676,17 +704,10 @@ public Graph initializeRootGraph(MersenneTwisterFast random, String graphLayout) * Updates hemodynamic properties for graph sites with root layouts. * * @param graph the graph instance - * @param arteries the list of artery edges - * @param veins the list of vein edges * @param level the graph resolution level * @param random the random number generator */ - private void updateRootGraph( - Graph graph, - ArrayList arteries, - ArrayList veins, - EdgeLevel level, - MersenneTwisterFast random) { + private void updateRootGraph(Graph graph, EdgeLevel level, MersenneTwisterFast random) { ArrayList list; ArrayList caps = new ArrayList<>(); @@ -762,7 +783,7 @@ private void updateRootGraph( Graph g2 = new Graph(); graph.getSubgraph(g1, e -> ((SiteEdge) e).level == EdgeLevel.LEVEL_1); graph.getSubgraph(g2, e -> ((SiteEdge) e).level == EdgeLevel.LEVEL_2); - mergeGraphs(g1, g2); + mergeGraphs(graph, g1, g2); break; default: break; @@ -811,10 +832,8 @@ private void updateRootGraph( * Refines the graph for graph sites with root layouts. * * @param graph the graph instance - * @param arteries the list of artery edges - * @param veins the list of vein edges */ - private void refineRootGraph(Graph graph, ArrayList arteries, ArrayList veins) { + private void refineRootGraph(Graph graph) { // Reverse edges that are veins and venules. ArrayList reverse = getEdgeByType(graph, new EdgeType[] {EdgeType.VEIN, EdgeType.VENULE}); diff --git a/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryRect.java b/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryRect.java index 75aa944eb..98d3722c2 100644 --- a/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryRect.java +++ b/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryRect.java @@ -155,6 +155,11 @@ int[] getOffset(EdgeDirection offset) { return OFFSETS.get(offset); } + @Override + EnumMap getOffsets() { + return OFFSETS; + } + /** * Checks if there is an edge in the cross diagonal. * @@ -446,7 +451,8 @@ && checkCross(graph, node0, node1, level)) { return bag; } - // Add the two leaves of the tripod if line is 0, otherwise add in the root line. + // Add the two leaves of the tripod if line is 0, otherwise add in the root + // line. if (offsets == null) { for (EdgeDirection offset : ROOT_OFFSETS.get(dir)) { SiteNode node2 = offsetNode(node1, offset, level); diff --git a/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryTri.java b/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryTri.java index 4935ec12a..6b6d0a750 100644 --- a/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryTri.java +++ b/src/arcade/patch/env/component/PatchComponentSitesGraphFactoryTri.java @@ -151,6 +151,11 @@ int[] getOffset(EdgeDirection offset) { return OFFSETS.get(offset); } + @Override + EnumMap getOffsets() { + return OFFSETS; + } + @Override int calcOffset(int k) { return (latticeHeight - k / 2 - 1) % 3; @@ -433,7 +438,8 @@ Bag addRoot( return bag; } - // Add the two leaves of the tripod if line is 0, otherwise add in the root line. + // Add the two leaves of the tripod if line is 0, otherwise add in the root + // line. if (offsets == null) { for (EdgeDirection offset : ROOT_OFFSETS.get(dir)) { SiteNode node2 = offsetNode(node1, offset, level); diff --git a/src/arcade/patch/env/component/PatchComponentSitesGraphUtilities.java b/src/arcade/patch/env/component/PatchComponentSitesGraphUtilities.java index 63c6bb050..134be2a31 100644 --- a/src/arcade/patch/env/component/PatchComponentSitesGraphUtilities.java +++ b/src/arcade/patch/env/component/PatchComponentSitesGraphUtilities.java @@ -1,8 +1,10 @@ package arcade.patch.env.component; import java.util.ArrayList; +import java.util.Collections; import java.util.HashSet; import java.util.LinkedHashSet; +import java.util.logging.Logger; import sim.util.Bag; import ec.util.MersenneTwisterFast; import arcade.core.util.Graph; @@ -19,6 +21,9 @@ /** Container for utility functions used by {@link PatchComponentSitesGraph}. */ abstract class PatchComponentSitesGraphUtilities { + private static final Logger LOGGER = + Logger.getLogger(PatchComponentSitesGraphUtilities.class.getName()); + /** Calculation types. */ enum CalculationType { /** Code for upstream radius calculation for all edge types. */ @@ -155,8 +160,19 @@ private static double calculateViscosity(double radius) { * @return the flow rate coefficient */ private static double getCoefficient(SiteEdge edge) { - double mu = PLASMA_VISCOSITY * calculateViscosity(edge.radius) / 60; - return (Math.PI * Math.pow(edge.radius, 4)) / (8 * mu * edge.length); + return getCoefficient(edge.radius, edge.length); + } + + /** + * Gets flow rate coefficient in units of um3/(mmHg min). + * + * @param radius the edge radius + * @param length the edge length + * @return the flow rate coefficient + */ + private static double getCoefficient(double radius, double length) { + double mu = PLASMA_VISCOSITY * calculateViscosity(radius) / 60; + return (Math.PI * Math.pow(radius, 4)) / (8 * mu * length); } /** @@ -348,10 +364,11 @@ static void checkPerfused(Graph graph, ArrayList arteries, ArrayList /** * Merges the nodes from one graph with another graph. * - * @param graph1 the first graph object - * @param graph2 the second graph object + * @param graph the original graph object + * @param graph1 the first subgraph object + * @param graph2 the second subgraph object */ - static void mergeGraphs(Graph graph1, Graph graph2) { + static void mergeGraphs(Graph graph, Graph graph1, Graph graph2) { // Merge nodes for subgraph. graph2.mergeNodes(); @@ -374,6 +391,7 @@ static void mergeGraphs(Graph graph1, Graph graph2) { } } } + graph.combine(graph1, graph2); } /** @@ -473,6 +491,10 @@ static void calculatePressures(Graph graph) { if (div != 0) { x0[id] /= div; } + + if (node.pressure > 0) { + x0[id] = node.pressure; + } } double[][] sA = Matrix.scale(mA, 1E-7); @@ -526,6 +548,36 @@ static void calculateStresses(Graph graph) { } } + /** + * Calculate the the flow rate for a given set of edges without any branches. + * + * @param radius the radius of the edges + * @param edges the list of edges + * @param deltaP the pressure change + * @return the flow rate (in um3/min) + */ + static double calculateLocalFlow(double radius, ArrayList edges, double deltaP) { + double length = 0; + + for (SiteEdge edge : edges) { + length += edge.length; + } + + return getCoefficient(radius, length) * (deltaP); + } + + /** + * Calculate the the flow rate for a given edge without any branches. + * + * @param radius the radius of the edge + * @param length the length of the edge + * @param deltaP the pressure change + * @return the flow rate (in um3/min) + */ + static double calculateLocalFlow(double radius, double length, double deltaP) { + return getCoefficient(radius, length) * deltaP; + } + /** * Calculates flow rate (in um3/min) and area (in um2) for all edges. * @@ -559,11 +611,21 @@ static void calculateFlows(Graph graph) { static void calculateThicknesses(Graph graph) { for (Object obj : graph.getAllEdges()) { SiteEdge edge = (SiteEdge) obj; - double d = 2 * edge.radius; - edge.wall = d * (0.267 - 0.084 * Math.log10(d)); + edge.wall = calculateThickness(edge); } } + /** + * Calculates the thickness of an edge. + * + * @param edge the edge object + * @return the thickness of the edge + */ + static double calculateThickness(SiteEdge edge) { + double d = 2 * edge.radius; + return d * (0.267 - 0.084 * Math.log10(d)); + } + /** * Gets in degree for edge in given calculation direction. * @@ -892,7 +954,8 @@ static void path(Graph graph, SiteNode start, SiteNode end) { settled.add(evalNode); // If end node found, exit from loop. - if (evalNode == end) { + if (evalNode.equals(end)) { + // end.prev = evalNode.prev; break; } @@ -923,6 +986,55 @@ static void path(Graph graph, SiteNode start, SiteNode end) { } } + /** + * Get the path between two nodes in the graph. + * + * @param graph the graph object + * @param start the start node + * @param end the end node + * @return the list of edges in the path + */ + static ArrayList getPath(Graph graph, SiteNode start, SiteNode end) { + path(graph, start, end); + ArrayList path = new ArrayList<>(); + SiteNode node = end; + if (node.prev == null) { + node = (SiteNode) graph.lookup(end); + } + while (node != null && !node.equals(start)) { + Bag b = graph.getEdgesIn(node); + if (b == null) { + LOGGER.info("NODETOINBAG IS NULL IN GETPATH FOR NODE " + node.toString()); + } + if (b.numObjs == 1) { + path.add((SiteEdge) b.objs[0]); + } else if (b.numObjs == 2) { + SiteEdge edgeA = ((SiteEdge) b.objs[0]); + SiteEdge edgeB = ((SiteEdge) b.objs[1]); + if (edgeA.getFrom().equals(node.prev)) { + path.add(edgeA); + } else { + path.add(edgeB); + } + } + + // if (node.prev == null) { + // LOGGER.info("START: " + start + " END: " + end); + // LOGGER.info("PREV IS NULL" + node); + // } + + node = node.prev; + } + + if (node == null) { + // LOGGER.info("Path in getPath is: " + path); + return null; + } + + Collections.reverse(path); + return path; + } + /** * Traverses through graph to find perfused paths. * @@ -1099,6 +1211,42 @@ static void updateRadii( * @param graph the graph object */ static void updateGraph(Graph graph) { + trimGraph(graph); + + calculateCurrentState(graph); + + // Set oxygen nodes. + for (Object obj : graph.getAllEdges()) { + SiteEdge edge = (SiteEdge) obj; + SiteNode to = edge.getTo(); + SiteNode from = edge.getFrom(); + if (Double.isNaN(to.pressure)) { + to.oxygen = Double.NaN; + } + if (Double.isNaN(from.pressure)) { + from.oxygen = Double.NaN; + } + } + } + + /** + * Updates hemodynamic properties based on the current state of the graph. + * + * @param graph the graph object + */ + static void calculateCurrentState(Graph graph) { + do { + calculatePressures(graph); + boolean reversed = reversePressures(graph); + if (reversed) { + calculatePressures(graph); + } + calculateFlows(graph); + calculateStresses(graph); + } while (checkForNegativeFlow(graph)); + } + + static void trimGraph(Graph graph) { ArrayList list; Graph gCurr = graph; @@ -1133,27 +1281,21 @@ static void updateGraph(Graph graph) { gCurr = gNew; } while (list.size() != 0); + } - calculatePressures(graph); - boolean reversed = reversePressures(graph); - if (reversed) { - calculatePressures(graph); - } - calculateFlows(graph); - calculateStresses(graph); - - // Set oxygen nodes. + // This *MIGHT* be a problem, I think we could revisit adding this check. + // I'm not sure why it would get to the point where there would be a negative flow in the graph? + static boolean checkForNegativeFlow(Graph graph) { + boolean negative = false; for (Object obj : graph.getAllEdges()) { SiteEdge edge = (SiteEdge) obj; - SiteNode to = edge.getTo(); - SiteNode from = edge.getFrom(); - if (Double.isNaN(to.pressure)) { - to.oxygen = Double.NaN; - } - if (Double.isNaN(from.pressure)) { - from.oxygen = Double.NaN; + if (edge.flow < 0) { + negative = true; + LOGGER.info("Negative flow detected, recalculating."); + break; } } + return negative; } /** @@ -1176,6 +1318,8 @@ static void updateTraverse(Graph graph, LinkedHashSet nodes, boolean r for (Object obj : out) { SiteEdge edge = (SiteEdge) obj; if (edge.flow < MINIMUM_FLOW_RATE || Double.isNaN(edge.flow)) { + LOGGER.info("Removing Edge. 1309, edge: " + edge); + LOGGER.info("Flow: " + edge.flow); graph.removeEdge(edge); edge.getFrom().pressure = Double.NaN; edge.getTo().pressure = Double.NaN; @@ -1191,6 +1335,8 @@ static void updateTraverse(Graph graph, LinkedHashSet nodes, boolean r for (Object obj : in) { SiteEdge edge = (SiteEdge) obj; if (edge.flow < MINIMUM_FLOW_RATE || Double.isNaN(edge.flow)) { + LOGGER.info("Removing Edge. 1325, edge: " + edge); + LOGGER.info("Flow: " + edge.flow); graph.removeEdge(edge); edge.getFrom().pressure = Double.NaN; edge.getTo().pressure = Double.NaN; @@ -1208,11 +1354,13 @@ static void updateTraverse(Graph graph, LinkedHashSet nodes, boolean r double totalFlow = edge1.flow + edge2.flow; if (edge1.flow / totalFlow < MINIMUM_FLOW_PERCENT) { + LOGGER.info("Removing Edge. 1343, edge: " + edge1); graph.removeEdge(edge1); edge1.getFrom().pressure = Double.NaN; edge1.getTo().pressure = Double.NaN; updateGraph(graph); } else if (edge2.flow / totalFlow < MINIMUM_FLOW_PERCENT) { + LOGGER.info("Removing Edge. 1349, edge: " + edge2); graph.removeEdge(edge2); edge2.getFrom().pressure = Double.NaN; edge2.getTo().pressure = Double.NaN; @@ -1223,6 +1371,7 @@ static void updateTraverse(Graph graph, LinkedHashSet nodes, boolean r } if (removeMin) { + LOGGER.info("Removing Edge. 1360, edge: " + minEdge); graph.removeEdge(minEdge); minEdge.getFrom().pressure = Double.NaN; minEdge.getTo().pressure = Double.NaN; diff --git a/src/arcade/patch/parameter.patch.xml b/src/arcade/patch/parameter.patch.xml index 40824712c..75876148b 100644 --- a/src/arcade/patch/parameter.patch.xml +++ b/src/arcade/patch/parameter.patch.xml @@ -94,7 +94,7 @@ - + @@ -155,6 +155,7 @@ + @@ -162,4 +163,11 @@ + + + + + + + diff --git a/src/arcade/patch/sim/PatchSimulationHex.java b/src/arcade/patch/sim/PatchSimulationHex.java index 6eddf28b5..24da8a3cd 100644 --- a/src/arcade/patch/sim/PatchSimulationHex.java +++ b/src/arcade/patch/sim/PatchSimulationHex.java @@ -10,6 +10,7 @@ import arcade.patch.agent.cell.PatchCellFactory; import arcade.patch.env.component.PatchComponentCycle; import arcade.patch.env.component.PatchComponentDegrade; +import arcade.patch.env.component.PatchComponentGrowth; import arcade.patch.env.component.PatchComponentPulse; import arcade.patch.env.component.PatchComponentRemodel; import arcade.patch.env.component.PatchComponentSitesGraphTri; @@ -82,6 +83,8 @@ public Component makeComponent(String componentClass, MiniBox parameters) { return new PatchComponentDegrade(series, parameters); case "remodel": return new PatchComponentRemodel(series, parameters); + case "growth": + return new PatchComponentGrowth(series, parameters); default: return null; } diff --git a/src/arcade/patch/sim/output/PatchOutputSerializer.java b/src/arcade/patch/sim/output/PatchOutputSerializer.java index cdfa02efa..44996416f 100644 --- a/src/arcade/patch/sim/output/PatchOutputSerializer.java +++ b/src/arcade/patch/sim/output/PatchOutputSerializer.java @@ -244,6 +244,8 @@ public JsonElement serialize( json.addProperty("shear", src.getShear()); json.addProperty("stress", src.getCircum()); json.addProperty("flow", src.getFlow()); + json.addProperty("fraction", src.getFraction()); + json.addProperty("transport", src.getTransport()); return json; } diff --git a/test/arcade/core/util/GraphTest.java b/test/arcade/core/util/GraphTest.java index 683577224..984f2ef45 100644 --- a/test/arcade/core/util/GraphTest.java +++ b/test/arcade/core/util/GraphTest.java @@ -797,7 +797,8 @@ public void update_updatesGraph() { () -> assertTrue(graph.contains(edge3)), () -> assertTrue(graph.contains(edge4)), () -> assertTrue(graph.contains(edge5)), - () -> assertTrue(graph.contains(edge6))); + () -> assertTrue(graph.contains(edge6)), + () -> assertEquals(graph.lookup(node3B), graph.lookup(new Node(-2, -2, 0)))); } @Test @@ -861,4 +862,66 @@ public void graph_toString_returnsString() { assertEquals(expected, graph.toString()); } + + @Test + public void lookup_returnsNode() { + Graph graph = new Graph(); + Node node0 = new Node(-1, 0, 0); + Node node1 = new Node(0, 0, 0); + Node node2 = new Node(3, 0, 0); + + Edge edge0 = new Edge(node0, node1); + Edge edge1 = new Edge(node1, node2); + + graph.addEdge(edge0); + graph.addEdge(edge1); + + Node node0Copy = new Node(-1, 0, 0); + Node node1Copy = new Node(0, 0, 0); + Node node2Copy = new Node(3, 0, 0); + + assertAll( + () -> assertFalse(node0Copy == graph.lookup(-1, 0, 0)), + () -> assertFalse(node1Copy == graph.lookup(0, 0, 0)), + () -> assertFalse(node2Copy == graph.lookup(3, 0, 0)), + () -> assertTrue(graph.lookup(node0Copy) == graph.lookup(-1, 0, 0)), + () -> assertTrue(graph.lookup(node1Copy) == graph.lookup(0, 0, 0)), + () -> assertTrue(graph.lookup(node2Copy) == graph.lookup(3, 0, 0))); + } + + @Test + public void lookup_returnsNull() { + Graph graph = new Graph(); + Node node0 = new Node(-1, 0, 0); + Node node1 = new Node(0, 0, 0); + Node node2 = new Node(3, 0, 0); + + Edge edge0 = new Edge(node0, node1); + Edge edge1 = new Edge(node1, node2); + + graph.addEdge(edge0); + graph.addEdge(edge1); + + assertAll( + () -> assertNull(graph.lookup(-2, 0, 0)), () -> assertNull(graph.lookup(4, 0, 0))); + } + + @Test + public void lookup_returnsNullAfterEdgeIsRemoved() { + Graph graph = new Graph(); + Node node0 = new Node(-1, 0, 0); + Node node1 = new Node(0, 0, 0); + Node node2 = new Node(3, 0, 0); + + Edge edge0 = new Edge(node0, node1); + Edge edge1 = new Edge(node1, node2); + + graph.addEdge(edge0); + graph.addEdge(edge1); + graph.removeEdge(edge1); + + assertAll( + () -> assertNotNull(graph.lookup(0, 0, 0)), + () -> assertNull(graph.lookup(3, 0, 0))); + } } diff --git a/test/arcade/core/util/SolverTest.java b/test/arcade/core/util/SolverTest.java index 81dcd7ef5..da07ea396 100644 --- a/test/arcade/core/util/SolverTest.java +++ b/test/arcade/core/util/SolverTest.java @@ -236,4 +236,21 @@ public void testBisection_quadraticFunctionAndSwappedInputs_returnsAnswer() { assertEquals(1.41421, result, 0.0001); } + + @Test + public void testBisection_exceedsMaxIterations_returnsNaN() { + Function f = (x) -> x * x - 2; + double result = Solver.bisection(f, 2, 0, 2); + + assertEquals(Double.NaN, result, 0.001); + } + + @Test + public void testBisection_givenPrecision_returnsAnswersWithinPrecision() { + Function f = (x) -> x * x - 2; + double result = Solver.bisection(f, 2, 0, 1E-3); + + assertEquals(1.41421, result, 1E-3); + assertNotEquals(1.41421, result, 1E-4); + } } diff --git a/test/arcade/patch/env/component/PatchComponentGrowthTest.java b/test/arcade/patch/env/component/PatchComponentGrowthTest.java new file mode 100644 index 000000000..e6a7575bf --- /dev/null +++ b/test/arcade/patch/env/component/PatchComponentGrowthTest.java @@ -0,0 +1,319 @@ +package arcade.patch.env.component; + +import java.lang.reflect.Field; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.EnumMap; +import java.util.HashMap; +import org.junit.jupiter.api.BeforeAll; +import org.junit.jupiter.api.Test; +import sim.engine.Schedule; +import ec.util.MersenneTwisterFast; +import arcade.core.util.MiniBox; +import arcade.core.util.Parameters; +import arcade.core.util.exceptions.IncompatibleFeatureException; +import arcade.core.util.exceptions.MissingSpecificationException; +import arcade.patch.agent.process.PatchProcessMetabolism; +import arcade.patch.agent.process.PatchProcessSignaling; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteEdge; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteNode; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeDirection; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeLevel; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeType; +import arcade.patch.env.grid.PatchGrid; +import arcade.patch.env.lattice.PatchLattice; +import arcade.patch.sim.PatchSeries; +import arcade.patch.sim.PatchSimulation; +import static org.junit.jupiter.api.Assertions.*; +import static org.mockito.Mockito.*; + +public class PatchComponentGrowthTest { + static final double EPSILON = 1E-8; + + static PatchSeries seriesMock; + + static PatchSimulation simMock; + + static Parameters parametersMock; + + static PatchProcessMetabolism metabolismMock; + + static PatchProcessSignaling signalingMock; + + static PatchGrid gridMock; + + static MersenneTwisterFast randomMock; + + static PatchLattice latticeMock; + + @BeforeAll + public static void setup() { + seriesMock = mock(PatchSeries.class, CALLS_REAL_METHODS); + simMock = mock(PatchSimulation.class); + doReturn(seriesMock).when(simMock).getSeries(); + parametersMock = spy(new Parameters(new MiniBox(), null, null)); + metabolismMock = mock(PatchProcessMetabolism.class); + signalingMock = mock(PatchProcessSignaling.class); + gridMock = mock(PatchGrid.class); + latticeMock = mock(PatchLattice.class); + doReturn(gridMock).when(simMock).getGrid(); + randomMock = mock(MersenneTwisterFast.class); + simMock.random = randomMock; + } + + @Test + public void schedule_calledWithMigrationRateLowerThanEdgeSize_setsScheduleCorrectly() { + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 10); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + Schedule schedule = mock(Schedule.class); + + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + + int ds = 30; + + try { + Field field = PatchComponentGrowth.class.getDeclaredField("edgeSize"); + field.setAccessible(true); + field.set(component, ds); + } catch (Exception ignored) { + } + + component.schedule(schedule); + + verify(schedule).scheduleRepeating(component, 10, 60); + } + + @Test + public void schedule_calledWithMigrationRateGreaterThanEdgeSize_setsScheduleCorrectly() { + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 60); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + Schedule schedule = mock(Schedule.class); + + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + + int ds = 30; + + try { + Field field = PatchComponentGrowth.class.getDeclaredField("edgeSize"); + field.setAccessible(true); + field.set(component, ds); + } catch (Exception ignored) { + } + + component.schedule(schedule); + + verify(schedule).scheduleRepeating(component, 10, 30); + } + + @Test + public void register_calledWithIncompatibleFeature_throwsException() { + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 60); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + PatchComponentSitesPattern pattern = mock(PatchComponentSitesPattern.class); + doReturn(pattern).when(simMock).getComponent("INCOMPATIBLE"); + + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + + assertThrows( + IncompatibleFeatureException.class, + () -> component.register(simMock, "INCOMPATIBLE")); + } + + @Test + public void register_calledWithoutVEGFLattice_throwsException() { + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 60); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + PatchComponentSitesGraph graph = mock(PatchComponentSitesGraph.class); + when(simMock.getLattice("VEGF")).thenReturn(null); + + try { + Field field = PatchComponentSitesGraph.class.getDeclaredField("graphFactory"); + field.setAccessible(true); + field.set(graph, mock(PatchComponentSitesGraphFactory.class)); + } catch (Exception ignored) { + } + + doReturn(new EnumMap(EdgeDirection.class)) + .when(graph.graphFactory) + .getOffsets(); + + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + assertThrows( + MissingSpecificationException.class, + () -> component.register(simMock, "COMPATIBLE")); + } + + @Test + public void register_calledWithSitesGraphObject_doesNotThrowException() { + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 60); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + PatchComponentSitesGraph graph = mock(PatchComponentSitesGraph.class); + doReturn(graph).when(simMock).getComponent("COMPATIBLE"); + when(simMock.getLattice("VEGF")).thenReturn(latticeMock); + + try { + Field field = PatchComponentSitesGraph.class.getDeclaredField("graphFactory"); + field.setAccessible(true); + field.set(graph, mock(PatchComponentSitesGraphFactory.class)); + } catch (Exception ignored) { + } + + doReturn(new EnumMap(EdgeDirection.class)) + .when(graph.graphFactory) + .getOffsets(); + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + assertDoesNotThrow(() -> component.register(simMock, "COMPATIBLE")); + } + + @Test + public void getDirectionalAverages_givenMap_returnsCorrectMap() { + EnumMap> map = new EnumMap<>(EdgeDirection.class); + map.put(EdgeDirection.UP, new ArrayList<>(Arrays.asList(1.0, 2.0, 3.0))); + map.put(EdgeDirection.DOWN, new ArrayList<>(Arrays.asList(4.0, 5.0, 6.0))); + map.put(EdgeDirection.LEFT, new ArrayList<>(Arrays.asList(7.0, 8.0, 9.0))); + map.put(EdgeDirection.RIGHT, new ArrayList<>(Arrays.asList(10.0, 11.0, 12.0))); + EnumMap averageMap = + PatchComponentGrowth.getDirectionalAverages(map); + assertEquals(averageMap.get(EdgeDirection.UP), 2.0); + assertEquals(averageMap.get(EdgeDirection.DOWN), 5.0); + assertEquals(averageMap.get(EdgeDirection.LEFT), 8.0); + assertEquals(averageMap.get(EdgeDirection.RIGHT), 11.0); + } + + @Test + public void getMaxKey_givenMap_returnsCorrectDirection() { + EnumMap map = new EnumMap<>(EdgeDirection.class); + map.put(EdgeDirection.UP, 1.0); + map.put(EdgeDirection.DOWN, 2.0); + map.put(EdgeDirection.LEFT, 3.0); + map.put(EdgeDirection.RIGHT, 4.0); + assertEquals(EdgeDirection.RIGHT, PatchComponentGrowth.getMaxKey(map)); + } + + @Test + public void normalizeDirectionalMap_givenMap_returnsNormalizedMap() { + EnumMap map = new EnumMap<>(EdgeDirection.class); + map.put(EdgeDirection.UP, 0.5); + map.put(EdgeDirection.DOWN, 0.5); + map.put(EdgeDirection.LEFT, 0.5); + map.put(EdgeDirection.RIGHT, 0.5); + EnumMap normalizedMap = + PatchComponentGrowth.normalizeDirectionalMap(map); + assertEquals(normalizedMap.get(EdgeDirection.UP), 0.25); + assertEquals(normalizedMap.get(EdgeDirection.RIGHT), 0.5); + assertEquals(normalizedMap.get(EdgeDirection.DOWN), 0.75); + assertEquals(normalizedMap.get(EdgeDirection.LEFT), 1.0); + assertEquals(normalizedMap.keySet().size(), 4); + assertNull(normalizedMap.get(EdgeDirection.UNDEFINED)); + } + + @Test + public void sumMap_givenMap_returnsCorrectSum() { + EnumMap map = new EnumMap<>(EdgeDirection.class); + map.put(EdgeDirection.UP, 1.0); + map.put(EdgeDirection.DOWN, 2.0); + map.put(EdgeDirection.LEFT, 3.0); + map.put(EdgeDirection.RIGHT, 4.0); + assertEquals(10.0, PatchComponentGrowth.sumMap(map)); + } + + @Test + public void averageDirectionalMap_givenMap_returnsCorrectAverage() { + EnumMap> map = new EnumMap<>(EdgeDirection.class); + map.put(EdgeDirection.UP, new ArrayList<>(Arrays.asList(1.0, 2.0, 3.0))); + map.put(EdgeDirection.DOWN, new ArrayList<>(Arrays.asList(4.0, 5.0, 6.0))); + map.put(EdgeDirection.LEFT, new ArrayList<>(Arrays.asList(7.0, 8.0, 9.0))); + map.put(EdgeDirection.RIGHT, new ArrayList<>(Arrays.asList(10.0, 11.0, 12.0))); + assertEquals(6.5, PatchComponentGrowth.averageDirectionalMap(map)); + } + + @Test + public void findKeyNodeInMap_givenTargetNode_returnsCorrectKeyNode() { + SiteNode targetNode = new SiteNode(0, 0, 0); + SiteNode originalNode = new SiteNode(1, 1, 0); + SiteNode connectingNode = new SiteNode(1, 0, 0); + + SiteNode node2 = new SiteNode(0, 2, 0); + SiteNode node3 = new SiteNode(0, 3, 0); + SiteNode node4 = new SiteNode(0, 4, 0); + SiteNode node5 = new SiteNode(0, 5, 0); + + // create maps with different key site nodes that contain list of edges where only one + // contains the target node + HashMap> angiogenicNodeMap = new HashMap<>(); + ArrayList edges = new ArrayList<>(); + SiteEdge edge = + new SiteEdge(targetNode, connectingNode, EdgeType.ANGIOGENIC, EdgeLevel.LEVEL_2); + edges.add(edge); + angiogenicNodeMap.put(targetNode, edges); + + edges = new ArrayList<>(); + edge = new SiteEdge(originalNode, connectingNode, EdgeType.ANGIOGENIC, EdgeLevel.LEVEL_2); + edges.add(edge); + angiogenicNodeMap.put(originalNode, edges); + + edges = new ArrayList<>(); + edge = new SiteEdge(node2, node3, EdgeType.ANGIOGENIC, EdgeLevel.LEVEL_2); + edges.add(edge); + angiogenicNodeMap.put(node2, edges); + + edges = new ArrayList<>(); + edge = new SiteEdge(node4, node5, EdgeType.ANGIOGENIC, EdgeLevel.LEVEL_2); + edges.add(edge); + angiogenicNodeMap.put(node4, edges); + + MiniBox parameters = new MiniBox(); + parameters.put("MIGRATION_RATE", 60); + parameters.put("VEGF_THRESHOLD", 0.5); + parameters.put("WALK_TYPE", "BIASED"); + parameters.put("MAX_LENGTH", 100); + parameters.put("CALCULATION_STRATEGY", "DIVERT"); + + PatchComponentSitesGraph graph = mock(PatchComponentSitesGraph.class); + doReturn(graph).when(simMock).getComponent("COMPATIBLE"); + when(simMock.getLattice("VEGF")).thenReturn(latticeMock); + + try { + Field field = PatchComponentSitesGraph.class.getDeclaredField("graphFactory"); + field.setAccessible(true); + field.set(graph, mock(PatchComponentSitesGraphFactory.class)); + } catch (Exception ignored) { + } + + PatchComponentGrowth component = new PatchComponentGrowth(seriesMock, parameters); + + try { + Field field = PatchComponentGrowth.class.getDeclaredField("angiogenicNodeMap"); + field.setAccessible(true); + field.set(component, angiogenicNodeMap); + } catch (Exception ignored) { + } + + assertEquals(targetNode, component.findKeyNodeInMap(connectingNode, originalNode)); + } +} diff --git a/test/arcade/patch/env/component/PatchComponentSitesGraphFactoryTriTest.java b/test/arcade/patch/env/component/PatchComponentSitesGraphFactoryTriTest.java new file mode 100644 index 000000000..271ea6155 --- /dev/null +++ b/test/arcade/patch/env/component/PatchComponentSitesGraphFactoryTriTest.java @@ -0,0 +1,157 @@ +package arcade.patch.env.component; + +import org.junit.jupiter.api.Test; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteEdge; +import arcade.patch.env.component.PatchComponentSitesGraph.SiteNode; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeDirection; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeLevel; +import arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeType; +import arcade.patch.sim.PatchSeries; +import static org.junit.jupiter.api.Assertions.*; +import static org.mockito.Mockito.*; +import static arcade.core.ARCADETestUtilities.*; + +public class PatchComponentSitesGraphFactoryTriTest { + @Test + public void getDirection_givenFromTo_returnsCorrectDirection() { + PatchSeries seriesMock = mock(PatchSeries.class); + + PatchComponentSitesGraphFactoryTri factory = + new PatchComponentSitesGraphFactoryTri(seriesMock); + + int x = randomIntBetween(10, 100); + int y = randomIntBetween(10, 100); + int z = randomIntBetween(0, 100); + + assertEquals( + EdgeDirection.LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 2, y, z), EdgeLevel.VARIABLE)); + assertEquals( + EdgeDirection.RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 2, y, z), EdgeLevel.VARIABLE)); + assertEquals( + EdgeDirection.DOWN_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 1, y + 1, z), EdgeLevel.VARIABLE)); + assertEquals( + EdgeDirection.DOWN_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 1, y + 1, z), EdgeLevel.VARIABLE)); + assertEquals( + EdgeDirection.UP_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 1, y - 1, z), EdgeLevel.VARIABLE)); + assertEquals( + EdgeDirection.UP_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 1, y - 1, z), EdgeLevel.VARIABLE)); + + assertEquals( + EdgeDirection.LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 4, y, z), EdgeLevel.LEVEL_2)); + assertEquals( + EdgeDirection.RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 4, y, z), EdgeLevel.LEVEL_2)); + assertEquals( + EdgeDirection.DOWN_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 2, y + 2, z), EdgeLevel.LEVEL_2)); + assertEquals( + EdgeDirection.DOWN_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 2, y + 2, z), EdgeLevel.LEVEL_2)); + assertEquals( + EdgeDirection.UP_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 2, y - 2, z), EdgeLevel.LEVEL_2)); + assertEquals( + EdgeDirection.UP_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 2, y - 2, z), EdgeLevel.LEVEL_2)); + + assertEquals( + EdgeDirection.LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 8, y, z), EdgeLevel.LEVEL_1)); + assertEquals( + EdgeDirection.RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 8, y, z), EdgeLevel.LEVEL_1)); + assertEquals( + EdgeDirection.DOWN_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 4, y + 4, z), EdgeLevel.LEVEL_1)); + assertEquals( + EdgeDirection.DOWN_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 4, y + 4, z), EdgeLevel.LEVEL_1)); + assertEquals( + EdgeDirection.UP_LEFT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x - 4, y - 4, z), EdgeLevel.LEVEL_1)); + assertEquals( + EdgeDirection.UP_RIGHT, + factory.getDirection( + new SiteNode(x, y, z), new SiteNode(x + 4, y - 4, z), EdgeLevel.LEVEL_1)); + } + + @Test + public void getOppositeDirection_givenEdge_returnsCorrectDirection() { + PatchSeries seriesMock = mock(PatchSeries.class); + + PatchComponentSitesGraphFactoryTri factory = + new PatchComponentSitesGraphFactoryTri(seriesMock); + + int x = randomIntBetween(0, 100); + int y = randomIntBetween(0, 100); + int z = randomIntBetween(0, 100); + + SiteEdge right = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x + 4, y, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + SiteEdge left = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x - 4, y, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + SiteEdge upLeft = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x - 2, y - 2, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + SiteEdge upRight = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x + 2, y - 2, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + SiteEdge downLeft = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x - 2, y + 2, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + SiteEdge downRight = + new SiteEdge( + new SiteNode(x, y, z), + new SiteNode(x + 2, y + 2, z), + EdgeType.CAPILLARY, + EdgeLevel.LEVEL_2); + + assertEquals(EdgeDirection.LEFT, factory.getOppositeDirection(right, right.level)); + assertEquals(EdgeDirection.RIGHT, factory.getOppositeDirection(left, right.level)); + assertEquals(EdgeDirection.DOWN_RIGHT, factory.getOppositeDirection(upLeft, right.level)); + assertEquals(EdgeDirection.DOWN_LEFT, factory.getOppositeDirection(upRight, right.level)); + assertEquals(EdgeDirection.UP_RIGHT, factory.getOppositeDirection(downLeft, right.level)); + assertEquals(EdgeDirection.UP_LEFT, factory.getOppositeDirection(downRight, right.level)); + } +} diff --git a/test/arcade/patch/env/component/PatchComponentSitesGraphTest.java b/test/arcade/patch/env/component/PatchComponentSitesGraphTest.java new file mode 100644 index 000000000..2592890fe --- /dev/null +++ b/test/arcade/patch/env/component/PatchComponentSitesGraphTest.java @@ -0,0 +1,3 @@ +package arcade.patch.env.component; + +public class PatchComponentSitesGraphTest {} diff --git a/test/arcade/patch/env/component/PatchComponentSitesGraphUtiltiesTest.java b/test/arcade/patch/env/component/PatchComponentSitesGraphUtiltiesTest.java new file mode 100644 index 000000000..664ac2bc5 --- /dev/null +++ b/test/arcade/patch/env/component/PatchComponentSitesGraphUtiltiesTest.java @@ -0,0 +1,63 @@ +package arcade.patch.env.component; + +import java.util.ArrayList; +import org.junit.jupiter.api.Test; +import arcade.core.util.Graph; +import static org.junit.jupiter.api.Assertions.*; +import static arcade.patch.env.component.PatchComponentSitesGraph.SiteEdge; +import static arcade.patch.env.component.PatchComponentSitesGraph.SiteNode; +import static arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeLevel; +import static arcade.patch.env.component.PatchComponentSitesGraphFactory.EdgeType; + +public class PatchComponentSitesGraphUtiltiesTest { + @Test + public void getPath_calledWithConnectedNodes_returnsPath() { + Graph graph = new Graph(); + SiteNode node0 = new SiteNode(0, 0, 0); + SiteNode node1 = new SiteNode(0, 1, 0); + SiteNode node2 = new SiteNode(0, 2, 0); + SiteNode node3 = new SiteNode(0, 3, 0); + SiteNode node4 = new SiteNode(0, 4, 0); + SiteNode node5 = new SiteNode(0, 5, 0); + SiteEdge edge0 = new SiteEdge(node0, node1, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge1 = new SiteEdge(node1, node2, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge2 = new SiteEdge(node2, node3, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge3 = new SiteEdge(node3, node4, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge4 = new SiteEdge(node4, node5, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + graph.addEdge(edge0); + graph.addEdge(edge1); + graph.addEdge(edge2); + graph.addEdge(edge3); + graph.addEdge(edge4); + + ArrayList expected = new ArrayList<>(); + expected.add(edge1); + expected.add(edge2); + expected.add(edge3); + + ArrayList actual = PatchComponentSitesGraphUtilities.getPath(graph, node1, node4); + assertIterableEquals(expected, actual); + } + + @Test + public void getPath_calledWithUnconnectedNodes_returnsNull() { + Graph graph = new Graph(); + SiteNode node0 = new SiteNode(0, 0, 0); + SiteNode node1 = new SiteNode(0, 1, 0); + SiteNode node2 = new SiteNode(0, 2, 0); + SiteNode node3 = new SiteNode(0, 3, 0); + SiteNode node4 = new SiteNode(0, 4, 0); + SiteNode node5 = new SiteNode(0, 5, 0); + SiteEdge edge0 = new SiteEdge(node0, node1, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge1 = new SiteEdge(node1, node2, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge2 = new SiteEdge(node2, node3, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + SiteEdge edge4 = new SiteEdge(node4, node5, EdgeType.CAPILLARY, EdgeLevel.LEVEL_1); + graph.addEdge(edge0); + graph.addEdge(edge1); + graph.addEdge(edge2); + graph.addEdge(edge4); + + ArrayList actual = PatchComponentSitesGraphUtilities.getPath(graph, node1, node5); + assertNull(actual); + } +} diff --git a/test/arcade/patch/sim/PatchSeriesTest.java b/test/arcade/patch/sim/PatchSeriesTest.java index 357a0e9d8..8d28989f9 100644 --- a/test/arcade/patch/sim/PatchSeriesTest.java +++ b/test/arcade/patch/sim/PatchSeriesTest.java @@ -80,9 +80,9 @@ public static void setupParameters() { PARAMETERS.addTag(PATCH_PARAMETER_NAMES[i], "PATCH"); PARAMETERS.addAtt(PATCH_PARAMETER_NAMES[i], "value", "" + PATCH_PARAMETER_VALUES[i]); } - MiniBox potts = PARAMETERS.getIdValForTag("PATCH"); - for (String key : potts.getKeys()) { - PATCH.put(key, potts.get(key)); + MiniBox patch = PARAMETERS.getIdValForTag("PATCH"); + for (String key : patch.getKeys()) { + PATCH.put(key, patch.get(key)); } // POPULATION