Version 1.1.0 - tallies, var. dec./opt., pre-samp. pts, sph. inc.

- Refactoring of flux capabilities to be easier on users and developers.
- New features of variance deconvolution (ClementsJQSRT2024) and optimization for mean quantities (OlsonANS2018)--these required enabling multiple histories per sample in 3D and multiple samples per partition for statistical estimates.
- New feature of pre-sampled points on a grid with grid lookups for 1D CoPS.
- Refactoring and new features for spherical inclusion geometry to enable faster realization generation and material lookups.

Author: aaronjeffreyolson

Core/1D/MonteCarloParticleSolverpy.py

@@ -184,7 +184,6 @@ def populateMarkovRealization(self,totxs=None,lam=None,s=None,scatxs=None,NaryTy
         b0 = [0]*self.nummats; b0.append(1)
         b=np.transpose(np.array(b0))
         self.chordfreq = list( np.linalg.solve(np.transpose(A).dot(A), np.transpose(A).dot(b)) )
-        #########print('self.chordfreq:',self.chordfreq)
         self.prob = []
         for i in range(0,self.nummats):
             self.prob.append( self.chordfreq[i]*self.lam[i] / np.dot(self.chordfreq,self.lam) )
@@ -203,8 +202,6 @@ def populateMarkovRealization(self,totxs=None,lam=None,s=None,scatxs=None,NaryTy
         else:
             while True: #sample next material either by abundance or with equal probability
                 local_probs = self.NaryType[curmat]
-#                local_probs[curmat] = 0.0                                      #give probability of zero to current material
-#                local_probs = np.divide( local_probs, np.sum(local_probs) )    #normalize probabilities
                 curmat = int( self.Rng.choice( self.nummats, p=local_probs ) ) #select material from any but current
 
                 self.mattype.append(curmat)
@@ -222,7 +219,8 @@ def populateMarkovRealization(self,totxs=None,lam=None,s=None,scatxs=None,NaryTy
     # \param[in] lam list of average material chord lengths
     # \param[in] s float, the length of the slab
     # \param     scatxs list, optional, list of scattering cross section values, if specified, absorption cross sections and scattering ratios will also be computed and stored as attributes
-    def populateMarkovRealizationPseudo(self,totxs=None,lam=None,s=None,scatxs=None):
+    def populateMarkovRealizationPseudo(self,totxs=None,lam=None,s=None,scatxs=None,NaryType='volume-fraction'):
+        assert NaryType=='volume-fraction' #the pseudo-interface approach is only a valid approach for 'volume-fraction'-type Markovian mixing; for other types of Markovian mixing, use the chord-based sampler
         assert isinstance(totxs,list)
         self.nummats = len(totxs)
         for i in range(0,self.nummats): assert isinstance(totxs[i],float) and totxs[i]>=0.0

Core/1D/OneDSlabpy.py

@@ -5,7 +5,6 @@
 from Particlepy import Particle
 from ClassToolspy import ClassTools
 import numpy as np
-import matplotlib.pyplot as plt
 
 ## \brief Multi-D base geometry class - useless by itself, but other multi-D geometry classes inherit from it
 # \author Aaron Olson, [email protected], [email protected]
@@ -40,9 +39,9 @@ def defineGeometryBoundaries(self,xbounds=None,ybounds=None,zbounds=None):
         self.ybounds = ybounds
         self.zbounds = zbounds
 
-        self.Volume = self.zbounds[1]-zbounds[0]
-        if self.Part.numDims>=2: self.Volume *= (self.xbounds[1]-xbounds[0])
-        if self.Part.numDims==3: self.Volume *= (self.ybounds[1]-ybounds[0])
+        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])
 
 
     ## \brief Define boundary conditions
@@ -188,20 +187,21 @@ def _generateRandomPoint(self):
 
 
 
+    ## \brief Initialize geometry parameters at start of each history
+    # \returns nothing
+    def _initializeHistoryGeometryMemory(self):
+        assert 1==0
 
-
-
+    ## \brief Initialize geometry parameters at start of each sample
+    # \returns nothing
+    def _initializeSampleGeometryMemory(self):
+        assert 1==0
 
     ## \brief Define mixing parameters (correlation length and material probabilities)
     # To be provided by non-base class
     def defineMixingParams(self,laminf=None,prob=None):
         assert 1==0
 
-    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
-    # To be provided by non-base class
-    def initializeGeometryMemory(self):
-        assert 1==0
-
     ## \brief Samples a new point according to the selected rules
     # To be provided by non-base class
     def samplePoint(self):

Core/1D/SpecialMonteCarloDriverspy.py

@@ -1,9 +1,7 @@
 #!usr/bin/env python
 from Geometry_Basepy import Geometry_Base
 import numpy as np
-import scipy.spatial
 import bisect as bisect
-import collections
 try:
     import pandas as pd
 except:
@@ -23,11 +21,27 @@ def __init__(self):
         self.flshowplot = False; self.flsaveplot = False
         self.kdtrees = False
         self.GeomType = 'BoxPoisson'
+        self.abundanceModel = 'ensemble'
 
     def __str__(self):
         return str(self.__dict__)
 
 
+    ## \brief Some geometries will need to overwrite this method definition, others won't.
+    #
+    # \returns nothing
+    def _initializeHistoryGeometryMemory(self):
+        pass
+
+    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
+    #
+    # 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
+    def _initializeSampleGeometryMemory(self):
+        self._sampleNumberOfHyperplanes()
+        self._sampleHyperplaneLocations()
+
     ## \brief Defines mixing parameters (correlation length and material probabilities)
     # Note: laminf and rhoB notation and formula following LarmierJQSRT2017 different mixing types paper
     #
@@ -45,16 +59,6 @@ def defineMixingParams(self,laminf=None,prob=None):
         self.rhoB   = 2.0 / 3.0 / self.laminf
 
 
-    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
-    #
-    # 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
-    def initializeGeometryMemory(self):
-        self._sampleNumberOfHyperplanes()
-        self._sampleHyperplaneLocations()
-        
-
     ## \brief compute number of pseudo-interfaces
     #
     # \computes number of pseudo-interfaces in x, y, and z direction for defined geometry 

Core/3D/Geometry_Basepy.py

@@ -2,8 +2,6 @@
 import sys
 from Geometry_Basepy import Geometry_Base
 sys.path.append('/../../Classes/Tools')
-from RandomNumberspy import RandomNumbers
-from Particlepy import Particle
 from ClassToolspy import ClassTools
 from MarkovianInputspy import MarkovianInputs
 import numpy as np
@@ -18,10 +16,22 @@ class Geometry_CLS(Geometry_Base,ClassTools,MarkovianInputs):
     def __init__(self):
         super(Geometry_CLS,self).__init__()
         self.GeomType = 'CLS'
+        self.abundanceModel = 'ensemble'
 
     def __str__(self):
         return str(self.__dict__)
 
+    ## \brief Some geometries will need to overwrite this method definition, others won't.
+    #
+    # \returns nothing
+    def _initializeHistoryGeometryMemory(self):
+        pass
+
+    ## \brief Some geometries will need to overwrite this method definition, others won't.
+    #
+    # \returns nothing
+    def _initializeSampleGeometryMemory(self):
+        pass
 
     ## \brief Defines mixing parameters (chord lengths and probabilities)
     #
@@ -52,12 +62,6 @@ def defineCLSGeometryType(self,CLSAlg="CLS",fl1DEmulation=False):
         assert isinstance(fl1DEmulation,bool)
         self.fl1DEmulation = fl1DEmulation
 
-    ## \brief Initializes geometry memory; CLS doesn't have any so only initialize if LRP
-    def initializeGeometryMemory(self):
-        if self.CLSAlg == "LRP": 
-            self.dip = 0.0
-            self.dim = 0.0
-
     ## \brief Samples material index at start of particle transport based on volume fractions
     #
     # \returns updates material index and point attributes

Core/3D/Geometry_BoxPoissonpy.py

@@ -2,8 +2,6 @@
 import sys
 from Geometry_Basepy import Geometry_Base
 sys.path.append('/../../Classes/Tools')
-from RandomNumberspy import RandomNumbers
-from Particlepy import Particle
 from ClassToolspy import ClassTools
 from MarkovianInputspy import MarkovianInputs
 import numpy as np
@@ -23,10 +21,30 @@ def __init__(self):
         self.flshowplot = False; self.flsaveplot = False
         self.flCollectPoints = False
         self.GeomType = 'CoPS'
+        self.abundanceModel = 'ensemble'
 
     def __str__(self):
         return str(self.__dict__)
 
+    ## \brief Initializes short-term points and material indices for CoPS (starts history)--erases memory
+    #
+    # 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.RecentPoints and self.RecentMatInds
+    def _initializeHistoryGeometryMemory(self):
+        self.RecentPoints    = []
+        self.RecentMatInds   = []
+
+    ## \brief Initializes long-term memory points and material indices for CoPS (starts cohort)--erases memory
+    #
+    # 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.LongTermPoints and self.LongTermMatInds
+    def _initializeSampleGeometryMemory(self):
+        self.LongTermPoints  = []
+        self.LongTermMatInds = []
 
     ## \brief Defines mixing parameters (chord lengths and probabilities)
     #
@@ -41,20 +59,6 @@ def defineMixingParams(self,lam=None):
         # Solve and store material probabilities and correlation length
         self.solveNaryMarkovianParamsBasedOnChordLengths(self.lam)
 
-    ## \brief Initializes points and material indices--erases any memory
-    #
-    # 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
-    # "LongTerm" refers to storage of point location and material index until the end of a cohort
-    # "Recent" refers to storage of the N most recent points not committed to long term memory
-    #
-    # \returns initializes self.LongTermPoints, self.LongTermMatInds, self.RecentPoints, and self.RecentMatInds
-    def initializeGeometryMemory(self):
-        self.LongTermPoints  = []
-        self.LongTermMatInds = []
-        self.RecentPoints    = []
-        self.RecentMatInds   = []
-
     ## \brief Samples material index of new point, and stores or forgets new point, according to selected rules
     #
     # Calls other method to solve conditional probabilities of material index at new point based on previously
@@ -464,25 +468,4 @@ def selectWhetherPrintPoints(self,flCollectPoints=False,pointsNpyFileName='CoPSP
         self.flCollectPoints = flCollectPoints
         if self.flCollectPoints:
             assert isinstance(pointsNpyFileName,str)
-            self.pointsNpyFileName = pointsNpyFileName
-
-    ## \brief Collect CoPS output points and store in accumulating list
-    #
-    # \param[in] ipart int, Particle number (same as realization number so far)
-    # \returns creates or contributes to master list of points self.OldPointData
-    def appendNewPointsToCollection(self,ipart):
-        if len(self.LongTermMatInds)>0:
-            #Create array of particle number, combine particle number, material indices, and x, y, z coordinates into one array
-            NewPointData = np.concatenate( ( np.stack( (np.full(len(self.LongTermMatInds),ipart), np.array(self.LongTermMatInds)),  axis=-1 ) , np.array(self.LongTermPoints) ) , axis=1)
-            #Add to master list of points defined in all realizations (or create that list)
-            if hasattr(self,'OldPointData'): self.OldPointData = np.concatenate( (self.OldPointData,NewPointData), axis=0)
-            else                           : self.OldPointData = NewPointData
-
-
-    ## \brief If flag set to collect points, print collected points to npy file
-    #
-    # \returns prints data to file
-    def printCollectedPointsToNumpyFile(self):
-        if self.flCollectPoints: np.save( self.pointsNpyFileName , self.OldPointData )
-
-
+            self.pointsNpyFileName = pointsNpyFileName

Core/3D/Geometry_CLSpy.py

@@ -19,10 +19,26 @@ def __init__(self):
         super(Geometry_Markovian,self).__init__()
         self.flshowplot = False; self.flsaveplot = False
         self.GeomType = 'Markovian'
+        self.abundanceModel = 'ensemble'
 
     def __str__(self):
         return str(self.__dict__)
-    
+
+    ## \brief Some geometries will need to overwrite this method definition, others won't.
+    #
+    # \returns nothing
+    def _initializeHistoryGeometryMemory(self):
+        pass
+
+    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
+    #
+    # 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
+    def _initializeSampleGeometryMemory(self):
+        self._sampleNumberOfHyperplanes()
+        self._sampleHyperplaneLocations()
+
     ## \brief Defines mixing parameters (correlation length and material probabilities)
     # Note: laminf and rhoP notation following LarmierJQSRT2017 different mixing types paper
     #
@@ -40,15 +56,6 @@ def defineMixingParams(self,laminf=None,prob=None):
         self.rhoP   = 1.0 / self.laminf
 
 
-    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
-    #
-    # 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
-    def initializeGeometryMemory(self):
-        self._sampleNumberOfHyperplanes()
-        self._sampleHyperplaneLocations()
-
     ## \brief compute average number of hyper-planes
     #
     # \computes average number of hyper-planes used to sample number of planes (Poisson)