diff --git a/Examples/PopulationTumor/PopTumor.java b/Examples/PopulationTumor/PopTumor.java index 2f1452e..035ef96 100644 --- a/Examples/PopulationTumor/PopTumor.java +++ b/Examples/PopulationTumor/PopTumor.java @@ -199,7 +199,7 @@ public void Step(int idraw) { } public static void main(String[] args) { - PopTumor ft = new PopTumor(150, 0.02, 0.001, 0.003, 100000, 2e-6, 2e-6, 0, 2, 100); + PopTumor ft = new PopTumor(150, 0.015, 0.001, 0.004, 100000, 2e-6, 2e-6, 0, 2, 100); ft.Init(100); int idraw = -1; for (int i = 0; i < 4001; i++) { diff --git a/Examples/_3OffLatticeExample/ExampleOffLattice.java b/Examples/_3OffLatticeExample/ExampleOffLattice.java index bf5885d..a7979ae 100644 --- a/Examples/_3OffLatticeExample/ExampleOffLattice.java +++ b/Examples/_3OffLatticeExample/ExampleOffLattice.java @@ -22,7 +22,7 @@ public void Init(int color){ if(overlap<0) { return 0;//if cells aren't actually overlapping, then there is no force response } - return Math.pow(G.FORCE_SCALER*overlap, G.FORCE_EXPONENT); + return G.FORCE_SCALER*overlap;//this constant scaling of the overlap is called Hooke's law! } public void CalcMove(){ //sets x and y velocity components of cell @@ -47,15 +47,14 @@ public void MoveDiv(){ public class ExampleOffLattice extends AgentGrid2D { static final int WHITE=RGB256(248,255,252), PURPLE =RGB256(77,0,170), PINK =RGB256(222,0,109),CYTOPLASM=RGB256(191,156,147); - double RADIUS=0.5; + double RADIUS=0.25; - double FORCE_EXPONENT=2;//these constants have been found to be rather stable, but tweak them and see what happens! - double FORCE_SCALER=0.7; + double FORCE_SCALER=0.25;//this constant was found to be rather stable, but tweak it and see what happens! double FRICTION=0.5; double PURP_DIV_BIAS =0.01;//grid holds onto phenotype differences, if drift were included these would probably be moved to individual cells. double PINK_DIV_BIAS =0.02; - double PURP_INHIB_WEIGHT =0.015; + double PURP_INHIB_WEIGHT =0.02; double PINK_INHIB_WEIGHT =0.05; ArrayList neighborList=new ArrayList<>(); ArrayList neighborInfo=new ArrayList<>(); @@ -72,7 +71,7 @@ public ExampleOffLattice(int x, int y,String outFileName) { out=new FileIO(outFileName,"w"); } public static void main(String[] args) { - int x=60,y=60; + int x=30,y=30; //to record output, call the constructor with an output filename ExampleOffLattice ex=new ExampleOffLattice(x,y,"PopOut.csv"); //ExampleOffLattice ex=new ExampleOffLattice(x,y); diff --git a/HAL/GridsAndAgents/Agent2DBase.java b/HAL/GridsAndAgents/Agent2DBase.java index bfe9f6c..c9286f1 100644 --- a/HAL/GridsAndAgents/Agent2DBase.java +++ b/HAL/GridsAndAgents/Agent2DBase.java @@ -138,5 +138,6 @@ public double DistSquared(double x,double y){ double dy = DispY(y); return NormSquared(dx,dy); } + } diff --git a/HAL/GridsAndAgents/AgentBaseSpatial.java b/HAL/GridsAndAgents/AgentBaseSpatial.java index db44ef7..e36c043 100644 --- a/HAL/GridsAndAgents/AgentBaseSpatial.java +++ b/HAL/GridsAndAgents/AgentBaseSpatial.java @@ -18,25 +18,6 @@ public int Isq() { return iSq; } - /** - * swaps the positions of two agents. useful mostly for the AgentSQ2unstackable and AgentSQ3unstackable classes, - * which don't allow stacking of agents, making this maneuver otherwise difficult. - */ - public void SwapPosition(AgentBaseSpatial other) { - if (!alive || !other.alive) { - throw new RuntimeException("attempting to move dead agent"); - } - if (other.G != G) { - throw new IllegalStateException("can't swap positions between agents on different grids!"); - } - int iOther = other.Isq(); - int iThis = Isq(); - other.RemSQ(); - this.RemSQ(); - other.MoveSQ(iThis); - this.MoveSQ(iOther); - } - /** * moves the agent to the middle of the square at the indices/index specified */ diff --git a/HAL/GridsAndAgents/AgentSQ1Dunstackable.java b/HAL/GridsAndAgents/AgentSQ1Dunstackable.java index 7a51a54..1a11207 100644 --- a/HAL/GridsAndAgents/AgentSQ1Dunstackable.java +++ b/HAL/GridsAndAgents/AgentSQ1Dunstackable.java @@ -120,6 +120,26 @@ void Setup(int i) { iSq=i; AddSQ(i); } + /** + * swaps the positions of two agents. useful for the AgentSQunstackable classes, + * which don't allow stacking of agents, making this maneuver otherwise impossible. + */ + public void SwapPosition(AgentSQ1Dunstackable other) { + if (!alive || !other.alive) { + throw new RuntimeException("attempting to move dead agent"); + } + if (other.G != G) { + throw new IllegalStateException("can't swap positions between agents on different grids!"); + } + other.RemSQ(); + this.RemSQ(); + int iThis = this.iSq; + this.iSq = other.iSq; + other.iSq=iThis; + other.AddSQ(other.iSq); + this.AddSQ(this.iSq); + } + @Override void Setup(int x, int y) { diff --git a/HAL/GridsAndAgents/AgentSQ2Dunstackable.java b/HAL/GridsAndAgents/AgentSQ2Dunstackable.java index 2428104..4205d80 100644 --- a/HAL/GridsAndAgents/AgentSQ2Dunstackable.java +++ b/HAL/GridsAndAgents/AgentSQ2Dunstackable.java @@ -164,6 +164,32 @@ public int Isq() { return iSq; } + /** + * swaps the positions of two agents. useful for the AgentSQunstackable classes, + * which don't allow stacking of agents, making this maneuver otherwise impossible. + */ + public void SwapPosition(AgentSQ2Dunstackable other) { + if (!alive || !other.alive) { + throw new RuntimeException("attempting to move dead agent"); + } + if (other.G != G) { + throw new IllegalStateException("can't swap positions between agents on different grids!"); + } + other.RemSQ(); + this.RemSQ(); + int xThis=this.xSq; + int yThis=this.ySq; + int iThis = this.iSq; + this.xSq = other.xSq; + this.ySq = other.ySq; + this.iSq = other.iSq; + other.xSq=xThis; + other.ySq=yThis; + other.iSq=iThis; + other.AddSQ(other.iSq); + this.AddSQ(this.iSq); + } + void Setup(double i) { Setup((int) i); } diff --git a/HAL/GridsAndAgents/AgentSQ3Dunstackable.java b/HAL/GridsAndAgents/AgentSQ3Dunstackable.java index 178b678..76782a1 100644 --- a/HAL/GridsAndAgents/AgentSQ3Dunstackable.java +++ b/HAL/GridsAndAgents/AgentSQ3Dunstackable.java @@ -181,6 +181,36 @@ void Setup(double x, double y) { } + /** + * swaps the positions of two agents. useful for the AgentSQunstackable classes, + * which don't allow stacking of agents, making this maneuver otherwise impossible. + */ + public void SwapPosition(AgentSQ3Dunstackable other) { + if (!alive || !other.alive) { + throw new RuntimeException("attempting to move dead agent"); + } + if (other.G != G) { + throw new IllegalStateException("can't swap positions between agents on different grids!"); + } + other.RemSQ(); + this.RemSQ(); + int xThis=this.xSq; + int yThis=this.ySq; + int zThis=this.zSq; + int iThis = this.iSq; + this.xSq = other.xSq; + this.ySq = other.ySq; + this.zSq=other.zSq; + this.iSq = other.iSq; + other.xSq=xThis; + other.ySq=yThis; + other.zSq=zThis; + other.iSq=iThis; + other.AddSQ(other.iSq); + this.AddSQ(this.iSq); + } + + @Override void Setup(int xSq, int ySq, int zSq) { this.xSq = xSq; diff --git a/HAL/GridsAndAgents/PDEGrid1D.java b/HAL/GridsAndAgents/PDEGrid1D.java index d1509ea..dc19e79 100644 --- a/HAL/GridsAndAgents/PDEGrid1D.java +++ b/HAL/GridsAndAgents/PDEGrid1D.java @@ -52,6 +52,9 @@ public PDEGrid1D(int xDim, boolean wrapX) { public double Get(int x) { return field[x]; } + public double Get(double x) { + return Get((int)x); + } /** * sets the prev field value at the specified coordinates @@ -59,6 +62,10 @@ public double Get(int x) { public void Set(int x, double val) { deltas[x] = val-field[x]; } + public void Set(double x, double val) { + Set((int)x,val); + } + /** * adds to the prev field value at the specified index @@ -66,6 +73,9 @@ public void Set(int x, double val) { public void Add(int x, double val) { deltas[x] += val; } + public void Add(double x, double val) { + Add((int)x,val); + } /** * multiplies a value in the “current field” and adds the change to the “delta field” @@ -73,6 +83,9 @@ public void Add(int x, double val) { public void Mul(int x, double val) { deltas[x] += field[x] * val; } + public void Mul(double x, double val) { + Mul((int)x,val); + } /** * adds the delta field into the current field, also increments the tick. diff --git a/HAL/GridsAndAgents/PDEGrid2D.java b/HAL/GridsAndAgents/PDEGrid2D.java index f1595db..0d12b0e 100644 --- a/HAL/GridsAndAgents/PDEGrid2D.java +++ b/HAL/GridsAndAgents/PDEGrid2D.java @@ -176,7 +176,9 @@ public static void Diffusion2DADIBetweenAction(double[]field, double[]scratch,do public double Get(int x, int y) { return field[x * yDim + y]; } - + public double Get(double x, double y) { + return Get((int)x,(int)y); + } /** * gets the prev field value at the specified index */ @@ -190,6 +192,9 @@ public double Get(int i) { public void Set(int x, int y, double val) { deltas[x * yDim + y] = val - field[x * yDim + y]; } + public void Set(double x, double y, double val) { + Set((int)x,(int)y,val); + } /** * sets the prev field value at the specified index @@ -204,6 +209,9 @@ public void Set(int i, double val) { public void Add(int x, int y, double val) { deltas[x * yDim + y] += val; } + public void Add(double x, double y, double val) { + Add((int)x,(int)y,val); + } /** * adds to the prev field value at the specified index @@ -218,6 +226,9 @@ public void Add(int i, double val) { public void Mul(int x, int y, double val) { deltas[x * yDim + y] += field[x * yDim + y] * val; } + public void Mul(double x, double y, double val) { + Mul((int)x,(int)y,val); + } /** * multiplies a value in the “current field” and adds the change to the “delta field” @@ -488,6 +499,9 @@ public double GradientX(int x, int y) { double right = PDEequations.DisplacedX2D(field,x + 1, y, xDim, yDim,wrapX,(X,Y)->Get(X-1,Y)); return right - left; } + public double GradientX(double x, double y) { + return GradientX((int)x,(int)y); + } /** * returns the gradient of the diffusible in the Y direction at the coordinates specified @@ -497,6 +511,9 @@ public double GradientY(int x, int y) { double up = PDEequations.DisplacedY2D(field,x, y+1, xDim, yDim,wrapX,(X,Y)->Get(X,Y-1)); return up - down; } + public double GradientY(double x, double y) { + return GradientY((int)x,(int)y); + } /** * returns the gradient of the diffusible in the X direction at the coordinates specified, will use the boundary @@ -507,6 +524,9 @@ public double GradientX(int x, int y, double boundaryCond) { double right = PDEequations.DisplacedX2D(field,x + 1, y, xDim, yDim,wrapX,(X,Y)->boundaryCond); return right - left; } + public double GradientX(double x, double y, double boundaryCond) { + return GradientX((int)x,(int)y,boundaryCond); + } /** * returns the gradient of the diffusible in the Y direction at the coordinates specified, will use the boundary @@ -517,6 +537,9 @@ public double GradientY(int x, int y, double boundaryCond) { double up = PDEequations.DisplacedY2D(field,x, y+1, xDim, yDim,wrapX,(X,Y)->boundaryCond); return up - down; } + public double GradientY(double x, double y, double boundaryCond) { + return GradientY((int)x,(int)y,boundaryCond); + } /** * ensures that all values will be non-negative on the next timestep, call before Update diff --git a/HAL/GridsAndAgents/PDEGrid3D.java b/HAL/GridsAndAgents/PDEGrid3D.java index a07ff1a..4d3633e 100644 --- a/HAL/GridsAndAgents/PDEGrid3D.java +++ b/HAL/GridsAndAgents/PDEGrid3D.java @@ -64,6 +64,9 @@ public PDEGrid3D(int xDim, int yDim, int zDim) { public double Get(int x, int y, int z) { return field[x * yDim * zDim + y * zDim + z]; } + public double Get(double x, double y, double z) { + return Get((int)x,(int)y,(int)z); + } /** * gets the delta field value at the specified index @@ -78,6 +81,9 @@ public double Get(int i) { public void Set(int x, int y, int z, double val) { deltas[x * yDim * zDim + y * zDim + z] = val - field[x * yDim * zDim + y * zDim + z]; } + public void Set(double x, double y, double z, double val) { + Set((int)x,(int)y,(int)z,val); + } /** * sets the current field value at the specified index @@ -93,6 +99,9 @@ public void Set(int i, double val) { public void Add(int x, int y, int z, double val) { deltas[x * yDim * zDim + y * zDim + z] += val; } + public void Add(double x, double y, double z, double val) { + Add((int)x,(int)y,(int)z,val); + } /** * adds to the delta field value at the specified index @@ -112,7 +121,11 @@ public void Mul(int i, double val) { * multiplies a value in the “current field” and adds the change to the “delta field” */ public void Mul(int x, int y, int z, double val) { - deltas[x * yDim * zDim + y * zDim + z] += field[x * yDim * zDim + y * zDim + z] * val; + int i = x * yDim * zDim + y * zDim + z; + deltas[i] += field[i] * val; + } + public void Mul(double x, double y, double z, double val) { + Mul((int)x,(int)y,(int)z,val); } /** @@ -409,6 +422,9 @@ public double GradientX(int x, int y,int z) { double right = PDEequations.DisplacedX3D(field,x + 1, y,z, xDim, yDim,zDim,wrapX,(X,Y,Z)->Get(X-1,Y,Z)); return right - left; } + public double GradientX(double x, double y,double z) { + return GradientX((int)x,(int)y,(int)z); + } /** * returns the gradient of the diffusible in the Y direction at the coordinates specified @@ -418,6 +434,9 @@ public double GradientY(int x, int y,int z) { double up = PDEequations.DisplacedY3D(field,x, y+1,z, xDim, yDim,zDim,wrapY,(X,Y,Z)->Get(X,Y-1,Z)); return up - down; } + public double GradientY(double x, double y,double z) { + return GradientY((int)x,(int)y,(int)z); + } /** * returns the gradient of the diffusible in the Z direction at the coordinates specified */ @@ -426,6 +445,9 @@ public double GradientZ(int x, int y,int z) { double up = PDEequations.DisplacedY3D(field,x, y,z+1, xDim, yDim,zDim,wrapZ,(X,Y,Z)->Get(X,Y,Z-1)); return up - down; } + public double GradientZ(double x, double y,double z) { + return GradientZ((int)x,(int)y,(int)z); + } /** * returns the gradient of the diffusible in the X direction at the coordinates specified, will use the boundary @@ -436,6 +458,9 @@ public double GradientX(int x, int y,int z,double boundaryCond) { double right = PDEequations.DisplacedX3D(field,x + 1, y,z, xDim, yDim,zDim,wrapX,(X,Y,Z)->boundaryCond); return right - left; } + public double GradientX(double x, double y,double z,double boundaryCond) { + return GradientX((int)x,(int)y,(int)z,boundaryCond); + } /** * returns the gradient of the diffusible in the Y direction at the coordinates specified, will use the boundary @@ -446,6 +471,9 @@ public double GradientY(int x, int y,int z,double boundaryCond) { double up = PDEequations.DisplacedY3D(field,x, y+1,z, xDim, yDim,zDim,wrapY,(X,Y,Z)->boundaryCond); return up - down; } + public double GradientY(double x, double y,double z,double boundaryCond) { + return GradientY((int)x,(int)y,(int)z,boundaryCond); + } /** * returns the gradient of the diffusible in the Z direction at the coordinates specified, will use the boundary @@ -456,6 +484,9 @@ public double GradientZ(int x, int y,int z,double boundaryCond) { double up = PDEequations.DisplacedY3D(field,x, y,z+1, xDim, yDim,zDim,wrapZ,(X,Y,Z)->boundaryCond); return up - down; } + public double GradientZ(double x, double y,double z,double boundaryCond) { + return GradientZ((int)x,(int)y,(int)z,boundaryCond); + } @Override diff --git a/HAL/Rand.java b/HAL/Rand.java index 5b60083..dbd0097 100644 --- a/HAL/Rand.java +++ b/HAL/Rand.java @@ -55,6 +55,9 @@ public int Int(int bound) { public double Double(double bound) { return rn.Double(bound); } + public double Double(int bound) { + return rn.Double(bound); + } /** * returns a random double from 0 to 1 diff --git a/HAL/Tools/FileIO.java b/HAL/Tools/FileIO.java index 82b360e..99df1d5 100644 --- a/HAL/Tools/FileIO.java +++ b/HAL/Tools/FileIO.java @@ -11,58 +11,108 @@ public class FileIO { public final String fileName; public final char ReadWriteAppend; - public final boolean binary; + public final char mode; public final BufferedReader reader; public final BufferedWriter writer; + //binary reader and writer public final DataOutputStream writerBin; public final DataInputStream readerBin; + //object reader and writer + public final FileInputStream serialfileReader; + public final ObjectInputStream serialobjectReader; + public final FileOutputStream serialfileWriter; + public final ObjectOutputStream serialobjectWriter; boolean isClosed = false; /** * @param fileName name of the file to read from or write to - * @param mode should be either "r":read, "w":write, "rb":readBinary, "wb":writeBinary + * @param mode should be either "r":read, "w":write, "rb":readBinary, "wb":writeBinary, "rs": readSerialization, "ws": writeSerialization */ public FileIO(String fileName, String mode) { char[] modeChars = mode.toCharArray(); - if (modeChars.length > 2 || (modeChars.length > 1 && modeChars[1] != 'b') || (modeChars[0] != 'w' && modeChars[0] != 'r' && modeChars[0] != 'a')) { - throw new IllegalArgumentException("inccorect mode argument! mode should be 'r' for read, 'w' for write, or 'a' for append, followed by optional 'b' for binary"); - } - this.fileName = fileName; - this.ReadWriteAppend = modeChars[0]; - this.binary = modeChars.length > 1; - boolean appendOut = false; - if (ReadWriteAppend == 'a') { - appendOut = true; - } - - BufferedReader reader = null; - BufferedWriter writer = null; - DataOutputStream writerBin = null; - DataInputStream readerBin = null; - try { - if (ReadWriteAppend == 'r') { - if (this.binary) { - readerBin = new DataInputStream(new BufferedInputStream(new FileInputStream(fileName))); - - } else { - reader = new BufferedReader(new FileReader(fileName)); + if(modeChars.length==2&&modeChars[1]=='s' && (modeChars[0]=='r'||modeChars[0]=='w')) { + reader=null; + writer=null; + writerBin=null; + readerBin=null; + this.fileName=fileName; + this.ReadWriteAppend=modeChars[0]; + this.mode='s'; + ObjectOutputStream serialobjectWriter=null; + FileOutputStream serialfileWriter=null; + ObjectInputStream serialobjectReader=null; + FileInputStream serialfileReader=null; + if (modeChars[0] == 'r') { + //setup serialization reader + serialobjectWriter=null; + serialfileWriter=null; + try { + serialfileReader=new FileInputStream(fileName); + serialobjectReader=new ObjectInputStream(serialfileReader); + } catch (IOException e) { + e.printStackTrace(); } - } else if (ReadWriteAppend == 'w' || ReadWriteAppend == 'a') { - if (this.binary) { - writerBin = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(fileName, appendOut))); + }else { + //setup serialization writer + serialobjectReader=null; + serialfileReader=null; + try { + serialfileWriter=new FileOutputStream(fileName); + serialobjectWriter=new ObjectOutputStream(serialfileWriter); + } catch (IOException e) { + e.printStackTrace(); + } + } + this.serialfileReader=serialfileReader; + this.serialobjectReader=serialobjectReader; + this.serialfileWriter=serialfileWriter; + this.serialobjectWriter=serialobjectWriter; + } + else { + if (modeChars.length > 2 || (modeChars.length > 1 && modeChars[1] != 'b') || (modeChars[0] != 'w' && modeChars[0] != 'r' && modeChars[0] != 'a')) { + throw new IllegalArgumentException("inccorect mode argument! mode should be 'r' for read, 'w' for write, or 'a' for append, followed by optional 'b' for binary or 's' for serialization"); + } + this.fileName = fileName; + this.ReadWriteAppend = modeChars[0]; + this.mode = modeChars.length > 1?'s':'n'; + boolean appendOut = false; + if (ReadWriteAppend == 'a') { + appendOut = true; + } + + BufferedReader reader = null; + BufferedWriter writer = null; + DataOutputStream writerBin = null; + DataInputStream readerBin = null; + try { + if (ReadWriteAppend == 'r') { + if (this.mode=='b') { + readerBin = new DataInputStream(new BufferedInputStream(new FileInputStream(fileName))); + + } else { + reader = new BufferedReader(new FileReader(fileName)); + } + } else if (ReadWriteAppend == 'w' || ReadWriteAppend == 'a') { + if (this.mode=='b') { + writerBin = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(fileName, appendOut))); + } else { + writer = new BufferedWriter(new FileWriter(fileName, appendOut)); + } } else { - writer = new BufferedWriter(new FileWriter(fileName, appendOut)); + throw new IllegalArgumentException("rwa character must be one of r(read) w(write) or a(append)"); } - } else { - throw new IllegalArgumentException("rwa character must be one of r(read) w(write) or a(append)"); + } catch (IOException e) { + e.printStackTrace(); } - } catch (IOException e) { - e.printStackTrace(); + this.reader = reader; + this.writer = writer; + this.writerBin = writerBin; + this.readerBin = readerBin; + this.serialfileReader=null; + this.serialfileWriter=null; + this.serialobjectReader=null; + this.serialobjectWriter=null; } - this.reader = reader; - this.writer = writer; - this.writerBin = writerBin; - this.readerBin = readerBin; } /** @@ -86,7 +136,7 @@ public String[] ReadLineDelimit(String delimiter) { String[] read = null; try { String line = reader.readLine(); - if (line != null) { + if (line != null && !line.equals("")) { read = line.split(delimiter); } } catch (IOException e) { @@ -253,6 +303,8 @@ public ArrayList Read() { } return ret; } + + //WRITE FUNCTIONS /** @@ -572,6 +624,24 @@ public void ReadBinBools(boolean[] putHere) { putHere[i] = ReadBinBool(); } } + //SERIALIZATON FUNCTIONS + public Object ReadObject(){ + try { + return serialobjectReader.readObject(); + } catch (IOException e) { + e.printStackTrace(); + } catch (ClassNotFoundException e) { + e.printStackTrace(); + } + return null; + } + public void WriteObject(Object toSave){ + try { + serialobjectWriter.writeObject(toSave); + } catch (IOException e) { + e.printStackTrace(); + } + } /** * returns the length of the file, or 0 if the file does not exist @@ -590,14 +660,20 @@ public void Close() { try { this.isClosed = true; if (ReadWriteAppend == 'r') { - if (binary) { + if (this.mode == 'b') { readerBin.close(); + } else if (this.mode == 's') { + serialfileReader.close(); + serialobjectReader.close(); } else { reader.close(); } } else if (ReadWriteAppend == 'w' || ReadWriteAppend == 'a') { - if (binary) { + if (this.mode == 'b') { writerBin.close(); + } else if (this.mode == 's') { + serialfileWriter.close(); + serialobjectWriter.close(); } else { writer.close(); } diff --git a/HAL/Util.java b/HAL/Util.java index 3374f97..c7706ac 100644 --- a/HAL/Util.java +++ b/HAL/Util.java @@ -772,6 +772,13 @@ public static String ArrToString(T[] arr, String delim, int start, int end) return out; } + public static double InterpolateLinear(double x, double startX, double endX, double startY, double endY){ + if(x>endX||x { + //model constants + public final static int RESISTANT = RGB(0, 1, 0), SENSITIVE = RGB(0, 0, 1); + public double TIMESTEP=2.0/24;//2 hours per timestep + public double SPACE_STEP=20;//um + public double DIV_PROB_SEN = ProbScale(0.5,TIMESTEP); + public double DIV_PROB_RES = ProbScale(0.2,TIMESTEP); + public double DEATH_PROB = ProbScale(0.02,TIMESTEP); + public double DRUG_START = 20/TIMESTEP; + public double DRUG_PERIOD = 15/TIMESTEP; + public double DRUG_DURATION = 3/TIMESTEP; + public double DRUG_DIFF_RATE = 0.02*60*60*24*(TIMESTEP/(SPACE_STEP*SPACE_STEP));//diffusion rate in um/seconds + public double DRUG_UPTAKE = -0.03 *TIMESTEP; + public double DRUG_DEATH = ProbScale(0.8,TIMESTEP); + public double DRUG_BOUNDARY_VAL = 1.0; + public long popTotal; + + double FORCE_SCALER=0.5; + + //public double DRUG_UPTAKE = 0; + //internal model objects + public PDEGrid2D drug; + public Rand rn; + public double[]scratch=new double[2]; + + public ExampleModel(int xDim, int yDim, Rand rn) { + super(xDim, yDim, ExampleCell.class); + this.rn = rn; + drug = new PDEGrid2D(xDim, yDim); + } + + @FunctionalInterface + interface ModelStep{ + void Run(ExampleModel m, int tick); + } + + public static void RunModel(int sideLen,ModelStep Step,boolean draw){ + int x = sideLen, y = sideLen; + ExampleModel m=new ExampleModel(x,y,new Rand(0)); + OpenGL2DWindow win =null; + if(draw) { + win = new OpenGL2DWindow(1000, 1000, m.xDim, m.yDim); + } + m.DRUG_START=0; + m.DRUG_DURATION=m.DRUG_PERIOD; + m.InitTumor(); + for (int tick = 0; tick < 10000; tick++) { + Step.Run(m,tick); + if(draw) { + win.Clear(BLACK); + for (ExampleCell cell : m) { + win.Circle(cell.Xpt(), cell.Ypt(), cell.radius/3, HeatMapGRB(m.drug.Get(cell.Xpt(),cell.Ypt())*0.8+0.2)); + } + win.Update(); + } + } + System.out.println("sideLen:"+sideLen+" AvgPop:"+m.popTotal*1.0/10000); + if(draw) { + win.Close(); + } + } + + public static void main(String[] args) { +// AwaitInput(); + long start=System.currentTimeMillis(); + RunModel(60, ExampleModel::ModelStep60,false); + RunModel(90, ExampleModel::ModelStep90,false); + RunModel(120, ExampleModel::ModelStep120,false); + RunModel(150, ExampleModel::ModelStep150,false); + RunModel(180, ExampleModel::ModelStep180,false); + System.out.println(System.currentTimeMillis()-start); + } + + public void InitTumor() { + for (int i = 0; i < drug.length; i++) { + ExampleCell c=NewAgentPT(rn.Double()*xDim,rn.Double()*yDim); + c.type=RESISTANT; + c.radius=0.5; + } + } + + public void DiffusionStep(int tick){ + if (tick > DRUG_START && (tick - DRUG_START) % DRUG_PERIOD < DRUG_DURATION) { + drug.DiffusionADI(DRUG_DIFF_RATE, DRUG_BOUNDARY_VAL); + } else { + drug.DiffusionADI(DRUG_DIFF_RATE); + } + drug.Update(); + } + public void StepAllCells(int tick){ + for (ExampleCell cell : this) { + cell.BirthDeath(); + } + for (ExampleCell cell : this) { + cell.Movement(); + } + } + + + public static void ModelStep60(ExampleModel m, int tick) { + m.StepAllCells(tick); + m.DiffusionStep(tick); + m.popTotal+=m.Pop(); + } + public static void ModelStep90(ExampleModel m, int tick) { + m.StepAllCells(tick); + m.DiffusionStep(tick); + m.popTotal+=m.Pop(); + } + public static void ModelStep120(ExampleModel m, int tick) { + m.StepAllCells(tick); + m.DiffusionStep(tick); + m.popTotal+=m.Pop(); + } + public static void ModelStep150(ExampleModel m, int tick) { + m.StepAllCells(tick); + m.DiffusionStep(tick); + m.popTotal+=m.Pop(); + } + public static void ModelStep180(ExampleModel m, int tick) { + m.StepAllCells(tick); + m.DiffusionStep(tick); + m.popTotal+=m.Pop(); + } + +} + +class ExampleCell extends SphericalAgent2D { + public int type; + + public void BirthDeath() { + //Consumption of Drug + G.drug.Mul(Xpt(),Ypt(), G.DRUG_UPTAKE); + //Chance of Death, depends on resistance and drug concentration + if (G.rn.Double() < G.DEATH_PROB + (type == RESISTANT ? 0 : G.drug.Get(Xpt(),Ypt()) * G.DRUG_DEATH)) { + Dispose(); + return; + } + double pressure=SumForces(1,(overlap, other) -> overlap*G.FORCE_SCALER); + //contact inhibition and division probability influence division event + if (G.rn.Double()*pressure*1000 < (type == RESISTANT ? G.DIV_PROB_RES : G.DIV_PROB_SEN)) { + ExampleCell c=Divide(radius*1.0/3,G.scratch,G.rn); + c.radius=radius; + c.xVel=0; + c.yVel=0; + } + + } + public void Movement() { + ForceMove(); + ApplyFriction(0); + } +} + diff --git a/Testing/PerformanceTestModels/PerformanceTestOnLattice/ExampleModel.java b/Testing/PerformanceTestModels/PerformanceTestOnLattice/ExampleModel.java index 27857b0..67e5895 100644 --- a/Testing/PerformanceTestModels/PerformanceTestOnLattice/ExampleModel.java +++ b/Testing/PerformanceTestModels/PerformanceTestOnLattice/ExampleModel.java @@ -53,12 +53,15 @@ public static void RunModel(int sideLen,ModelStep Step){ } public static void main(String[] args) { - AwaitInput(); + long start=System.currentTimeMillis(); RunModel(60, ExampleModel::ModelStep60); RunModel(90, ExampleModel::ModelStep90); RunModel(120, ExampleModel::ModelStep120); RunModel(150, ExampleModel::ModelStep150); RunModel(180, ExampleModel::ModelStep180); + System.out.println(System.currentTimeMillis()-start); + + } public void InitTumor() { diff --git a/Testing/PerformanceTestModels/PerformanceTestSpherical/ExampleModel.java b/Testing/PerformanceTestModels/PerformanceTestSpherical/ExampleModel.java index 27e3673..32d0878 100644 --- a/Testing/PerformanceTestModels/PerformanceTestSpherical/ExampleModel.java +++ b/Testing/PerformanceTestModels/PerformanceTestSpherical/ExampleModel.java @@ -26,8 +26,6 @@ public class ExampleModel extends AgentGrid2D { public double DRUG_BOUNDARY_VAL = 1.0; public long popTotal; - double FORCE_EXPONENT=2;//these constants have been found to be rather stable, but tweak them and see what happens! - //double FORCE_SCALER=0.7; double FORCE_SCALER=0.5; //public double DRUG_UPTAKE = 0; @@ -74,12 +72,14 @@ public static void RunModel(int sideLen,ModelStep Step,boolean draw){ } public static void main(String[] args) { - AwaitInput(); + long start=System.currentTimeMillis(); +// AwaitInput(); RunModel(60, ExampleModel::ModelStep60,false); RunModel(90, ExampleModel::ModelStep90,false); RunModel(120, ExampleModel::ModelStep120,false); RunModel(150, ExampleModel::ModelStep150,false); RunModel(180, ExampleModel::ModelStep180,false); + System.out.println(System.currentTimeMillis()-start); } public void InitTumor() { @@ -150,7 +150,7 @@ public void BirthDeath() { double pressure=SumForces(0.5,(overlap, other) -> overlap*G.FORCE_SCALER); //contact inhibition and division probability influence division event if (G.rn.Double()*pressure*1000 < (type == RESISTANT ? G.DIV_PROB_RES : G.DIV_PROB_SEN)) { - ExampleCell c=Divide(radius*2.0/3,G.scratch,G.rn); + ExampleCell c=Divide(radius*1.0/3,G.scratch,G.rn); c.radius=radius; c.xVel=0; c.yVel=0;