Version 1.2.0--Voxel geometry improvements, CPF refactoring

 - Streamline capabilities for creating and visualizing voxel versions of 3D stochastic geometries
 - Enable particle simulations on voxel geometries
 - New features to compute spatial statistics quantities from voxel geometries
 - Refactor conditional probability function (CPF) functionality to enable new CPFs

Author: aaronjeffreyolson

.gitignore

@@ -6,5 +6,6 @@
 **/*pdf
 **/*npy
 **/*csv
+**/*vtk
 **/*~
 

Core/3D/CPF_Basepy.py

@@ -0,0 +1,25 @@
+#!usr/bin/env python
+import sys
+sys.path.append('/../../Classes/Tools')
+from ClassToolspy import ClassTools
+import numpy as np
+
+## \brief Base class for conditional probability functions which return probability of material type at one point based on known material type at others
+# \author Aaron Olson, [email protected], [email protected]
+#
+class CPF_Base(ClassTools):
+    def __init__(self):
+        super(CPF_Base,self).__init__()
+
+    def __str__(self):
+        return str(self.__dict__)
+        
+    ## \brief Defines mixing parameters to use in conditional probability function evaluations
+    # To be provided by non-base class
+    def defineMixingParams(self):
+        assert 1==0
+
+    ## \brief Computes conditional probabilities from governing points
+    # To be provided by non-base class
+    def returnConditionalProbabilities(self):
+        assert 1==0

Core/3D/CPF_MarkovianIndependentContribspy.py

@@ -0,0 +1,63 @@
+#!usr/bin/env python
+from CPF_Basepy import CPF_Base
+from MarkovianInputspy import MarkovianInputs
+import numpy as np
+
+## \brief Conditional probability function evaluator for model where each point is assumed to contribute to the new point's material likelihoods independently
+# \author Aaron Olson, [email protected], [email protected]
+#
+class CPF_MarkovianIndependentContribs(CPF_Base,MarkovianInputs):
+    def __init__(self):
+        super(CPF_MarkovianIndependentContribs,self).__init__()
+
+    def __str__(self):
+        return str(self.__dict__)
+        
+    ## \brief Defines mixing parameters (chord lengths and probabilities)
+    #
+    # \param[in] lam, list of floats, list of chord lengths
+    # \returns sets chord length, probabilities, correlation length, seed locations
+    def defineMixingParams(self,*args):
+        # Assert material chord lengths and slab length and store in object
+        lam = args[0][0]
+        assert isinstance(lam,list)
+        self.nummats = len(lam)
+        for i in range(0,self.nummats): assert isinstance(lam[i],float) and lam[i]>0.0
+        self.lam = lam
+
+        # Solve and store material probabilities and correlation length
+        self.solveNaryMarkovianParamsBasedOnChordLengths(self.lam)
+
+    ## \brief Computes conditional probabilities from governing points
+    #
+    # \returns sets condprob list of floats, probabilities of selecting each material type
+    def returnConditionalProbabilities(self,locpoint,govpoints,govmatinds):
+        govpoints = np.asarray(govpoints)
+        self.govdists =                                       np.power(np.subtract(govpoints[:,0],locpoint[0]),2)
+        self.govdists =                np.add( self.govdists, np.power(np.subtract(govpoints[:,1],locpoint[1]),2) )
+        self.govdists = list( np.sqrt( np.add( self.govdists, np.power(np.subtract(govpoints[:,2],locpoint[2]),2) ) ) )
+
+        compliments = np.ones( self.nummats )
+        #for loop through governing points, collect compliments of new point in same region as each material type
+        for igovpt in range(0,len(govpoints)):
+            imat = govmatinds[igovpt]
+            compliments[imat] *= self.freqAtLeastOnePseudointerface_Markovian(igovpt)
+        #compute fraction of combination space for which new point is not in the same cell as any governing points
+        frac_newregion = np.prod( compliments )
+        #compute fraction of combination space for which new point is in the same cell as governing points of each type
+        frac_sameregion = []
+        for imat in range(0,self.nummats):
+            frac_sameregion.append( frac_newregion/compliments[imat] - frac_newregion )
+        #compute probability of sharing cell with no governing point or governing point(s) of each type
+        frac_valid = np.sum(frac_sameregion) + frac_newregion
+        prob_newregion = frac_newregion / frac_valid
+        prob_sameregion= np.divide(frac_sameregion,frac_valid)
+        #compute probabilities of sampling each material type
+        return np.add( prob_sameregion, np.multiply(prob_newregion,self.prob) )
+
+    ## \brief Returns frequency of at least one psuedo-interface between two points in a Markovian-mixed medium
+    #
+    # \param[in] igovpt int, which governing point to compute the probability with regard to
+    # \returns float, probability of having at least one pseudo-interface
+    def freqAtLeastOnePseudointerface_Markovian(self,igovpt):
+        return 1.0 - np.exp( - self.govdists[igovpt] / self.lamc )

Core/3D/Geometry_Basepy.py

@@ -27,9 +27,17 @@ def defineGeometryBoundaries(self,xbounds=None,ybounds=None,zbounds=None):
         assert isinstance(xbounds,list) and isinstance(ybounds,list) and isinstance(zbounds,list)
         assert len(xbounds)==2 and len(ybounds)==2 and len(zbounds)==2
         assert isinstance(xbounds[0],float) and isinstance(xbounds[1],float) and xbounds[0]<=xbounds[1]
-        if self.Part.numDims<2: assert xbounds[0]==-np.inf and xbounds[1]==np.inf
-        assert isinstance(ybounds[0],float) and isinstance(ybounds[1],float) and ybounds[0]<=ybounds[1]
-        if self.Part.numDims<3: assert ybounds[0]==-np.inf and ybounds[1]==np.inf
+        if xbounds[0]==-np.inf and xbounds[1]==np.inf and ybounds[0]==-np.inf and ybounds[1]==np.inf:
+            self.numDims = 1
+        elif ybounds[0]==-np.inf and ybounds[1]==np.inf:
+            self.numDims = 2
+        else:
+            self.numDims = 3        
+
+        if hasattr(self, 'Part'):
+            assert self.numDims == self.Part.numDims
+        
+        assert isinstance(ybounds[0],float) and isinstance(ybounds[1],float) and ybounds[0]<=ybounds[1]        
         assert isinstance(zbounds[0],float) and isinstance(zbounds[1],float) and zbounds[0]<=zbounds[1]
         if self.GeomType=='Markovian':
             if not ( xbounds[0]==ybounds[0] and xbounds[0]==zbounds[0] and xbounds[1]==ybounds[1] and xbounds[1]==zbounds[1] and -xbounds[0]==xbounds[1] ):
@@ -40,8 +48,9 @@ def defineGeometryBoundaries(self,xbounds=None,ybounds=None,zbounds=None):
         self.zbounds = zbounds
 
         self.Volume = self.zbounds[1]-self.zbounds[0]
-        if self.Part.numDims>=2: self.Volume *= (self.xbounds[1]-self.xbounds[0])
-        if self.Part.numDims==3: self.Volume *= (self.ybounds[1]-self.ybounds[0])
+        
+        if self.numDims>=2: self.Volume *= (self.xbounds[1]-self.xbounds[0])
+        if self.numDims==3: self.Volume *= (self.ybounds[1]-self.ybounds[0])
 
 
     ## \brief Define boundary conditions
@@ -52,7 +61,7 @@ def defineBoundaryConditions(self,xBCs,yBCs,zBCs):
         assert isinstance(xBCs,list) and isinstance(yBCs,list) and isinstance(zBCs,list)
         assert len(xBCs)==2 and len(yBCs)==2 and len(zBCs)==2
 
-        if self.Part.numDims<2: assert xBCs[0]==None and xBCs[1]==None
+        if self.numDims<2: assert xBCs[0]==None and xBCs[1]==None
         else                  :
             assert xBCs[0]=='vacuum' or xBCs[0]=='reflective'
             assert xBCs[1]=='vacuum' or xBCs[1]=='reflective'
@@ -77,7 +86,8 @@ def defineBoundaryConditions(self,xBCs,yBCs,zBCs):
     # \returns sets cross sections
     def defineCrossSections(self,totxs=None,scatxs=None):
         assert isinstance(totxs,list)
-        self.nummats = len(totxs)
+        if not hasattr(self, "nummats"): self.nummats = len(totxs)
+        assert self.nummats==len(totxs)
         for i in range(0,self.nummats):
             assert isinstance(totxs[i],float) and totxs[i]>=0.0
         self.totxs = totxs
@@ -107,6 +117,8 @@ def associateRng(self,Rng):
     # \returns initializes self.Part
     def associatePart(self,Part):
         assert isinstance(Part,Particle)
+        if hasattr(self, "numDims"):
+            assert self.numDims == Part.numDims
         self.Part = Part
 
     ## \brief Tests for a leakage or boundary reflection event, flags the first, evaluates the second
@@ -169,7 +181,7 @@ def calculateDistanceToBoundary(self):
     # pseudocode: if woodcock, sample whether or not it is accepted. If it's not woodcock,
     # we just want to sample absorb or scatter. Will also do that if collision is accepted for woodcock.
     def evaluateCollision(self):
-        if self.trackingType == "Woodcock": # is there a way to tell that we are using woodcock tracking at this point in the code??? If not, I think I need two functions
+        if self.trackingType == "Woodcock":
             self.samplePoint() #define new point
             if self.Rng.rand()>self.totxs[self.CurrentMatInd]/self.Majorant: return 'reject' #accept or reject potential collision
         return 'absorb' if self.Rng.rand()>self.scatrat[self.CurrentMatInd] else 'scatter' #choose collision type
@@ -205,4 +217,4 @@ def defineMixingParams(self,laminf=None,prob=None):
     ## \brief Samples a new point according to the selected rules
     # To be provided by non-base class
     def samplePoint(self):
-        assert 1==0
+        assert 1==0

Core/3D/Geometry_BoxPoissonpy.py

@@ -1,5 +1,5 @@
 #!usr/bin/env python
-from Geometry_Basepy import Geometry_Base
+from Geometry_Voxelpy import Geometry_Voxel
 import numpy as np
 import bisect as bisect
 try:
@@ -15,7 +15,7 @@
 #
 # Class for multi-D Box-Poisson geometry 
 #
-class Geometry_BoxPoisson(Geometry_Base):
+class Geometry_BoxPoisson(Geometry_Voxel):
     def __init__(self):
         super(Geometry_BoxPoisson,self).__init__()
         self.flshowplot = False; self.flsaveplot = False
@@ -33,14 +33,14 @@ def __str__(self):
     def _initializeHistoryGeometryMemory(self):
         pass
 
-    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
+    ## \brief Samples hyperplanes to define geometry, voxelizes if selected
     #
-    # Note: Numpy arrays don't start blank and concatenate, where lists do, therefore plan to use lists, but convert to arrays using np.asarray(l) if needed to use array format
-    #
-    # \returns initializes self.MatInds and self.xBoundaries, self.yBoundaries, self.zBoundaries
+    # \returns initializes self.MatInds, self.xBoundaries, self.yBoundaries, self.zBoundaries; sets self.samplePoint; and voxelizes if selected
     def _initializeSampleGeometryMemory(self):
         self._sampleNumberOfHyperplanes()
         self._sampleHyperplaneLocations()
+        self.samplePoint = self.samplePoint_BoxPoisson
+        if self.flVoxelize: self.voxelizeGeometry()
 
     ## \brief Defines mixing parameters (correlation length and material probabilities)
     # Note: laminf and rhoB notation and formula following LarmierJQSRT2017 different mixing types paper
@@ -51,6 +51,8 @@ def _initializeSampleGeometryMemory(self):
     def defineMixingParams(self,laminf=None,prob=None):
         # Assert material chord lengths and slab length and store in object
         assert isinstance(laminf,float) and laminf>0.0
+        if not hasattr(self, "nummats"):
+            self.nummats = len(prob)
         assert isinstance(prob,list) and len(prob)==self.nummats
         for i in range(0,self.nummats): assert isinstance(prob[i],float) and prob[i]>=0.0 and prob[i]<=1.0
         assert self.isclose( np.sum(prob), 1.0 )
@@ -71,7 +73,6 @@ def _sampleNumberOfHyperplanes(self):
         self.yPInterfaces = self.Rng.poisson(yAve)
         self.zPInterfaces = self.Rng.poisson(zAve)
         self.MatInds = np.full( ( self.xPInterfaces+1 , self.yPInterfaces+1 , self.zPInterfaces+1 ),None )
-        #print(self.xPInterfaces,self.yPInterfaces,self.zPInterfaces, "planes sampled in x, y, and z")
 
 
     ## \brief sample location of pseudo-interfaces and constructs Box-Poisson geometry
@@ -98,24 +99,20 @@ def _sampleHyperplaneLocations(self):
         self.yBoundaries.sort()
         self.zBoundaries.sort()
 
-    ## \brief Samples a new point according to the selected rules
+    ## \brief Finds (or samples) material type at current location
     #
-    # \returns scattering ratio of cell material type
-    def samplePoint(self):
-        #find index of nearest point
-        x = self.Part.x
-        y = self.Part.y
-        z = self.Part.z
-
-        xIndex = bisect.bisect(self.xBoundaries,x) - 1
-        yIndex = bisect.bisect(self.yBoundaries,y) - 1
-        zIndex = bisect.bisect(self.zBoundaries,z) - 1
-        #check if material has been sampled. If not sample material type.
+    # \returns sets self.CurrentMatInd
+    def samplePoint_BoxPoisson(self):
+        #find indices of current box
+        xIndex = bisect.bisect(self.xBoundaries,self.Part.x) - 1
+        yIndex = bisect.bisect(self.yBoundaries,self.Part.y) - 1
+        zIndex = bisect.bisect(self.zBoundaries,self.Part.z) - 1
+        #check if material has been sampled in this box; if not, sample material type
+        if self.MatInds[xIndex,yIndex,zIndex]==None:
+            self.MatInds[xIndex,yIndex,zIndex] = int(self.Rng.choice(p=self.prob))
         self.CurrentMatInd = self.MatInds[xIndex,yIndex,zIndex]
-        if self.CurrentMatInd==None:
-            self.CurrentMatInd = int(self.Rng.choice(p=self.prob))
-            self.MatInds[xIndex,yIndex,zIndex] = self.CurrentMatInd
-        
+
+
     def samplePointsFast(self, points):
         #find index of nearest point 
         xinds = np.searchsorted(self.xBoundaries, points[:,0], side='right')-1
@@ -144,11 +141,7 @@ def samplePointsFast(self, points):
             xi = np.int16(xi)
             yi = np.int16(yi)
             zi = np.int16(zi)   
-            #uniques = np.array(list(uniques.keys()))
-            #uniques = np.unique(unassigned_inds, axis=0)
             new_assignments = self.Rng.RngObj.choice(self.nummats, size=len(uniques), p=self.prob)            
-            #self.MatInds[uniques[:,0], uniques[:,1], uniques[:,2]] = new_assignments
             self.MatInds[xi, yi, zi] = new_assignments
             matTypes = self.MatInds[xinds, yinds, zinds]
-        return matTypes.astype(np.int16)
-
+        return matTypes.astype(np.int16)

Core/3D/Geometry_CLSpy.py

@@ -10,8 +10,6 @@
 ## \brief Multi-D CLS geometry class
 # \author Dominic Lioce, [email protected], [email protected]
 #
-# A good chunk of this class is copied from the CoPS geometry class. Only a few things are different
-#
 class Geometry_CLS(Geometry_Base,ClassTools,MarkovianInputs):
     def __init__(self):
         super(Geometry_CLS,self).__init__()