diff --git a/demos/MPI/cornell_box_hybrid.py b/demos/MPI/cornell_box_hybrid.py new file mode 100644 index 00000000..9fe547c1 --- /dev/null +++ b/demos/MPI/cornell_box_hybrid.py @@ -0,0 +1,213 @@ +""" +Cornell Box MPI Demo +==================== + +This demo renders a variant of the classic Cornell Box scene. + +For the original Cornell Box see: + + http://www.graphics.cornell.edu/online/box/data.html + +The wall colours and light spectrum used in this demo are the values measured +for the physical Cornell Box. + +The MPI (message passing interface) render engine is used for this demo. +To run, use `mpirun -n python cornell_box_mpi.py`. An +MPI library must be installed on your system, and the mpi4py Python +package must be installed in the Python environment. + +There are some small differences between this demo and the shared-memory +equilvalent cornell_box.py, due to details of the MPI implementation meaning +we need to communicate the results of one render pass to all the worker +processes before the next pass. Compare the last section of the two demo +scripts to see how this is done. +""" + +from numpy import array + +from raysect.primitive import Sphere, Box +from raysect.optical import World, Node, translate, rotate, Point3D +from raysect.optical.material import Lambert, UniformSurfaceEmitter +from raysect.optical.library import InterpolatedSF, schott +from raysect.optical.observer import PinholeCamera +from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, PowerPipeline2D +from raysect.optical.observer import RGBAdaptiveSampler2D +from raysect.core.workflow import HybridEngine, MulticoreEngine + + +# define reflectivity for box surfaces +wavelengths = array( + [400, 404, 408, 412, 416, 420, 424, 428, 432, 436, 440, 444, 448, 452, 456, 460, 464, 468, 472, 476, 480, 484, 488, + 492, 496, 500, 504, 508, 512, 516, 520, 524, 528, 532, 536, 540, 544, 548, 552, 556, 560, 564, 568, 572, 576, 580, + 584, 588, 592, 596, 600, 604, 608, 612, 616, 620, 624, 628, 632, 636, 640, 644, 648, 652, 656, 660, 664, 668, 672, + 676, 680, 684, 688, 692, 696, 700]) + +white = array( + [0.343, 0.445, 0.551, 0.624, 0.665, 0.687, 0.708, 0.723, 0.715, 0.71, 0.745, 0.758, 0.739, 0.767, 0.777, 0.765, + 0.751, 0.745, 0.748, 0.729, 0.745, 0.757, 0.753, 0.75, 0.746, 0.747, 0.735, 0.732, 0.739, 0.734, 0.725, 0.721, + 0.733, 0.725, 0.732, 0.743, 0.744, 0.748, 0.728, 0.716, 0.733, 0.726, 0.713, 0.74, 0.754, 0.764, 0.752, 0.736, + 0.734, 0.741, 0.74, 0.732, 0.745, 0.755, 0.751, 0.744, 0.731, 0.733, 0.744, 0.731, 0.712, 0.708, 0.729, 0.73, + 0.727, 0.707, 0.703, 0.729, 0.75, 0.76, 0.751, 0.739, 0.724, 0.73, 0.74, 0.737]) + +green = array( + [0.092, 0.096, 0.098, 0.097, 0.098, 0.095, 0.095, 0.097, 0.095, 0.094, 0.097, 0.098, 0.096, 0.101, 0.103, 0.104, + 0.107, 0.109, 0.112, 0.115, 0.125, 0.14, 0.16, 0.187, 0.229, 0.285, 0.343, 0.39, 0.435, 0.464, 0.472, 0.476, 0.481, + 0.462, 0.447, 0.441, 0.426, 0.406, 0.373, 0.347, 0.337, 0.314, 0.285, 0.277, 0.266, 0.25, 0.23, 0.207, 0.186, + 0.171, 0.16, 0.148, 0.141, 0.136, 0.13, 0.126, 0.123, 0.121, 0.122, 0.119, 0.114, 0.115, 0.117, 0.117, 0.118, 0.12, + 0.122, 0.128, 0.132, 0.139, 0.144, 0.146, 0.15, 0.152, 0.157, 0.159]) + +red = array( + [0.04, 0.046, 0.048, 0.053, 0.049, 0.05, 0.053, 0.055, 0.057, 0.056, 0.059, 0.057, 0.061, 0.061, 0.06, 0.062, 0.062, + 0.062, 0.061, 0.062, 0.06, 0.059, 0.057, 0.058, 0.058, 0.058, 0.056, 0.055, 0.056, 0.059, 0.057, 0.055, 0.059, + 0.059, 0.058, 0.059, 0.061, 0.061, 0.063, 0.063, 0.067, 0.068, 0.072, 0.08, 0.09, 0.099, 0.124, 0.154, 0.192, + 0.255, 0.287, 0.349, 0.402, 0.443, 0.487, 0.513, 0.558, 0.584, 0.62, 0.606, 0.609, 0.651, 0.612, 0.61, 0.65, 0.638, + 0.627, 0.62, 0.63, 0.628, 0.642, 0.639, 0.657, 0.639, 0.635, 0.642]) + +white_reflectivity = InterpolatedSF(wavelengths, white) +red_reflectivity = InterpolatedSF(wavelengths, red) +green_reflectivity = InterpolatedSF(wavelengths, green) + +# define light spectrum +light_spectrum = InterpolatedSF(array([400, 500, 600, 700]), array([0.0, 8.0, 15.6, 18.4])) + +# set-up scenegraph +world = World() + +# enclosing box +enclosure = Node(world) + +e_back = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, 0, 1) * rotate(0, 0, 0), + material=Lambert(white_reflectivity)) + +e_bottom = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, -1, 0) * rotate(0, -90, 0), + # material=m) + material=Lambert(white_reflectivity)) + +e_top = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, 1, 0) * rotate(0, 90, 0), + material=Lambert(white_reflectivity)) + +e_left = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(1, 0, 0) * rotate(-90, 0, 0), + material=Lambert(red_reflectivity)) + +e_right = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(-1, 0, 0) * rotate(90, 0, 0), + material=Lambert(green_reflectivity)) + +# ceiling light +light = Box(Point3D(-0.4, -0.4, -0.01), Point3D(0.4, 0.4, 0.0), + parent=enclosure, + transform=translate(0, 1, 0) * rotate(0, 90, 0), + material=UniformSurfaceEmitter(light_spectrum, 2)) + +# alternate light #1 +# light = Box(Point3D(-0.4, -0.4, -0.01), Point3D(0.4, 0.4, 0.0), +# parent=enclosure, +# transform=translate(0, 1, 0) * rotate(0, 90, 0), +# material=UniformSurfaceEmitter(d65_white, 2)) + +# alternate light #2 +# back_light = Sphere(0.1, +# parent=enclosure, +# transform=translate(0.80, -0.85, 0.80)*rotate(0, 0, 0), +# material=UniformSurfaceEmitter(light_spectrum, 10.0)) + +# objects in enclosure +box = Box(Point3D(-0.4, 0, -0.4), Point3D(0.3, 1.4, 0.3), + parent=world, + transform=translate(0.4, -1 + 1e-6, 0.4)*rotate(30, 0, 0), + material=schott("N-BK7")) + +sphere = Sphere(0.4, + parent=world, + transform=translate(-0.4, -0.6 + 1e-6, -0.4)*rotate(0, 0, 0), + material=schott("N-BK7")) + + +filter_red = InterpolatedSF([100, 650, 660, 670, 680, 800], [0, 0, 1, 1, 0, 0]) +filter_green = InterpolatedSF([100, 530, 540, 550, 560, 800], [0, 0, 1, 1, 0, 0]) +filter_blue = InterpolatedSF([100, 480, 490, 500, 510, 800], [0, 0, 1, 1, 0, 0]) + +# create and setup the camera +power_unfiltered = PowerPipeline2D(display_unsaturated_fraction=0.96, name="Unfiltered") +power_unfiltered.display_update_time = 15 + +power_green = PowerPipeline2D(filter=filter_green, display_unsaturated_fraction=0.96, name="Green Filter") +power_green.display_update_time = 15 + +power_red = PowerPipeline2D(filter=filter_red, display_unsaturated_fraction=0.96, name="Red Filter") +power_red.display_update_time = 15 + +rgb = RGBPipeline2D(display_unsaturated_fraction=0.96, name="sRGB") + +bayer = BayerPipeline2D(filter_red, filter_green, filter_blue, display_unsaturated_fraction=0.96, name="Bayer Filter") +bayer.display_update_time = 15 + +pipelines = [rgb, power_unfiltered, power_green, power_red, bayer] + +sampler = RGBAdaptiveSampler2D(rgb, ratio=10, fraction=0.2, min_samples=500, cutoff=0.01) + + +camera = PinholeCamera((1024, 1024), parent=world, transform=translate(0, 0, -3.3) * rotate(0, 0, 0), pipelines=pipelines) +camera.frame_sampler = sampler +camera.spectral_rays = 1 +camera.spectral_bins = 15 +camera.pixel_samples = 250 +camera.ray_importance_sampling = True +camera.ray_important_path_weight = 0.25 +camera.ray_max_depth = 500 +camera.ray_extinction_min_depth = 3 +camera.ray_extinction_prob = 0.01 + +# Speedups for testing purposes. +# Fewer camera pixels for faster runtime per pass. +# camera.pixels = (256, 256) +# Less strict cutoff to finish with fewer passes. +sampler.cutoff = 0.05 + +# Get the available parallelism for this MPI process automatically, +# or hard code it if your scheduler isn't supported. +nworkers = HybridEngine.estimate_subworker_count() +# nworkers = 16 +camera.render_engine = HybridEngine(MulticoreEngine(nworkers)) + +# Don't make plots during the render as MPI is typically run non-interactively. +for pipeline in camera.pipelines: + pipeline.display_progress = False + +# Worker processes do not have accurate statistics, only the root process has +# all the render results. So don't bother outputting statistics on the workers. +if camera.render_engine.rank != 0: + camera.quiet = True + +# start ray tracing. +print(f"Starting ray tracing on rank {camera.render_engine.rank}") +render_pass = 1 +while not camera.render_complete: + if camera.render_engine.rank == 0: + print(f"Rendering pass {render_pass}...") + camera.observe() + # Rank 0 processes the ray tracing results so is the only process which + # knows the true progress of the render. In order to know when the render + # is complete, we need to ensure the frame sampler on each of the worker + # processes has an up-to-date copy of the pipeline used for sampling. + # We also need to use a different variable name for the received broadcast + # object else it is not broadcast properly. + root_rgb = camera.render_engine.comm.bcast(rgb, root=0) + camera.frame_sampler.pipeline = root_rgb + camera.render_engine.comm.Barrier() + render_pass += 1 +print(f"Finished ray tracing on rank {camera.render_engine.rank}") + +# Again, only rank 0 has all the ray tracing results so is the only one which +# can produce a correct image. +if camera.render_engine.rank == 0: + rgb.save('CornellHybrid_rgb.png') diff --git a/demos/MPI/cornell_box_mpi.py b/demos/MPI/cornell_box_mpi.py new file mode 100644 index 00000000..950b5708 --- /dev/null +++ b/demos/MPI/cornell_box_mpi.py @@ -0,0 +1,210 @@ +""" +Cornell Box MPI Demo +==================== + +This demo renders a variant of the classic Cornell Box scene. + +For the original Cornell Box see: + + http://www.graphics.cornell.edu/online/box/data.html + +The wall colours and light spectrum used in this demo are the values measured +for the physical Cornell Box. + +The MPI (message passing interface) render engine is used for this demo. +To run, use `mpirun -n python cornell_box_mpi.py`. An +MPI library must be installed on your system, and the mpi4py Python +package must be installed in the Python environment. + +There are some small differences between this demo and the shared-memory +equilvalent cornell_box.py, due to details of the MPI implementation meaning +we need to communicate the results of one render pass to all the worker +processes before the next pass. Compare the last section of the two demo +scripts to see how this is done. +""" + +from numpy import array + +from raysect.primitive import Sphere, Box +from raysect.optical import World, Node, translate, rotate, Point3D +from raysect.optical.material import Lambert, UniformSurfaceEmitter +from raysect.optical.library import InterpolatedSF, schott +from raysect.optical.observer import PinholeCamera +from raysect.optical.observer import RGBPipeline2D, BayerPipeline2D, PowerPipeline2D +from raysect.optical.observer import RGBAdaptiveSampler2D +from raysect.core.workflow import MPIEngine + + +# define reflectivity for box surfaces +wavelengths = array( + [400, 404, 408, 412, 416, 420, 424, 428, 432, 436, 440, 444, 448, 452, 456, 460, 464, 468, 472, 476, 480, 484, 488, + 492, 496, 500, 504, 508, 512, 516, 520, 524, 528, 532, 536, 540, 544, 548, 552, 556, 560, 564, 568, 572, 576, 580, + 584, 588, 592, 596, 600, 604, 608, 612, 616, 620, 624, 628, 632, 636, 640, 644, 648, 652, 656, 660, 664, 668, 672, + 676, 680, 684, 688, 692, 696, 700]) + +white = array( + [0.343, 0.445, 0.551, 0.624, 0.665, 0.687, 0.708, 0.723, 0.715, 0.71, 0.745, 0.758, 0.739, 0.767, 0.777, 0.765, + 0.751, 0.745, 0.748, 0.729, 0.745, 0.757, 0.753, 0.75, 0.746, 0.747, 0.735, 0.732, 0.739, 0.734, 0.725, 0.721, + 0.733, 0.725, 0.732, 0.743, 0.744, 0.748, 0.728, 0.716, 0.733, 0.726, 0.713, 0.74, 0.754, 0.764, 0.752, 0.736, + 0.734, 0.741, 0.74, 0.732, 0.745, 0.755, 0.751, 0.744, 0.731, 0.733, 0.744, 0.731, 0.712, 0.708, 0.729, 0.73, + 0.727, 0.707, 0.703, 0.729, 0.75, 0.76, 0.751, 0.739, 0.724, 0.73, 0.74, 0.737]) + +green = array( + [0.092, 0.096, 0.098, 0.097, 0.098, 0.095, 0.095, 0.097, 0.095, 0.094, 0.097, 0.098, 0.096, 0.101, 0.103, 0.104, + 0.107, 0.109, 0.112, 0.115, 0.125, 0.14, 0.16, 0.187, 0.229, 0.285, 0.343, 0.39, 0.435, 0.464, 0.472, 0.476, 0.481, + 0.462, 0.447, 0.441, 0.426, 0.406, 0.373, 0.347, 0.337, 0.314, 0.285, 0.277, 0.266, 0.25, 0.23, 0.207, 0.186, + 0.171, 0.16, 0.148, 0.141, 0.136, 0.13, 0.126, 0.123, 0.121, 0.122, 0.119, 0.114, 0.115, 0.117, 0.117, 0.118, 0.12, + 0.122, 0.128, 0.132, 0.139, 0.144, 0.146, 0.15, 0.152, 0.157, 0.159]) + +red = array( + [0.04, 0.046, 0.048, 0.053, 0.049, 0.05, 0.053, 0.055, 0.057, 0.056, 0.059, 0.057, 0.061, 0.061, 0.06, 0.062, 0.062, + 0.062, 0.061, 0.062, 0.06, 0.059, 0.057, 0.058, 0.058, 0.058, 0.056, 0.055, 0.056, 0.059, 0.057, 0.055, 0.059, + 0.059, 0.058, 0.059, 0.061, 0.061, 0.063, 0.063, 0.067, 0.068, 0.072, 0.08, 0.09, 0.099, 0.124, 0.154, 0.192, + 0.255, 0.287, 0.349, 0.402, 0.443, 0.487, 0.513, 0.558, 0.584, 0.62, 0.606, 0.609, 0.651, 0.612, 0.61, 0.65, 0.638, + 0.627, 0.62, 0.63, 0.628, 0.642, 0.639, 0.657, 0.639, 0.635, 0.642]) + +white_reflectivity = InterpolatedSF(wavelengths, white) +red_reflectivity = InterpolatedSF(wavelengths, red) +green_reflectivity = InterpolatedSF(wavelengths, green) + +# define light spectrum +light_spectrum = InterpolatedSF(array([400, 500, 600, 700]), array([0.0, 8.0, 15.6, 18.4])) + +# set-up scenegraph +world = World() + +# enclosing box +enclosure = Node(world) + +e_back = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, 0, 1) * rotate(0, 0, 0), + material=Lambert(white_reflectivity)) + +e_bottom = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, -1, 0) * rotate(0, -90, 0), + # material=m) + material=Lambert(white_reflectivity)) + +e_top = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(0, 1, 0) * rotate(0, 90, 0), + material=Lambert(white_reflectivity)) + +e_left = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(1, 0, 0) * rotate(-90, 0, 0), + material=Lambert(red_reflectivity)) + +e_right = Box(Point3D(-1, -1, 0), Point3D(1, 1, 0), + parent=enclosure, + transform=translate(-1, 0, 0) * rotate(90, 0, 0), + material=Lambert(green_reflectivity)) + +# ceiling light +light = Box(Point3D(-0.4, -0.4, -0.01), Point3D(0.4, 0.4, 0.0), + parent=enclosure, + transform=translate(0, 1, 0) * rotate(0, 90, 0), + material=UniformSurfaceEmitter(light_spectrum, 2)) + +# alternate light #1 +# light = Box(Point3D(-0.4, -0.4, -0.01), Point3D(0.4, 0.4, 0.0), +# parent=enclosure, +# transform=translate(0, 1, 0) * rotate(0, 90, 0), +# material=UniformSurfaceEmitter(d65_white, 2)) + +# alternate light #2 +# back_light = Sphere(0.1, +# parent=enclosure, +# transform=translate(0.80, -0.85, 0.80)*rotate(0, 0, 0), +# material=UniformSurfaceEmitter(light_spectrum, 10.0)) + +# objects in enclosure +box = Box(Point3D(-0.4, 0, -0.4), Point3D(0.3, 1.4, 0.3), + parent=world, + transform=translate(0.4, -1 + 1e-6, 0.4)*rotate(30, 0, 0), + material=schott("N-BK7")) + +sphere = Sphere(0.4, + parent=world, + transform=translate(-0.4, -0.6 + 1e-6, -0.4)*rotate(0, 0, 0), + material=schott("N-BK7")) + + +filter_red = InterpolatedSF([100, 650, 660, 670, 680, 800], [0, 0, 1, 1, 0, 0]) +filter_green = InterpolatedSF([100, 530, 540, 550, 560, 800], [0, 0, 1, 1, 0, 0]) +filter_blue = InterpolatedSF([100, 480, 490, 500, 510, 800], [0, 0, 1, 1, 0, 0]) + +# create and setup the camera +power_unfiltered = PowerPipeline2D(display_unsaturated_fraction=0.96, name="Unfiltered") +power_unfiltered.display_update_time = 15 + +power_green = PowerPipeline2D(filter=filter_green, display_unsaturated_fraction=0.96, name="Green Filter") +power_green.display_update_time = 15 + +power_red = PowerPipeline2D(filter=filter_red, display_unsaturated_fraction=0.96, name="Red Filter") +power_red.display_update_time = 15 + +rgb = RGBPipeline2D(display_unsaturated_fraction=0.96, name="sRGB") + +bayer = BayerPipeline2D(filter_red, filter_green, filter_blue, display_unsaturated_fraction=0.96, name="Bayer Filter") +bayer.display_update_time = 15 + +pipelines = [rgb, power_unfiltered, power_green, power_red, bayer] + +sampler = RGBAdaptiveSampler2D(rgb, ratio=10, fraction=0.2, min_samples=500, cutoff=0.01) + + +camera = PinholeCamera((1024, 1024), parent=world, transform=translate(0, 0, -3.3) * rotate(0, 0, 0), pipelines=pipelines) +camera.frame_sampler = sampler +camera.spectral_rays = 1 +camera.spectral_bins = 15 +camera.pixel_samples = 250 +camera.ray_importance_sampling = True +camera.ray_important_path_weight = 0.25 +camera.ray_max_depth = 500 +camera.ray_extinction_min_depth = 3 +camera.ray_extinction_prob = 0.01 + +# Speedups for testing purposes. +# Fewer camera pixels for faster runtime per pass. +camera.pixels = (256, 256) +# Less strict cutoff to finish with fewer passes. +sampler.cutoff = 0.05 + + +camera.render_engine = MPIEngine() + +# Don't make plots during the render as MPI is typically run non-interactively. +for pipeline in camera.pipelines: + pipeline.display_progress = False + +# Worker processes do not have accurate statistics, only the root process has +# all the render results. So don't bother outputting statistics on the workers. +if camera.render_engine.rank != 0: + camera.quiet = True + +# start ray tracing. +print(f"Starting ray tracing on rank {camera.render_engine.rank}") +render_pass = 1 +while not camera.render_complete: + if camera.render_engine.rank == 0: + print(f"Rendering pass {render_pass}...") + camera.observe() + # Rank 0 processes the ray tracing results so is the only process which + # knows the true progress of the render. In order to know when the render + # is complete, we need to ensure the frame sampler on each of the worker + # processes has an up-to-date copy of the pipeline used for sampling. + # We also need to use a different variable name for the received broadcast + # object else it is not broadcast properly. + root_rgb = camera.render_engine.comm.bcast(rgb, root=0) + camera.frame_sampler.pipeline = root_rgb + camera.render_engine.comm.Barrier() + render_pass += 1 +print(f"Finished ray tracing on rank {camera.render_engine.rank}") + +# Again, only rank 0 has all the ray tracing results so is the only one which +# can produce a correct image. +if camera.render_engine.rank == 0: + rgb.save('CornellMPI_rgb.png') diff --git a/demos/MPI/raysect_logo_hybrid.py b/demos/MPI/raysect_logo_hybrid.py new file mode 100644 index 00000000..a1327dde --- /dev/null +++ b/demos/MPI/raysect_logo_hybrid.py @@ -0,0 +1,90 @@ +""" +Renders the raysect logo: a top-down view of 6 coloured slabs. + +MPI (message passing interface) is used to parallelise the workflow. +Running this demo therefore requires an MPI library, such as MPICH or +OpenMPI. The mpi4py Python package must also be installed. + +To run this demo, call `mpirun -n python raysect_logo_mpi.py`. +If a batch job scheduler is used on a cluster, the render may be spread +across multiple nodes of the cluster. If run interactively at the command +line, all rendering will be done on the same machine as the command line. +""" + +from matplotlib.pyplot import * +from numpy import array + +from raysect.primitive import Sphere, Box + +from raysect.optical import World, Node, translate, rotate, Point3D, d65_white, ConstantSF, InterpolatedSF +from raysect.optical.observer import PinholeCamera +from raysect.optical.material.emitter import UniformSurfaceEmitter +from raysect.optical.material.dielectric import Dielectric +from raysect.core.workflow import MulticoreEngine, HybridEngine + + +world = World() + +wavelengths = array([300, 490, 510, 590, 610, 800]) +red_attn = array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0]) * 0.98 +green_attn = array([0.0, 0.0, 1.0, 1.0, 0.0, 0.0]) * 0.85 +blue_attn = array([1.0, 1.0, 0.0, 0.0, 0.0, 0.0]) * 0.98 +yellow_attn = array([0.0, 0.0, 1.0, 1.0, 1.0, 1.0]) * 0.85 +cyan_attn = array([1.0, 1.0, 1.0, 1.0, 0.0, 0.0]) * 0.85 +purple_attn = array([1.0, 1.0, 0.0, 0.0, 1.0, 1.0]) * 0.95 + +red_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, red_attn)) +green_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, green_attn)) +blue_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, blue_attn)) +yellow_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, yellow_attn)) +cyan_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, cyan_attn)) +purple_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, purple_attn)) + +Sphere(1000, world, material=UniformSurfaceEmitter(d65_white, 1.0)) + +node = Node(parent=world, transform=rotate(0, 0, 90)) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 0) * translate(0, 1, -0.500001), red_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 60) * translate(0, 1, -0.500001), yellow_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 120) * translate(0, 1, -0.500001), green_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 180) * translate(0, 1, -0.500001), cyan_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 240) * translate(0, 1, -0.500001), blue_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 300) * translate(0, 1, -0.500001), purple_glass) + +camera = PinholeCamera((256, 256), fov=45, parent=world, transform=translate(0, 0, -6.5) * rotate(0, 0, 0)) + +camera.ray_max_depth = 500 +camera.ray_extinction_prob = 0.01 +camera.pixel_samples = 100 +camera.spectral_rays = 1 +camera.spectral_bins = 21 + +# Work a bit harder across a cluster. +camera.pixels = (512, 512) +camera.pixel_samples = 500 + +# Either hard-code the number of sub-engine workers, or infer it from the run +# time environment if using a supported job scheduler like Slurm or Grid Engine. +# nworkers = 8 +nworkers = HybridEngine.estimate_subworker_count() +camera.render_engine = HybridEngine(MulticoreEngine(nworkers)) +print(f"MPI process {camera.render_engine.name} running with {nworkers} workers.") + +# With the hybrid engine, only the rank 0 process has all the results available +# to generate sampling statistics. So don't output stats on the workers. +if camera.render_engine.rank != 0: + camera.quiet = True + +# Don't try making plots during the render: this will likely fail when run +# through a non-interactive job scheduler. +for pipeline in camera.pipelines: + pipeline.display_progress = False + +camera.observe() + +# Again, only rank 0 has all the results, so this is the only process that +# should produce any output. +if camera.render_engine.rank == 0: + # Comment out if you don't want to save the result to a file: + camera.pipelines[0].save("raysect_logo.png") + # Uncomment if running interactively for a plot of the result: + # camera.pipelines[0].display() diff --git a/demos/MPI/raysect_logo_mpi.py b/demos/MPI/raysect_logo_mpi.py new file mode 100644 index 00000000..d70c745d --- /dev/null +++ b/demos/MPI/raysect_logo_mpi.py @@ -0,0 +1,81 @@ +""" +Renders the raysect logo: a top-down view of 6 coloured slabs. + +MPI (message passing interface) is used to parallelise the workflow. +Running this demo therefore requires an MPI library, such as MPICH, +OpenMPI or MicrosoftMPI. The mpi4py Python package must also +be installed. + +To run this demo, call `mpirun -n python raysect_logo_mpi.py`. +If a batch job scheduler is used on a cluster, the render may be spread +across multiple nodes of the cluster. If run interactively at the command +line, all rendering will be done on the same machine as the command line. +""" + +from matplotlib.pyplot import * +from numpy import array + +from raysect.primitive import Sphere, Box + +from raysect.optical import World, Node, translate, rotate, Point3D, d65_white, ConstantSF, InterpolatedSF +from raysect.optical.observer import PinholeCamera +from raysect.optical.material.emitter import UniformSurfaceEmitter +from raysect.optical.material.dielectric import Dielectric +from raysect.core.workflow import MPIEngine + + +world = World() + +wavelengths = array([300, 490, 510, 590, 610, 800]) +red_attn = array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0]) * 0.98 +green_attn = array([0.0, 0.0, 1.0, 1.0, 0.0, 0.0]) * 0.85 +blue_attn = array([1.0, 1.0, 0.0, 0.0, 0.0, 0.0]) * 0.98 +yellow_attn = array([0.0, 0.0, 1.0, 1.0, 1.0, 1.0]) * 0.85 +cyan_attn = array([1.0, 1.0, 1.0, 1.0, 0.0, 0.0]) * 0.85 +purple_attn = array([1.0, 1.0, 0.0, 0.0, 1.0, 1.0]) * 0.95 + +red_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, red_attn)) +green_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, green_attn)) +blue_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, blue_attn)) +yellow_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, yellow_attn)) +cyan_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, cyan_attn)) +purple_glass = Dielectric(index=ConstantSF(1.4), transmission=InterpolatedSF(wavelengths, purple_attn)) + +Sphere(1000, world, material=UniformSurfaceEmitter(d65_white, 1.0)) + +node = Node(parent=world, transform=rotate(0, 0, 90)) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 0) * translate(0, 1, -0.500001), red_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 60) * translate(0, 1, -0.500001), yellow_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 120) * translate(0, 1, -0.500001), green_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 180) * translate(0, 1, -0.500001), cyan_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 240) * translate(0, 1, -0.500001), blue_glass) +Box(Point3D(-0.5, 0, -2.5), Point3D(0.5, 0.25, 0.5), node, rotate(0, 0, 300) * translate(0, 1, -0.500001), purple_glass) + +camera = PinholeCamera((256, 256), fov=45, parent=world, transform=translate(0, 0, -6.5) * rotate(0, 0, 0)) + +camera.ray_max_depth = 500 +camera.ray_extinction_prob = 0.01 +camera.pixel_samples = 100 +camera.spectral_rays = 1 +camera.spectral_bins = 21 +camera.render_engine = MPIEngine() + +# With the MPI engine, only the rank 0 process has all the results available +# to generate sampling statistics. So don't output stats on the workers. +if camera.render_engine.rank != 0: + camera.quiet = True + +# Don't try making plots during the render: this will likely fail when run +# through a non-interactive job scheduler. +for pipeline in camera.pipelines: + pipeline.display_progress = False + +camera.observe() + +# Again, only rank 0 has all the results, so this is the only process that +# should produce any output. +if camera.render_engine.rank == 0: + # Comment out if you don't want to save the result to a file: + camera.pipelines[0].save("raysect_logo.png") + # Uncomment if running interactively for a plot of the result: + # camera.pipelines[0].display() diff --git a/docs/source/api_reference/core/render_engines.rst b/docs/source/api_reference/core/render_engines.rst index c0a8154a..c4ad8eef 100644 --- a/docs/source/api_reference/core/render_engines.rst +++ b/docs/source/api_reference/core/render_engines.rst @@ -11,5 +11,8 @@ Render Engines .. autoclass:: raysect.core.workflow.MulticoreEngine :show-inheritance: +.. autoclass:: raysect.core.workflow.MPIEngine + :show-inheritance: - +.. autoclass:: raysect.core.workflow.HybridEngine + :show-inheritance: diff --git a/raysect/core/meson.build b/raysect/core/meson.build index ab291f7a..a419221b 100644 --- a/raysect/core/meson.build +++ b/raysect/core/meson.build @@ -4,7 +4,7 @@ target_path = 'raysect/core' # source files -py_files = ['__init__.py', 'constants.py', 'workflow.py'] +py_files = ['__init__.py', 'constants.py'] pyx_files = ['boundingbox.pyx', 'boundingsphere.pyx', 'containers.pyx', 'intersection.pyx', 'material.pyx', 'ray.pyx'] pxd_files = ['__init__.pxd', 'boundingbox.pxd', 'boundingsphere.pxd', 'containers.pxd', 'intersection.pxd', 'material.pxd', 'ray.pxd'] data_files = [] @@ -31,3 +31,4 @@ subdir('acceleration') subdir('math') subdir('scenegraph') subdir('tests') +subdir('workflow') diff --git a/raysect/core/workflow/__init__.py b/raysect/core/workflow/__init__.py new file mode 100644 index 00000000..7e3b5aa6 --- /dev/null +++ b/raysect/core/workflow/__init__.py @@ -0,0 +1,33 @@ +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from .base import RenderEngine +from .serial import SerialEngine +from .multicore import MulticoreEngine +from .mpi import MPIEngine, HybridEngine diff --git a/raysect/core/workflow/base.py b/raysect/core/workflow/base.py new file mode 100644 index 00000000..784151ff --- /dev/null +++ b/raysect/core/workflow/base.py @@ -0,0 +1,93 @@ +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + + +class RenderEngine: + """ + Provides a common rendering workflow interface. + + This is a base class, its functionality must be implemented fully by the deriving class. + + This class provides a rendering workflow that abstracts away the underlying + system performing the work. It is intended that render engines may be built + that provide rendering on single cores, multi-cores (SMP) and clusters. + + The basic workflow is as follows. The render task is split into small, + self-contained chunks of work - 'tasks'. These tasks are passed to the + render engine which distributes the work to the available computing + resources. These discrete computing resources are know as "workers". + Workers process one task at a time and return their result to the render + engine. When results are received the render engine assembles them into + the final result. + + This workflow is implemented by supplying a set of tasks and two methods to + the render engines' run() method which processes those tasks. The functions + supplied to the run() method may be given additional args and kwargs. + + A worker calls render for each task object received. render has the following signature: :: + + def render(task, *render_args, **render_kwargs) + + where args and kwargs are additional arguments supplied by the user. + + Similarly, the worker calls update() for the results generated by a call to + render(). Update() has the following signature: :: + + def update(results, *update_args, **update_kwargs) + + where args and kwargs are additional arguments supplied by the user. + + The render() function must return an object representing the results, + this must be a picklable python object. + + The execution order of tasks is not guaranteed to be in order. If the order + is critical, an identifier should be passed as part of the task definition + and returned in the result. This will permit the order to be reconstructed. + """ + + def run(self, tasks, render, update, render_args=(), render_kwargs={}, update_args=(), update_kwargs={}): + """ + Starts the render engine executing the requested tasks. + + :param list tasks: List of user defined tuples that describe the task to execute. + :param object render: Callable python object that executes the tasks. + :param object update: Callable python object that is called following a render task and must be + used to update the internal state of the object requesting work. + :param tuple render_args: Additional arguments to pass to user defined render function. + :param tuple render_kwargs: Additional keyword arguments to pass to user defined render function. + :param tuple update_args: Additional arguments to pass to user defined update function. + :param tuple update_kwargs: Additional keyword arguments to pass to user defined update function. + """ + raise NotImplementedError("Virtual method must be implemented in sub-class.") + + def worker_count(self): + """ + Returns the number of workers in use by this engine. + """ + raise NotImplementedError("Virtual method must be implemented in sub-class.") diff --git a/raysect/core/workflow/meson.build b/raysect/core/workflow/meson.build new file mode 100644 index 00000000..ba0105fd --- /dev/null +++ b/raysect/core/workflow/meson.build @@ -0,0 +1,30 @@ +# WARNING: This file is automatically generated by dev/generate_meson_files.py. +# The template file used to generate this file is dev/subdir-meson.build. + +target_path = 'raysect/core/workflow' + +# source files +py_files = ['__init__.py', 'base.py', 'mpi.py', 'multicore.py', 'serial.py'] +pyx_files = [] +pxd_files = [] +data_files = [] + +# compile cython +foreach pyx_file: pyx_files + py.extension_module( + fs.replace_suffix(pyx_file, ''), + pyx_file, + dependencies: cython_dependencies, + install: true, + subdir: target_path, + cython_args: cython_args + ) +endforeach + +# add python, pxd and data files to the build +py.install_sources( + py_files + pxd_files + data_files, + subdir: target_path +) + +subdir('tests') diff --git a/raysect/core/workflow/mpi.py b/raysect/core/workflow/mpi.py new file mode 100644 index 00000000..582ce5f1 --- /dev/null +++ b/raysect/core/workflow/mpi.py @@ -0,0 +1,397 @@ +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from collections import defaultdict +from multiprocessing import get_context +import platform +import os + +from . import RenderEngine, SerialEngine + + +try: + from mpi4py import MPI + HAVE_MPI = True +except ImportError: + HAVE_MPI = False + + +class MPIEngine(RenderEngine): + """ + Render engine for running in an MPI context. + + This engine is useful for distributed memory systems, and for shared + memory systems where the overhead of inter-process communication of + the scenegraph is large compared with the time taken for a single + process to produce the scenegraph (e.g. on Windows). + + This render engine requires mpi4py to be installed, and the program + to be run using ``mpirun -n python ``. + + >>> from raysect.core import MPIEngine + >>> from raysect.optical.observer import PinholeCamera + >>> + >>> camera = PinholeCamera((512, 512)) + >>> camera.render_engine = MPIEngine() + + The render engine uses the single process, multiple data (SPMD) paradigm. + Each process is treated as a separate serial renderer. It is assumed that + the scene graph is created in every process, and each process processes a + subset of the total rendering tasks sequentially. The results are then + gathered back to the root (rank 0) process. + + The SPMD paradigm means there are many copies of the program running and + each copy runs the same instructions, so programs must be written with + this in mind. While all copies build the scene to render, only one copy + (the process with rank 0) actually receives all the render results + and calls the update function. This means that after an observe, + only the rank 0 process contains the results of the call to the render + function for all tasks, and so any post processing (including plotting + or saving images) should only be done on the rank 0 process. + + Also, any further renders which depend on the results of a previous render + (for example, using adapive samplers) will require communication of the + results from rank 0 to all the other processes: + + >>> pipeline = RGBPipeline2d() + >>> camera.pipelines = [pipeline] + >>> camera.sampler = RGBAdaptivesampler2d(pipeline) + >>> camera.observe() + >>> # Update the sampler in each worker process with the results + >>> # of the previous render for the next pass. + >>> root_pipeline = camera.render_engine.comm.bcast(pipeline, root=0) + >>> camera.sampler.pipeline = root_pipeline + >>> # Now subsequent observes will have the correct statistics. + >>> camera.observe() + + The class contains some attributes relevant to the MPI environment: + + :ivar comm: the MPI communicator (``MPI_COMM_WORLD``). + :ivar rank: the process rank in the communicator. Only rank 0 contains the + full results after a call ``to RenderEngine.run()``. + :ivar size: the number of processes in the communicator. + """ + def __init__(self): + if not HAVE_MPI: + raise RuntimeError("The mpi4py package is required to use this engine.") + comm = MPI.COMM_WORLD + self.comm = comm + self.rank = comm.Get_rank() + self.size = comm.Get_size() + if self.size < 2: + raise RuntimeError("At least 2 separate processes are required to use this engine.") + + def run(self, tasks, render, update, render_args=(), render_kwargs=None, update_args=(), update_kwargs=None): + # Avoid mutable default arguments. + if render_kwargs is None: + render_kwargs = {} + if update_kwargs is None: + update_kwargs = {} + # All processes must have the same tasks, in the same order, so that + # all processes agree on which subset of tasks each is to work on. + tasks = self.comm.bcast(tasks, root=0) + ntasks = len(tasks) + nworkers = self.size - 1 + worker_tasks = defaultdict(list) + for i, task in enumerate(tasks): + worker_tasks[i % nworkers].append(task) + if self.rank == 0: # The root node processes the results. + remaining = ntasks + while remaining: + result = self.comm.recv() + if isinstance(result, Exception): + raise result + update(result, *update_args, **update_kwargs) + remaining -= 1 + else: # Each worker renders the subset of tasks assigned to them. + # Unlike MulticoreEngine, there is no need to re-seed the random + # number generator to prevent all workers inheriting the same sequence + # since all processes are independent. + for task in worker_tasks[self.rank - 1]: + try: + result = render(task, *render_args, **render_kwargs) + except Exception as e: + result = e + self.comm.send(result, 0) + self.comm.Barrier() + + def worker_count(self): + return self.size + + +class HybridEngine(RenderEngine): + """ + Render engine for combined shared and distributed memory systems. + + This render engine requires mpi4py to be installed, and the program + to be run using ``mpirun -n python ``. + + >>> from raysect.core import HybridEngine, MulticoreEngine + >>> from raysect.optical.observer import PinholeCamera + >>> + >>> nworkers = 4 + >>> camera = PinholeCamera((512, 512)) + >>> camera.render_engine = HybridEngine(MulticoreEngine(nworkers)) + + The render engine uses a variation of the single process, multiple + data (SPMD) paradigm. In this paradigm, "tasks" are independent + processes potentially running on separate computers. Each task will + have one or more "workers" which it spawns in order to do its share + of the rendering. Each worker will be on the same computer as its + parent task and should share the parent task's memory. + + When the program is started with `mpirun`, there will be `ntasks` + separate processes, all of which run concurrently and build their + own copies of the scenegraph. Then when this render engine's `run` + method is called in each process, the sub-engine will create + `nworkers` worker subprograms to perform a subset of the render: how + this is done depends entirely on which sub-engine is used. For + example, if the `MulticoreEngine` is used then `nworkers` + subprocesses will be forked from the parent task to perform the + render subset. If the `SerialEngine` is used (the default) then the + render subset will be computed in serial in the task's own + process. The rendering results for all workers in all tasks are + gathered back to the root (rank 0) task. + + The SPMD paradigm means there are many copies of the program running + and each copy runs the same instructions, so programs must be + written with this in mind. While all copies build the scene to + render, only one copy (the MPI process with rank 0) actually + receives all the render results and calls the update function. This + means that after an observe, only the rank 0 process contains the + results of the call to the render function for all tasks, and so any + post processing (including plotting or saving images) should only be + done on the rank 0 process. + + Also, any further renders which depend on the results of a previous + render (for example, using adapive samplers) will require + communication of the results from rank 0 to all the other processes: + + >>> pipeline = RGBPipeline2d() + >>> camera.pipelines = [pipeline] + >>> camera.sampler = RGBAdaptivesampler2d(pipeline) + >>> camera.observe() + >>> # Update the sampler in each worker process with the results + >>> # of the previous render for the next pass. + >>> root_pipeline = camera.render_engine.comm.bcast(pipeline, root=0) + >>> camera.sampler.pipeline = root_pipeline + >>> # Now subsequent observes will have the correct statistics. + >>> camera.observe() + + It is the end user's responsibility to ensure that the number of + workers is configured appropriately, and this will strongly depend + on the environment in which the program is launched. Examples of + running the hybrid engine for 3 popular schedulers are given here. + There is also a helper method, ``estimate_subworker_count``, which + can assist for some schedulers. + + Slurm: + + $ sbatch -n -c + $ # Within script.sh: + $ mpirun -n --bind-to none + + >>> # Within application.py: + >>> camera.render_engine = HybridEngine(MulticoreEngine() + + PSB/Torque: + + $ qsub -l nodes=:ppn= + $ # Within script.sh: + $ mpirun -n --map-by node --bind-to none + + >>> # Within application.py: + >>> camera.render_engine = HybridEngine(MulticoreEngine() + + Grid engine is more complicated as there is no portable way to + specify the number of slots per node, though there is an environment + variable for the number of separate nodes the job is being run on: + + $ qsub -pe + $ # Within script.sh: + $ mpirun -n $NHOSTS --map-by node --bind-to none + + >>> # Within application.py: + >>> # Find out how many slots are allocated this MPI process's node. + >>> import os, platform + >>> host = platform.node() + >>> pe_hostfile = os.environ["PE_HOSTFILE"] + >>> with open(pe_hostfile, "r", encoding="UTF-8") as f: + >>> for line in f: + >>> hostname, slots, *_ = line.strip().split() + >>> if host in hostname: + >>> NSUB = int(slots) + >>> break + >>> + >>> camera.render_engine = HybridEngine(MulticoreEngine(NSUB)) + + This class contains some attributes relevant to the MPI environment. + The sub-engine attributes are accessible through the sub-engine + directly. + + :ivar comm: the MPI communicator (``MPI_COMM_WORLD``). + :ivar rank: the mpi rank in the communicator. Only rank 0 contains the + full results after a call ``to RenderEngine.run()``. + :ivar nmpi: the number of MPI processes. + :ivar name: the MPI processor name. + :ivar subengine: The render engine used by each individual task. + :ivar total_size: Total number of workers summed over all MPI processes. + """ + def __init__(self, subengine=None): + if not HAVE_MPI: + raise RuntimeError("The mpi4py package is required to use this engine.") + if subengine is None: + subengine = SerialEngine() + self.subengine = subengine + comm = MPI.COMM_WORLD + self.comm = comm + self.rank = comm.Get_rank() + self.nmpi = comm.Get_size() + self.name = MPI.Get_processor_name() + _local_size = subengine.worker_count() + self._all_sizes = comm.allgather(_local_size) + self.total_size = sum(self._all_sizes) + + def run(self, tasks, render, update, + render_args=(), render_kwargs=None, + update_args=(), update_kwargs=None): + # Avoid mutable default arguments. + if render_kwargs is None: + render_kwargs = {} + if update_kwargs is None: + update_kwargs = {} + # All processes must have the same tasks, in the same order, so that + # all processes agree on which subset of tasks each is to work on. + tasks = self.comm.bcast(tasks, root=0) + ntasks = len(tasks) + # Split tasks evenly per individual worker, as different MPI processes + # may have different numbers of workers. Then re-group the tasks into + # batches for each MPI worker. + nworkers = self.total_size + worker_tasks = [[] for _ in range(nworkers)] + for i, task in enumerate(tasks): + worker_tasks[i % nworkers].append(task) + rank_tasks = [] + idx = 0 + for rank in range(self.nmpi): + size = self._all_sizes[rank] + rank_tasks.append(sum(worker_tasks[idx:idx+size], [])) + idx += size + # On rank 0, spawn a separate process to do the rendering while gathering + # results in the main thread. We can't use mpi4py to send results from rank + # 0 to rank 0: here be segfaults. So use a multiprocessing queue instead. + # This does mean we have to gather data differently on rank 0. + # We use the 'fork' context here for memory efficiency, which makes this + # implementation unix-only (and even then, dodgy on MacOS). + if self.rank == 0: + ctx = get_context('fork') + queue = ctx.SimpleQueue() + worker = ctx.Process( + target=self.subengine.run, + kwargs=dict( + tasks=rank_tasks[self.rank], + render=render, render_args=render_args, render_kwargs=render_kwargs, + update=queue.put, update_args=(), update_kwargs={}, + ) + ) + worker.start() + remaining = ntasks + local_remaining = len(rank_tasks[self.rank]) + remote_remaining = remaining - local_remaining + while remaining: + # Data may come in either from the local queue or from remote MPI. + if local_remaining and not queue.empty(): + result = queue.get() + local_remaining -= 1 + elif remote_remaining: + result = self.comm.recv() + remote_remaining -= 1 + else: # Nothing local or remote yet. Try again. + continue + if isinstance(result, Exception): + raise result + update(result, *update_args, **update_kwargs) + remaining = local_remaining + remote_remaining + worker.join() + else: + # All other ranks render the subset of tasks assigned to them and + # send the results through MPI to rank 0 for collating. + self.subengine.run( + tasks=rank_tasks[self.rank], + render=render, render_args=render_args, render_kwargs=render_kwargs, + update=self.comm.send, update_args=(0,), update_kwargs={}, + ) + self.comm.Barrier() + + def worker_count(self): + return self.totalsize + + @staticmethod + def estimate_subworker_count(scheduler=None): + """ + Estimate the number of processes the sub-worker should use. + + This routine should be called in processes launched from + job schedulers. It will use environment variables exported + by the job scheduler to work out the degree of parallelism + available to the sub-worker. + + Currently, Slurm and Grid Engine are supported. For other + schedulers it is the user's responsibility to work out (or hard + code) the sub-worker count. + + :param scheduler: the name of the scheduler in use. + One of ("slurm", "gridengine") + :return: The number of sub-workers available in this MPI task. + """ + if scheduler is None: + # Try to work out which scheduler is in use. + if "SLURM_JOB_ID" in os.environ: + scheduler = "slurm" + elif "SGE_ROOT" in os.environ: + scheduler = "gridengine" + else: + raise RuntimeError("Can't find a supported scheduler.") + scheduler = scheduler.lower() + if scheduler == "slurm": + nsub = int(os.getenv("SLURM_CPUS_PER_TASK", "1")) + elif scheduler == "gridengine": + # Parse the parallel environment file to get the number of slots + # allocated to this node. + node = platform.node() + pe_hostfile = os.environ["PE_HOSTFILE"] + with open(pe_hostfile, "r", encoding="utf-8") as f: + for line in f: + hostname, slots, *_ = line.strip().split() + if node in hostname: + nsub = int(slots) + break + else: + raise RuntimeError(f"{scheduler} is not supported by this function.") + return nsub diff --git a/raysect/core/workflow.py b/raysect/core/workflow/multicore.py similarity index 65% rename from raysect/core/workflow.py rename to raysect/core/workflow/multicore.py index 859d9123..a803a078 100644 --- a/raysect/core/workflow.py +++ b/raysect/core/workflow/multicore.py @@ -1,4 +1,4 @@ -# Copyright (c) 2014-2025, Dr Alex Meakins, Raysect Project +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project # All rights reserved. # # Redistribution and use in source and binary forms, with or without @@ -28,96 +28,9 @@ # POSSIBILITY OF SUCH DAMAGE. from multiprocessing import get_context, cpu_count -from raysect.core.math import random +import random import time - - -class RenderEngine: - """ - Provides a common rendering workflow interface. - - This is a base class, its functionality must be implemented fully by the deriving class. - - This class provides a rendering workflow that abstracts away the underlying - system performing the work. It is intended that render engines may be built - that provide rendering on single cores, multi-cores (SMP) and clusters. - - The basic workflow is as follows. The render task is split into small, - self-contained chunks of work - 'tasks'. These tasks are passed to the - render engine which distributes the work to the available computing - resources. These discrete computing resources are know as "workers". - Workers process one task at a time and return their result to the render - engine. When results are received the render engine assembles them into - the final result. - - This workflow is implemented by supplying a set of tasks and two methods to - the render engines' run() method which processes those tasks. The functions - supplied to the run() method may be given additional args and kwargs. - - A worker calls render for each task object received. render has the following signature: :: - - def render(task, *render_args, **render_kwargs) - - where args and kwargs are additional arguments supplied by the user. - - Similarly, the worker calls update() for the results generated by a call to - render(). Update() has the following signature: :: - - def update(results, *update_args, **update_kwargs) - - where args and kwargs are additional arguments supplied by the user. - - The render() function must return an object representing the results, - this must be a picklable python object. - - The execution order of tasks is not guaranteed to be in order. If the order - is critical, an identifier should be passed as part of the task definition - and returned in the result. This will permit the order to be reconstructed. - """ - - def run(self, tasks, render, update, render_args=(), render_kwargs={}, update_args=(), update_kwargs={}): - """ - Starts the render engine executing the requested tasks. - - :param list tasks: List of user defined tuples that describe the task to execute. - :param object render: Callable python object that executes the tasks. - :param object update: Callable python object that is called following a render task and must be - used to update the internal state of the object requesting work. - :param tuple render_args: Additional arguments to pass to user defined render function. - :param tuple render_kwargs: Additional keyword arguments to pass to user defined render function. - :param tuple update_args: Additional arguments to pass to user defined update function. - :param tuple update_kwargs: Additional keyword arguments to pass to user defined update function. - """ - raise NotImplementedError("Virtual method must be implemented in sub-class.") - - def worker_count(self): - """ - Returns the number of workers in use by this engine. - """ - raise NotImplementedError("Virtual method must be implemented in sub-class.") - - -class SerialEngine(RenderEngine): - """ - Render engine for running on a single CPU processor. - - This engine is useful for debugging. - - >>> from raysect.core import SerialEngine - >>> from raysect.optical.observer import PinholeCamera - >>> - >>> camera = PinholeCamera((512, 512)) - >>> camera.render_engine = SerialEngine() - """ - - def run(self, tasks, render, update, render_args=(), render_kwargs={}, update_args=(), update_kwargs={}): - - for task in tasks: - result = render(task, *render_args, **render_kwargs) - update(result, *update_args, **update_kwargs) - - def worker_count(self): - return 1 +from .base import RenderEngine class MulticoreEngine(RenderEngine): @@ -126,7 +39,7 @@ class MulticoreEngine(RenderEngine): The number of processes spawned by this render engine is controlled via the processes attribute. This can also be set at object initialisation. - + If the processes attribute is set to None (the default), the render engine will automatically set the number of processes to be equal to the number of CPU cores detected on the machine. @@ -324,36 +237,3 @@ def _worker(self, render, args, kwargs, job_queue, result_queue): # hand back results result_queue.put(results) - - -if __name__ == '__main__': - - class Job: - - def __init__(self, engine=None): - self.total = 0 - self.engine = engine if engine else MulticoreEngine() - - def run(self, v): - self.total = 0 - self.engine.run(list(range(v)), self.render, self.update, render_args=(10000,)) - return self.total - - def render(self, task, count): - sum = 0 - for i in range(count): - sum += 1 / count - return sum - - def update(self, result): - self.total += result - - n = 20000 - - t = time.time() - j = Job(SerialEngine()) - print(j.run(n), time.time() - t) - - t = time.time() - j = Job(MulticoreEngine()) - print(j.run(n), time.time() - t) diff --git a/raysect/core/workflow/serial.py b/raysect/core/workflow/serial.py new file mode 100644 index 00000000..90cdf2e8 --- /dev/null +++ b/raysect/core/workflow/serial.py @@ -0,0 +1,53 @@ +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from .base import RenderEngine + + +class SerialEngine(RenderEngine): + """ + Render engine for running on a single CPU processor. + + This engine is useful for debugging. + + >>> from raysect.core import SerialEngine + >>> from raysect.optical.observer import PinholeCamera + >>> + >>> camera = PinholeCamera((512, 512)) + >>> camera.render_engine = SerialEngine() + """ + + def run(self, tasks, render, update, render_args=(), render_kwargs={}, update_args=(), update_kwargs={}): + + for task in tasks: + result = render(task, *render_args, **render_kwargs) + update(result, *update_args, **update_kwargs) + + def worker_count(self): + return 1 diff --git a/raysect/core/workflow/tests/__init__.py b/raysect/core/workflow/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/raysect/core/workflow/tests/meson.build b/raysect/core/workflow/tests/meson.build new file mode 100644 index 00000000..a260fbe6 --- /dev/null +++ b/raysect/core/workflow/tests/meson.build @@ -0,0 +1,29 @@ +# WARNING: This file is automatically generated by dev/generate_meson_files.py. +# The template file used to generate this file is dev/subdir-meson.build. + +target_path = 'raysect/core/workflow/tests' + +# source files +py_files = ['__init__.py', 'test_workflow.py'] +pyx_files = [] +pxd_files = [] +data_files = [] + +# compile cython +foreach pyx_file: pyx_files + py.extension_module( + fs.replace_suffix(pyx_file, ''), + pyx_file, + dependencies: cython_dependencies, + install: true, + subdir: target_path, + cython_args: cython_args + ) +endforeach + +# add python, pxd and data files to the build +py.install_sources( + py_files + pxd_files + data_files, + subdir: target_path +) + diff --git a/raysect/core/workflow/tests/readme.txt b/raysect/core/workflow/tests/readme.txt new file mode 100644 index 00000000..2bd1fb2f --- /dev/null +++ b/raysect/core/workflow/tests/readme.txt @@ -0,0 +1,14 @@ +This directory contains tests for the render engines. +The Serial and MultiCore render engine tests can be run with the standard test command `./dev/test.sh`. +However, the MPI-based render engines (MPIEngine and HybridEngine) require an MPI runtime for proper testing. +They are therefore skipped if this is not available. + +To run the MPI render engine tests, it's necessary to run the test suite using MPI. This has 2 prerequisites: + +* Install the mpi4py Python package and an MPI runtime (impi-rt is a suitable manylinux runtime for Intel processors). +* Run the test suite using `mpirun -np 2 ./dev/test.sh -k mpi -k hybrid`. + +It's recommended to use the `-k mpi -k hybrid` filters to only run the tests requiring MPI: all other tests should be run separately. +Note that when using `mpirun -np 2` there will be 2 outputs for every test. +It's therefore best to limit to 2 processes to minimise the duplicated output: 2 is the minimum number required to test the MPI functionality properly. +By limiting to only the MPI tests the amount of duplicated output will be further reduced to only these tests. diff --git a/raysect/core/workflow/tests/test_workflow.py b/raysect/core/workflow/tests/test_workflow.py new file mode 100644 index 00000000..cbfd4cf9 --- /dev/null +++ b/raysect/core/workflow/tests/test_workflow.py @@ -0,0 +1,140 @@ +# Copyright (c) 2014-2026, Dr Alex Meakins, Raysect Project +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Raysect Project nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import unittest +from raysect.core.workflow import SerialEngine, MulticoreEngine, MPIEngine, HybridEngine +from raysect.core.workflow.mpi import HAVE_MPI + + +if HAVE_MPI: + from mpi4py import MPI + MPI_WORLD_SIZE = MPI.COMM_WORLD.size +else: + MPI_WORLD_SIZE = 0 + + +class TestWorkflow(unittest.TestCase): + """Test the render engines.""" + + def setUp(self): + self.total = 0 + # A list of numbers [1, N] to perform a weighted sum of. + self.N = 100 + self.numbers = list(range(1, self.N + 1)) + # Each task represents an index into the list of numbers. + self.tasks = [i for i, _ in enumerate(self.numbers)] + + def render(self, task, weight=1): + index = task + number = self.numbers[index] + scaled = number * weight + return scaled + + def update(self, result, neg=False): + if neg: + result = -result + self.total += result + + def test_serial(self): + engine = SerialEngine() + engine.run(self.tasks, self.render, self.update) + self.assertEqual(self.total, self.N * (self.N + 1) / 2) + + def test_serial_args(self): + engine = SerialEngine() + engine.run(self.tasks, self.render, self.update, render_args=(2,), update_args=(True,)) + self.assertEqual(self.total, -self.N * (self.N + 1)) + + def test_serial_kwargs(self): + engine = SerialEngine() + engine.run(self.tasks, self.render, self.update, render_kwargs={'weight': 2}, update_args={'neg': True}) + self.assertEqual(self.total, -self.N * (self.N + 1)) + + def test_multicore(self): + engine = MulticoreEngine(2) + engine.run(self.tasks, self.render, self.update) + self.assertEqual(self.total, self.N * (self.N + 1) / 2) + + def test_multicore_args(self): + engine = MulticoreEngine() + engine.run(self.tasks, self.render, self.update, render_args=(2,), update_args=(True,)) + self.assertEqual(self.total, -self.N * (self.N + 1)) + + def test_multicore_kwargs(self): + engine = MulticoreEngine() + engine.run(self.tasks, self.render, self.update, render_kwargs={'weight': 2}, update_args={'neg': True}) + self.assertEqual(self.total, -self.N * (self.N + 1)) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_mpi(self): + engine = MPIEngine() + engine.run(self.tasks, self.render, self.update) + if engine.rank == 0: + self.assertEqual(self.total, self.N * (self.N + 1) / 2) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_mpi_args(self): + engine = MPIEngine() + engine.run(self.tasks, self.render, self.update, render_args=(2,), update_args=(True,)) + if engine.rank == 0: + self.assertEqual(self.total, -self.N * (self.N + 1)) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE == 1, "Requires MPI and running with mpirun -np 1.") + def test_mpi_disallow_1proc_world(self): + with self.assertRaises(RuntimeError): + MPIEngine() + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_hybrid_serial(self): + engine = HybridEngine(SerialEngine()) + engine.run(self.tasks, self.render, self.update) + if engine.rank == 0: + self.assertEqual(self.total, self.N * (self.N + 1) / 2) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_hybrid_serial_args(self): + engine = HybridEngine(SerialEngine()) + engine.run(self.tasks, self.render, self.update, render_args=(2,), update_args=(True,)) + if engine.rank == 0: + self.assertEqual(self.total, -self.N * (self.N + 1)) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_hybrid_multicore(self): + engine = HybridEngine(MulticoreEngine(2)) + engine.run(self.tasks, self.render, self.update) + if engine.rank == 0: + self.assertEqual(self.total, self.N * (self.N + 1) / 2) + + @unittest.skipUnless(HAVE_MPI and MPI_WORLD_SIZE > 1, "Requires MPI and running with mpirun -np >=2.") + def test_hybrid_multicore_args(self): + engine = HybridEngine(MulticoreEngine(2)) + engine.run(self.tasks, self.render, self.update, render_args=(2,), update_args=(True,)) + if engine.rank == 0: + self.assertEqual(self.total, -self.N * (self.N + 1)) +