import from dgcode ncrystaldev git: candidate for version 3.9.6

CHANGELOG

@@ -1,3 +1,10 @@
+v3.9.6 2024-08-30
+     * Add several new convenience methods to MiniMC results and histograms.
+     * Fix bug in MiniMC histogram error propagation.
+     * NCMATComposer gets new expert-only methods for adding @CUSTOM_ sections
+       or raw text to NCMAT data.
+     * NCMATComposer.from_hfg sets plotlabel from title.
+
 v3.9.5 2024-08-28
      * Make ncrystal_hfg2ncmat functionality available in Python API as well in
        the new NCrystal.hfg2ncmat module. Also add NCMATComposer.from_hfg(..)

CMakeLists.txt

@@ -86,7 +86,7 @@ endif()
 
 cmake_policy( SET CMP0048 NEW )#Not sure if this is really needed
 
-project( NCrystal VERSION 3.9.5 ${_project_metadata} )
+project( NCrystal VERSION 3.9.6 ${_project_metadata} )
 
 unset( _project_metadata )
 

NCrystal/__init__.py

@@ -54,7 +54,7 @@
 
 #NB: Synchronize meta-data below with fields in setup.py+template_setup.py.in meta data:
 __license__ = "Apache 2.0, http://www.apache.org/licenses/LICENSE-2.0"
-__version__ = '3.9.5'
+__version__ = '3.9.6'
 __status__ = "Production"
 __author__ = "NCrystal developers (Thomas Kittelmann, Xiao Xiao Cai)"
 __copyright__ = "Copyright 2015-2024 %s"%__author__

NCrystal/_cli_hfg2ncmat.py

@@ -72,13 +72,13 @@ def _parseArgs( default_debye_temp ):
     parser = argparse.ArgumentParser(description=descr,
                                      formatter_class=RawTextHelpFormatter)
     parser.add_argument("--output",'-o',default='autogen.ncmat',type=str,
-                        help=(f"Output file name (defaults to autogen.ncmat)."
+                        help=("Output file name (defaults to autogen.ncmat)."
                               " Can be stdout."))
     parser.add_argument('--force',action='store_true',
-                        help=(f"Will overwrite existing file "
+                        help=("Will overwrite existing file "
                               "if it already exists."))
     parser.add_argument("--spec",'-s',metavar='SPEC',type=str,required=True,
-                        help=f"Hydrogen binding specification (see above).")
+                        help="Hydrogen binding specification (see above).")
     parser.add_argument("--formula",'-f',metavar='FORMULA',
                         type=str,
                         required=True,
@@ -87,7 +87,7 @@ def _parseArgs( default_debye_temp ):
     parser.add_argument("--density",'-d',metavar='DENSITY',
                         type=float,
                         required=True,
-                        help=f"Material density in g/cm3.")
+                        help="Material density in g/cm3.")
     parser.add_argument("--debyetemp",metavar='VALUE',
                         type=float,
                         default=default_debye_temp,
@@ -98,7 +98,7 @@ def _parseArgs( default_debye_temp ):
                               ' comment near top of output file). Use \\n '
                               'for line-breaks.'))
     parser.add_argument('--notrim',action='store_true',
-                        help=f"No trimming of resulting VDOS curve.")
+                        help="No trimming of resulting VDOS curve.")
     args=parser.parse_args()
     return args
 

NCrystal/_mmc.py

@@ -111,6 +111,96 @@ def clone( self, rebin_factor = 1 ):
                 c.rebin( rebin_factor )
             return c
 
+        def integrate( self, xlow, xhigh, tolerance = 1e-5 ):
+            """
+            Returns integrated contents of the histogram over the area
+            [xlow,xhigh] along with the error of that value in a tuple
+            (content,error).
+
+            This is done translating xlow and xhigh to exact bin edges and then
+            calling integrate_bins. If that is not possible within the
+            tolerance, an exception is raised.
+
+            """
+            if not ( xhigh >= xlow ):
+                from .exceptions import NCBadInput
+                raise NCBadInput('Invalid integration range requested.')
+
+            bw = self.binwidth
+            def _findedge(x):
+                if x <= self.__xmin:
+                    return 0
+                if x >= self.__xmax:
+                    return self.nbins
+                r = ( x - self.__xmin ) / bw
+                ir = int(r+0.5)
+                if abs(r-ir) > tolerance:
+                    from .exceptions import NCBadInput
+                    raise NCBadInput(f'Value {x} does not correspond exactly'
+                                     ' to a bin edge within the tolerance.')
+                return ir
+            e_low = _findedge(xlow)
+            e_high = _findedge(xhigh)
+            if e_low == e_high:
+                return ( 0.0, 0.0 )
+            assert e_low >= 0
+            assert e_low < e_high < self.nbins
+            return self.integrate_bins( e_low, e_high - 1 )
+
+        def integrate_bins( self, bin_low = None, bin_up = None ):
+            """
+            Returns integrated contents of the bins [bin_low,bin_up[ along with
+            the error of that value in a tuple (content,error).
+
+            If bin_low is None the integration will start at the first bin and
+            include the underflow bin.
+
+            If bin_up is None the integration will end at the last bin and
+            include the overflow bin.
+            """
+
+            add_overflow, add_underflow = False, False
+            if bin_low is None:
+                add_underflow = True
+                bin_low = 0
+                underflow_c = self.stats.get('underflow')
+                underflow_e2 = self.stats.get('underflow_errorsq')
+                if bool(underflow_c is None) != bool(underflow_e2 is None):
+                    from .exceptions import NCBadInput
+                    raise NCBadInput('Inconsistent underflow info')
+                if underflow_c is None:
+                    add_underflow = False
+
+            if bin_up is None:
+                add_overflow = True
+                bin_up = self.__nbins
+                overflow_c = self.stats.get('overflow')
+                overflow_e2 = self.stats.get('overflow_errorsq')
+                if bool(overflow_c is None) != bool(overflow_e2 is None):
+                    from .exceptions import NCBadInput
+                    raise NCBadInput('Inconsistent overflow info')
+                if overflow_c is None:
+                    add_overflow = False
+
+            bin_low, bin_up = int(bin_low), int(bin_up)
+            if bin_up < bin_low or bin_low<0 or bin_up > self.__nbins:
+                from .exceptions import NCBadInput
+                raise NCBadInput('Invalid bin range requested')
+            content_integral = self.__y[bin_low:bin_up].sum()
+            if add_underflow:
+                content_integral += underflow_c
+            if add_overflow:
+                content_integral += overflow_c
+            if self.__yerrsq is None:
+                #unweighted, just base erros on contents:
+                return ( content_integral, _np.sqrt(content_integral) )
+            errorsq_integral = self.__yerrsq[bin_low:bin_up].sum()
+            if add_underflow:
+                errorsq_integral += underflow_e2
+            if add_overflow:
+                errorsq_integral += overflow_e2
+            return ( content_integral, _np.sqrt(errorsq_integral) )
+
         def add_contents( self, other_hist ):
             o = other_hist
             assert self.__xmin == o.__xmin
@@ -124,12 +214,12 @@ def add_contents( self, other_hist ):
                 if o.__yerrsq is None:
                     pass#done
                 else:
-                    self.__yerrsq = o.__yerrsq
+                    self.__yerrsq = self.__y + o.__yerrsq
             else:
                 if o.__yerrsq is None:
-                    self.__yerrsq = o.__y
+                    self.__yerrsq += o.__y
                 else:
-                    self.__yerrsq = o.__yerrsq
+                    self.__yerrsq += o.__yerrsq
 
         def rebin( self, rebin_factor ):
             assert self.__nbins % rebin_factor == 0
@@ -151,11 +241,15 @@ def stats( self ):
         @property
         def errors( self ):
             if self.__yerr is None:
-                self.__yerr = _np.sqrt( self.__yerrsq
-                                        if self.__yerrsq is not None
-                                        else self.__y )
+                self.__yerr = _np.sqrt( self.errors_squared )
             return self.__yerr
 
+        @property
+        def errors_squared( self ):
+            return ( self.__yerrsq
+                     if self.__yerrsq is not None
+                     else self.__y )
+
         @property
         def content( self ):
             return self.__y
@@ -286,6 +380,33 @@ def __init__( self,
         def histograms( self ):
             return self.__hists
 
+        @property
+        def histogram_main( self ):
+            return self.__hists[0]
+
+        @property
+        def histogram_breakdown( self ):
+            return dict( (h.title,h) for h in self.__hists[1:] )
+
+        def histogram_sum( self, *, select=None, exclude=None ):
+            if isinstance(exclude,str):
+                exclude=[exclude]
+            if isinstance(select,str):
+                select=[select]
+            hl = self.__hists[1:]
+            if not exclude and not select:
+                return self.histogram_main
+            if select:
+                hl = [ h for h in hl if h.title in select ]
+            if exclude:
+                hl = [ h for h in hl if h.title not in exclude ]
+            if len(hl) <= 1:
+                return hl[0] if hl else None
+            h = hl[0].clone()
+            for o in hl[1:]:
+                h.add_contents( o )
+            return h
+
         @property
         def cfgstr( self ):
             return self.__cfgstr
@@ -606,8 +727,6 @@ def quick_diffraction_pattern( cfgstr, *,
     mfp_scatter = 1/macroxs_scatter if macroxs_scatter else float('inf')
 
     sphere_diameter = _parse_length(material_thickness, mfp=mfp_scatter)
-    #print(f'mfp_scatter = {mfp_scatter/unit_cm:g} cm')
-    #print(f'sphere_diameter = {sphere_diameter/unit_cm:g} cm')
 
     def simfct( n, cfgstr ):
         import time