diff --git a/README.md b/README.md index 3a386727..7b8c6103 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,8 @@ Monte Carlo eXtreme (MCX) - CUDA Edition - Author: Qianqian Fang (q.fang at neu.edu) - License: GNU General Public License version 3 (GPLv3) - Version: 2.2.pre (v2024.1, Interstellar Ion) -- Website: +- Website: +- Download: ![Mex and Binaries](https://github.com/fangq/mcx/actions/workflows/build_all.yml/badge.svg)\ ![Linux Python Module](https://github.com/fangq/mcx/actions/workflows/build_linux_manywheel.yml/badge.svg)\ diff --git a/README.txt b/README.txt index ab8fbf7a..7e48af46 100644 --- a/README.txt +++ b/README.txt @@ -6,8 +6,8 @@ *Author: Qianqian Fang *License: GNU General Public License version 3 (GPLv3) *Version: 2.2.pre (v2024.1, Interstellar Ion) -*Website: http://mcx.space - +*Website: https://mcx.space +*Download: https://mcx.space/wiki/?Get --------------------------------------------------------------------- Table of Content: diff --git a/example/multisrc/run_eachsrc.bat b/example/multisrc/run_eachsrc.bat new file mode 100644 index 00000000..6478168b --- /dev/null +++ b/example/multisrc/run_eachsrc.bat @@ -0,0 +1 @@ +..\..\bin\mcx.exe -f multisrc.json -D P -s 'eachsrc' -j '{"Optode":{"Source":{"ID":-1}}}' %* diff --git a/example/multisrc/run_eachsrc.sh b/example/multisrc/run_eachsrc.sh new file mode 100755 index 00000000..4f2ce022 --- /dev/null +++ b/example/multisrc/run_eachsrc.sh @@ -0,0 +1,3 @@ + +#!/bin/sh +../../bin/mcx -f multisrc.json -D P -s 'eachsrc' --srcid -1 $@ diff --git a/example/multisrc/run_src9.bat b/example/multisrc/run_src9.bat new file mode 100644 index 00000000..79605b83 --- /dev/null +++ b/example/multisrc/run_src9.bat @@ -0,0 +1 @@ +..\..\bin\mcx.exe -f multisrc.json -D P -s 'source9' --srcid 9 %* diff --git a/example/multisrc/run_src9.sh b/example/multisrc/run_src9.sh new file mode 100755 index 00000000..c1c752eb --- /dev/null +++ b/example/multisrc/run_src9.sh @@ -0,0 +1,3 @@ + +#!/bin/sh +../../bin/mcx -f multisrc.json -D P -s 'source9' --srcid 9 $@ diff --git a/mcxlab/examples/demo_multisrc.m b/mcxlab/examples/demo_multisrc.m index a0b9548d..6b4f4614 100644 --- a/mcxlab/examples/demo_multisrc.m +++ b/mcxlab/examples/demo_multisrc.m @@ -80,3 +80,21 @@ flux=mcxlab(cfg); mcxplotvol(log10(flux.data)); + +%% setting cfg.srcid to a positive number 1,2,3,.. specifies which src to simulate +cfg.srcid=8; + +flux=mcxlab(cfg); + +mcxplotvol(log10(flux.data)); + +%% setting cfg.srcid to -1, solution of each src is stored separately + +cfg.srcid=-1; + +flux=mcxlab(cfg); + +% use up-down button to shift between different sources, there are 9 of +% them + +mcxplotvol(log10(squeeze(flux.data))); diff --git a/mcxlab/mcxlab.m b/mcxlab/mcxlab.m index 62c61aac..cf62e76c 100644 --- a/mcxlab/mcxlab.m +++ b/mcxlab/mcxlab.m @@ -246,6 +246,9 @@ % simultaneously simulated; only works for 'pattern' % source, see cfg.srctype='pattern' for details % Example +% cfg.srcid: when multiple sources are defined, if srcid is -1, each source is separately +% simulated; if set to 0, all source solution are summed; if set to a positive +% number starting from 1, only the specified source is simulated; default to 0 % cfg.omega: source modulation frequency (rad/s) for RF replay, 2*pi*f % cfg.srciquv: 1x4 vector [I,Q,U,V], Stokes vector of the incident light % I: total light intensity (I >= 0) diff --git a/pmcx/README.md b/pmcx/README.md index add70604..5dd0c814 100644 --- a/pmcx/README.md +++ b/pmcx/README.md @@ -4,7 +4,7 @@ - Copyright: (C) Matin Raayai Ardakani (2022-2023) , Qianqian Fang (2019-2023) , Fan-Yu Yen (2023) - License: GNU Public License V3 or later -- Version: 0.2.9 +- Version: 0.2.10 - URL: https://pypi.org/project/pmcx/ - Github: https://github.com/fangq/mcx diff --git a/pmcx/pmcx/__init__.py b/pmcx/pmcx/__init__.py index 44c1aa38..21c3844d 100644 --- a/pmcx/pmcx/__init__.py +++ b/pmcx/pmcx/__init__.py @@ -49,7 +49,7 @@ # from .files import loadmc2, loadmch, load, save from .bench import bench -__version__ = "0.2.9" +__version__ = "0.2.10" __all__ = ( "gpuinfo", diff --git a/pmcx/setup.py b/pmcx/setup.py index df050c39..fd04c75c 100644 --- a/pmcx/setup.py +++ b/pmcx/setup.py @@ -123,7 +123,7 @@ def build_extension(self, ext): setup( name="pmcx", packages=["pmcx"], - version="0.2.9", + version="0.2.10", requires=["numpy"], license="GPLv3+", author="Matin Raayai Ardakani, Qianqian Fang, Fan-Yu Yen", diff --git a/src/Makefile b/src/Makefile index f7f86d60..c4feaceb 100644 --- a/src/Makefile +++ b/src/Makefile @@ -74,7 +74,7 @@ NVCCOMP=$(OMP) CUDA_STATIC=--cudart static -Xcompiler "-static-libgcc -static-libstdc++" CFLAGS+=-std=c99 -CPPFLAGS+=-g -Wall -pedantic #-DNO_LZMA # -DUSE_OS_TIMER +CPPFLAGS+=-g -Wall #-pedantic #-DNO_LZMA # -DUSE_OS_TIMER OBJSUFFIX=.o EXESUFFIX= diff --git a/src/mcx_core.cu b/src/mcx_core.cu index fb38f06b..4abbb5e9 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -375,6 +375,10 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float* if (SAVE_DETID(gcfg->savedetflag)) { n_det[baseaddr++] = detid; + + if (gcfg->extrasrclen && gcfg->srcid <= 0) { + n_det[baseaddr++] = (((int)ppath[gcfg->w0offset - 1]) << 16); + } } for (i = 0; i < gcfg->partialdata; i++) { @@ -392,7 +396,7 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float* } if (SAVE_W0(gcfg->savedetflag)) { - n_det[baseaddr++] = ppath[gcfg->w0offset - 1]; + n_det[baseaddr++] = ppath[gcfg->w0offset - 2]; } if (SAVE_IQUV(gcfg->savedetflag)) { @@ -415,7 +419,7 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float* * @param[in] gdebugdata: pointer to the global-memory buffer to store the trajectory info */ -__device__ inline void savedebugdata(MCXpos* p, uint id, float* gdebugdata) { +__device__ inline void savedebugdata(MCXpos* p, uint id, float* gdebugdata, int srcid) { uint pos = atomicAdd(gjumpdebug, 1); if (pos < gcfg->maxjumpdebug) { @@ -425,7 +429,7 @@ __device__ inline void savedebugdata(MCXpos* p, uint id, float* gdebugdata) { gdebugdata[pos++] = p->y; gdebugdata[pos++] = p->z; gdebugdata[pos++] = p->w; - gdebugdata[pos++] = 0; + gdebugdata[pos++] = srcid; } } @@ -1039,12 +1043,17 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* *w0 = 1.f; //< reuse to count for launchattempt int canfocus = 1; //< non-zero: focusable, zero: not focusable MCXSrc* launchsrc = &(gcfg->src); + ppath[gcfg->w0offset - 1] = 0.f; - if (gcfg->extrasrclen) { - int srcid = rand_uniform01(t) * (gcfg->extrasrclen + 1); + if (gcfg->extrasrclen && gcfg->srcid != 1) { + if (gcfg->srcid > 1) { + launchsrc = (MCXSrc*)(gproperty + gcfg->maxmedia + 1 + gcfg->detnum + ((gcfg->srcid - 2) * 4)); + } else { // gcfg->srcid = 0 or -1: simulate all sources; = 0 merge all solutions; = -1 separately store each source + ppath[gcfg->w0offset - 1] = (int)(rand_uniform01(t) * JUST_BELOW_ONE * (gcfg->extrasrclen + 1)) + 1; // borrow initial weight section of photon-sharing for storing launch src id - if (srcid) { - launchsrc = (MCXSrc*)(gproperty + gcfg->maxmedia + 1 + gcfg->detnum + ((srcid - 1) << 2)); + if ((int)ppath[gcfg->w0offset - 1] > 1) { + launchsrc = (MCXSrc*)(gproperty + gcfg->maxmedia + 1 + gcfg->detnum + ((int)(ppath[gcfg->w0offset - 1] - 2) * 4)); + } } } @@ -1065,13 +1074,17 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* ppath[gcfg->partialdata] += p->w; //< sum all the remaining energy if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) { - savedebugdata(p, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata); + savedebugdata(p, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]); } if (*mediaid == 0 && *idx1d != OUTSIDE_VOLUME_MIN && *idx1d != OUTSIDE_VOLUME_MAX && gcfg->issaveref) { if (gcfg->issaveref == 1) { int tshift = MIN(gcfg->maxgate - 1, (int)(floorf((f->t - gcfg->twin0) * gcfg->Rtstep))); + if (gcfg->extrasrclen && gcfg->srcid < 0) { + tshift += ((int)ppath[gcfg->w0offset - 1] - 1) * gcfg->maxgate; + } + if (gcfg->srctype != MCX_SRC_PATTERN && gcfg->srctype != MCX_SRC_PATTERN3D) { #ifdef USE_ATOMIC #ifdef USE_DOUBLE @@ -1206,12 +1219,12 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* if (gcfg->srctype == MCX_SRC_PATTERN) { // need to prevent rx/ry=1 here if (gcfg->srcnum <= 1) { p->w = launchsrc->pos.w * srcpattern[(int)(ry * JUST_BELOW_ONE * launchsrc->param2.w) * (int)(launchsrc->param1.w) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.w)]; - ppath[3] = p->w; + ppath[4] = p->w; } else { *((uint*)(ppath + 2)) = ((int)(ry * JUST_BELOW_ONE * launchsrc->param2.w) * (int)(launchsrc->param1.w) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.w)); for (int i = 0; i < gcfg->srcnum; i++) { - ppath[i + 3] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; + ppath[i + 4] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; } p->w = 1.f; @@ -1220,13 +1233,13 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* if (gcfg->srcnum <= 1) { p->w = launchsrc->pos.w * srcpattern[(int)(rz * JUST_BELOW_ONE * launchsrc->param1.z) * (int)(launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(ry * JUST_BELOW_ONE * launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.x)]; - ppath[3] = p->w; + ppath[4] = p->w; } else { *((uint*)(ppath + 2)) = ((int)(rz * JUST_BELOW_ONE * launchsrc->param1.z) * (int)(launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(ry * JUST_BELOW_ONE * launchsrc->param1.y) * (int)(launchsrc->param1.x) + (int)(rx * JUST_BELOW_ONE * launchsrc->param1.x)); for (int i = 0; i < gcfg->srcnum; i++) { - ppath[i + 3] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; + ppath[i + 4] = srcpattern[(*((uint*)(ppath + 2))) * gcfg->srcnum + i]; } p->w = 1.f; @@ -1558,7 +1571,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* updateproperty(prop, *mediaid, t, *idx1d, media, (float3*)p, nuvox, flipdir); if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) { - savedebugdata(p, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata); + savedebugdata(p, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[3]); } /** @@ -1574,7 +1587,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* if (gcfg->outputtype == otRF) { // if run RF replay f->pathlen = photontof[(threadid * gcfg->threadphoton + min(threadid, gcfg->oddphotons - 1) + (int)f->ndone)]; - sincosf(gcfg->omega * f->pathlen, ppath + 4 + gcfg->srcnum, ppath + 3 + gcfg->srcnum); + sincosf(gcfg->omega * f->pathlen, ppath + 5 + gcfg->srcnum, ppath + 4 + gcfg->srcnum); } f->pathlen = 0.f; @@ -1896,6 +1909,11 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], tmp0 = (gcfg->outputtype == otDCS) ? (1.f - ctheta) : 1.f; tshift = (int)(floorf((photontof[tshift] - gcfg->twin0) * gcfg->Rtstep)) + ( (gcfg->replaydet == -1) ? ((photondetid[tshift] - 1) * gcfg->maxgate) : 0); + + if (gcfg->extrasrclen && gcfg->srcid < 0) { + tshift += ((int)ppath[gcfg->w0offset - 1] - 1) * gcfg->maxgate; + } + #ifdef USE_ATOMIC if (!gcfg->isatomic) { @@ -1924,7 +1942,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], } if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) { - savedebugdata(&p, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata); + savedebugdata(&p, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]); } } @@ -2074,6 +2092,10 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], weight = w0 * f.pathlen; } + if (gcfg->extrasrclen && gcfg->srcid < 0) { + tshift += ((int)ppath[gcfg->w0offset - 1] - 1) * gcfg->maxgate; + } + GPUDEBUG(("deposit to [%d] %e, w=%f\n", idx1dold, weight, p.w)); if (fabsf(weight) > 0.f || gcfg->outputtype == otRF) { @@ -2607,7 +2629,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { uint sharedbuf = 0; /** \c dimxyz - output volume variable \c field voxel count, Nx*Ny*Nz*Ns where Ns=cfg.srcnum is the pattern number for photon sharing */ - int dimxyz = cfg->dim.x * cfg->dim.y * cfg->dim.z * ((cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) ? cfg->srcnum : 1); + int dimxyz = cfg->dim.x * cfg->dim.y * cfg->dim.z * ((cfg->srctype == MCX_SRC_PATTERN || cfg->srctype == MCX_SRC_PATTERN3D) ? cfg->srcnum : (cfg->srcid == -1) ? (cfg->extrasrclen + 1) : 1); /** \c media - input volume representing the simulation domain, format specified in cfg.mediaformat, read-only */ uint* media = (uint*)(cfg->vol); @@ -2661,15 +2683,15 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { * |-----------------------------------------------> hostdetreclen <--------------------------------------| * |------------------------> partialdata <-------------------| *host detected photon buffer: detid (1), partial_scat (#media), partial_path (#media), momontum (#media), p_exit (3), v_exit(3), w0 (1) - * |---------------------------------------------> w0offset <-------------------------------------||<----- w0 (#srcnum) ----->||<- RF replay (2)->| - *gpu detected photon buffer: partial_scat (#media), partial_path (#media), momontum (#media), E_escape (1), E_launch (1), w0 (1), w0_photonsharing (#srcnum) cos(w*T),sin(w*T) + * |---------------------------------------------> w0offset <-----------------------------------------------||<----- w0 (#srcnum) ----->||<- RF replay (2)->| + *gpu detected photon buffer: partial_scat (#media), partial_path (#media), momontum (#media), E_escape (1), E_launch (1), w0 (1), srcid(1), w0_photonsharing (#srcnum) cos(w*T),sin(w*T) */ //< \c partialdata: per-photon buffer length for media-specific data, copy from GPU to host unsigned int partialdata = (cfg->medianum - 1) * (SAVE_NSCAT(cfg->savedetflag) + SAVE_PPATH(cfg->savedetflag) + SAVE_MOM(cfg->savedetflag)); //< \c w0offset - offset in the per-photon buffer to the start of the photon sharing related data - unsigned int w0offset = partialdata + 3; + unsigned int w0offset = partialdata + 4; //< the extra 4 numbers are total-escaped-energy, total-launched-energy, initial-weight, source_id //< \c hostdetreclen - host-side det photon data buffer per-photon length unsigned int hostdetreclen = partialdata + SAVE_DETID(cfg->savedetflag) + 3 * (SAVE_PEXIT(cfg->savedetflag) + SAVE_VEXIT(cfg->savedetflag)) + SAVE_W0(cfg->savedetflag) + 4 * SAVE_IQUV(cfg->savedetflag); @@ -2680,8 +2702,8 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { /** \c param - constants to be used in the GPU, copied to GPU as \c gcfg, stored in the constant memory */ MCXParam param = {cfg->steps, minstep, 0, 0, cfg->tend, R_C0* cfg->unitinmm, (uint)cfg->issave2pt, (uint)cfg->isreflect, (uint)cfg->isrefint, (uint)cfg->issavedet, 1.f / cfg->tstep, - p0, c0, cfg->srcparam1, cfg->srcparam2, cfg->extrasrclen, s0, maxidx, uint4(0, 0, 0, 0), cp0, cp1, uint2(0, 0), cfg->minenergy, - cfg->sradius* cfg->sradius, minstep* R_C0* cfg->unitinmm, cfg->srctype, + p0, c0, cfg->srcparam1, cfg->srcparam2, cfg->extrasrclen, cfg->srcid, s0, maxidx, uint4(0, 0, 0, 0), + cp0, cp1, uint2(0, 0), cfg->minenergy, cfg->sradius* cfg->sradius, minstep* R_C0* cfg->unitinmm, cfg->srctype, cfg->voidtime, cfg->maxdetphoton, cfg->medianum - 1, cfg->detnum, cfg->polmedianum, cfg->maxgate, ABS(cfg->sradius + 2.f) < EPS /*isatomic*/, (uint)cfg->maxvoidstep, cfg->issaveseed > 0, (uint)cfg->issaveref, cfg->isspecular > 0, @@ -2782,6 +2804,8 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { param.maxgate = gpu[gpuid].maxgate; + printf("allocate %d floats [%d %d %d]\n", dimxyz * gpu[gpuid].maxgate * 2, dimxyz, gpu[gpuid].maxgate, cfg->detnum); + /** If cfg.respin is positive, the output data have to be accummulated, so we use a double-buffer to retrieve and then accummulate */ if (ABS(cfg->respin) > 1) { if (cfg->seed == SEED_FROM_FILE && cfg->replaydet == -1) { @@ -3051,9 +3075,9 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { /** Inside the GPU kernel, volume is always assumbed to be col-major (like those generated by MATLAB or FORTRAN) */ cachebox.x = (cp1.x - cp0.x + 1); cachebox.y = (cp1.y - cp0.y + 1) * (cp1.x - cp0.x + 1); + dimlen.x = cfg->dim.x; dimlen.y = cfg->dim.y * cfg->dim.x; - dimlen.z = cfg->dim.x * cfg->dim.y * cfg->dim.z; dimlen.w = fieldlen; @@ -3700,6 +3724,10 @@ is more than what your have specified (%d), please use the -H option to specify } } + if (cfg->extrasrclen && cfg->srcid < 0) { // when multiple sources are simulated, the total photons are evenly divided + scale[0] *= (cfg->extrasrclen + 1); + } + cfg->normalizer = scale[0]; cfg->his.normalizer = scale[0]; diff --git a/src/mcx_core.h b/src/mcx_core.h index b23be05d..bbde57f4 100644 --- a/src/mcx_core.h +++ b/src/mcx_core.h @@ -145,6 +145,7 @@ typedef struct __align__(16) KernelParams { float Rtstep; /**< reciprocal of the step size */ MCXSrc src; /**< additional source data, including pos, dir, param1, param2 */ unsigned int extrasrclen; /**< number of additional sources */ + int srcid; /**< 0: simulate all sources combined, -1: simulate all sources separately; >0: only simulate a single source */ float4 s0; /**< initial stokes parameters, for polarized photon simulation */ float3 maxidx; /**< maximum index in x/y/z directions for out-of-bound tests */ uint4 dimlen; /**< maximum index used to convert x/y/z to 1D array index */ diff --git a/src/mcx_utils.c b/src/mcx_utils.c index a556e1f8..5607e048 100644 --- a/src/mcx_utils.c +++ b/src/mcx_utils.c @@ -104,7 +104,7 @@ const char shortopt[] = {'h', 'i', 'f', 'n', 't', 'T', 's', 'a', 'g', 'b', '-', 'z', 'u', 'H', 'P', 'd', 'r', 'S', 'p', 'e', 'U', 'R', 'l', 'L', '-', 'I', '-', 'G', 'M', 'A', 'E', 'v', 'D', 'k', 'q', 'Y', 'O', 'F', '-', '-', 'x', 'X', '-', 'K', 'm', 'V', 'B', 'W', 'w', '-', - '-', '-', 'Z', 'j', '-', '\0' + '-', '-', 'Z', 'j', '-', '-', '\0' }; /** @@ -123,7 +123,8 @@ const char* fullopt[] = {"--help", "--interactive", "--input", "--photon", "--replaydet", "--outputtype", "--outputformat", "--maxjumpdebug", "--maxvoidstep", "--saveexit", "--saveref", "--gscatter", "--mediabyte", "--momentum", "--specular", "--bc", "--workload", "--savedetflag", - "--internalsrc", "--bench", "--dumpjson", "--zip", "--json", "--atomic", "" + "--internalsrc", "--bench", "--dumpjson", "--zip", "--json", "--atomic", + "--srcid", "" }; /** @@ -332,7 +333,9 @@ void mcx_initcfg(Config* cfg) { cfg->invcdf = NULL; cfg->nangle = 0; cfg->angleinvcdf = NULL; + cfg->srcid = 0; cfg->extrasrclen = 0; + cfg->his.totalsource = cfg->extrasrclen + 1; cfg->srcdata = NULL; memset(cfg->jsonfile, 0, MAX_PATH_LENGTH); memset(cfg->bc, 0, 13); @@ -583,7 +586,7 @@ void mcx_savenii(float* dat, size_t len, char* name, int type32bit, int outputfo * @param[in] cfg: simulation configuration */ -void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, Config* cfg) { +void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, int iscol, Config* cfg) { FILE* fp; char fname[MAX_FULL_PATH] = {'\0'}; int affine[] = {0, 0, 1, 0, 0, 0}; @@ -704,7 +707,7 @@ void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name /* the "NIFTIData" section stores volumetric data */ ubjw_begin_object(root, UBJ_MIXED, 0); - if (mcx_jdataencode(vol, ndim, dims, (isfloat ? "single" : "uint32"), 4, cfg->zipid, root, 1, cfg)) { + if (mcx_jdataencode(vol, ndim, dims, (isfloat ? "single" : "uint32"), 4, cfg->zipid, root, 1, iscol, cfg)) { MCX_ERROR(-1, "error when converting to JSON"); } @@ -746,7 +749,7 @@ void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name * @param[in] cfg: simulation configuration */ -void mcx_savejnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, Config* cfg) { +void mcx_savejnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, int iscol, Config* cfg) { FILE* fp; char fname[MAX_FULL_PATH] = {'\0'}; int affine[] = {0, 0, 1, 0, 0, 0}; @@ -831,7 +834,7 @@ void mcx_savejnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name /* the "NIFTIData" section stores volumetric data */ cJSON_AddItemToObject(root, "NIFTIData", dat = cJSON_CreateObject()); - if (mcx_jdataencode(vol, ndim, dims, (isfloat ? "single" : "uint32"), 4, cfg->zipid, dat, 0, cfg)) { + if (mcx_jdataencode(vol, ndim, dims, (isfloat ? "single" : "uint32"), 4, cfg->zipid, dat, 0, iscol, cfg)) { MCX_ERROR(-1, "error when converting to JSON"); } @@ -887,48 +890,25 @@ void mcx_savedata(float* dat, size_t len, Config* cfg) { mcx_savenii(dat, len * (1 + (cfg->outputtype == otRF)), name, NIFTI_TYPE_FLOAT32, cfg->outputformat, cfg); return; } else if (cfg->outputformat == ofJNifti || cfg->outputformat == ofBJNifti) { - int d1 = (cfg->maxgate == 1); + uint dims[6] = {cfg->dim.x, cfg->dim.y, cfg->dim.z, cfg->maxgate, cfg->srcnum, 1}; + float voxelsize[6] = {cfg->steps.x, cfg->steps.y, cfg->steps.z, cfg->tstep, 1, 1}; - if (cfg->seed == SEED_FROM_FILE && ((cfg->replaydet == -1 && cfg->detnum > 1) || cfg->srcnum > 1 || cfg->outputtype == otRF)) { - if (cfg->outputtype == otRF) { - uint dims[6] = {2, cfg->detnum* cfg->srcnum, cfg->maxgate, cfg->dim.z, cfg->dim.y, cfg->dim.x}; - float voxelsize[] = {1, 1, cfg->tstep, cfg->steps.z, cfg->steps.y, cfg->steps.x}; - - if (cfg->outputformat == ofJNifti) { - mcx_savejnii(dat, 6, dims, voxelsize, name, 1, cfg); - } else { - mcx_savebnii(dat, 6, dims, voxelsize, name, 1, cfg); - } - } else { - uint dims[5] = {cfg->detnum* cfg->srcnum, cfg->maxgate, cfg->dim.z, cfg->dim.y, cfg->dim.x}; - float voxelsize[] = {1, cfg->tstep, cfg->steps.z, cfg->steps.y, cfg->steps.x}; + if (cfg->extrasrclen && cfg->srcid == -1) { + dims[5] = cfg->extrasrclen + 1; + } - if (cfg->outputformat == ofJNifti) { - mcx_savejnii(dat, 5, dims, voxelsize, name, 1, cfg); - } else { - mcx_savebnii(dat, 5, dims, voxelsize, name, 1, cfg); - } - } - } else { - uint dims[] = {cfg->dim.x, cfg->dim.y, cfg->dim.z, cfg->maxgate}; - float voxelsize[] = {cfg->steps.x, cfg->steps.y, cfg->steps.z, cfg->tstep}; - size_t datalen = cfg->dim.x * cfg->dim.y * cfg->dim.z * cfg->maxgate; - uint* buf = (uint*)malloc(datalen * sizeof(float)); - memcpy(buf, dat, datalen * sizeof(float)); - - if (d1) { - mcx_convertcol2row(&buf, (uint3*)dims); - } else { - mcx_convertcol2row4d(&buf, (uint4*)dims); - } + if (cfg->seed == SEED_FROM_FILE && (cfg->replaydet == -1 && cfg->detnum > 1)) { + dims[5] *= cfg->detnum; + } - if (cfg->outputformat == ofJNifti) { - mcx_savejnii((float*)buf, 4 - d1, dims, voxelsize, name, 1, cfg); - } else { - mcx_savebnii((float*)buf, 4 - d1, dims, voxelsize, name, 1, cfg); - } + if (cfg->seed == SEED_FROM_FILE && cfg->outputtype == otRF) { + dims[5] *= 2; + } - free(buf); + if (cfg->outputformat == ofJNifti) { + mcx_savejnii(dat, 5 + (dims[5] > 1), dims, voxelsize, name, 1, 1, cfg); + } else { + mcx_savebnii(dat, 5 + (dims[5] > 1), dims, voxelsize, name, 1, 1, cfg); } return; @@ -1032,6 +1012,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c cJSON_AddNumberToObject(hdr, "Repeat", cfg->his.respin); cJSON_AddNumberToObject(hdr, "SrcNum", cfg->his.srcnum); cJSON_AddNumberToObject(hdr, "SaveDetFlag", cfg->his.savedetflag); + cJSON_AddNumberToObject(hdr, "TotalSource", cfg->his.totalsource); cJSON_AddItemToObject(hdr, "Media", sub = cJSON_CreateArray()); for (int i = 0; i < cfg->medianum; i++) { @@ -1059,7 +1040,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c cJSON_AddItemToObject(dat, dname[id], sub = cJSON_CreateObject()); - if (mcx_jdataencode(buf, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, cfg)) { + if (mcx_jdataencode(buf, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) { MCX_ERROR(-1, "error when converting to JSON"); } @@ -1101,7 +1082,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c cJSON_AddItemToObject(dat, dname[id], sub = cJSON_CreateObject()); - if (mcx_jdataencode(val, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, cfg)) { + if (mcx_jdataencode(val, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) { MCX_ERROR(-1, "error when converting to JSON"); } @@ -1115,7 +1096,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c uint dims[2] = {count, cfg->his.seedbyte}; cJSON_AddItemToObject(dat, "seed", sub = cJSON_CreateObject()); - if (mcx_jdataencode(seeds, 2, dims, "uint8", 1, cfg->zipid, sub, 0, cfg)) { + if (mcx_jdataencode(seeds, 2, dims, "uint8", 1, cfg->zipid, sub, 0, 0, cfg)) { MCX_ERROR(-1, "error when converting to JSON"); } } @@ -1610,6 +1591,15 @@ void mcx_preprocess(Config* cfg) { } if (cfg->extrasrclen) { + if (cfg->srcnum > 1) { + MCX_ERROR(-4, "simulating multiple sources currently can not be used with photon-sharing"); + } + + if (cfg->srcid > (int)cfg->extrasrclen + 1) { + printf("cfg->srcid=%d\n", cfg->srcid); + MCX_ERROR(-4, "srcid exceeds total defined source count"); + } + for (int i = 0; i < cfg->extrasrclen; i++) { if (cfg->srcdata[i].srcpos.w == 0.f) { cfg->srcdata[i].srcpos.w = 1.f; @@ -2093,6 +2083,7 @@ void mcx_loadconfig(FILE* in, Config* cfg) { cfg->his.detnum = cfg->detnum; cfg->his.srcnum = cfg->srcnum; cfg->his.savedetflag = cfg->savedetflag; + cfg->his.totalsource = ((cfg->srcid <= 0) ? cfg->extrasrclen + 1 : 1); //cfg->his.colcount=2+(cfg->medianum-1)*(2+(cfg->ismomentum>0))+(cfg->issaveexit>0)*6; /*column count=maxmedia+2*/ if (in == stdin) { @@ -2667,7 +2658,12 @@ int mcx_loadjson(cJSON* root, Config* cfg) { cfg->srcpos.z--; /*convert to C index, grid center*/ } - cfg->srctype = mcx_keylookup((char*)FIND_JSON_KEY("Type", "Optode.Source.Type", src, "pencil", valuestring), srctypeid); + subitem = FIND_JSON_OBJ("Type", "Optode.Source.Type", src); + + if (subitem) { + cfg->srctype = mcx_keylookup((char*)FIND_JSON_KEY("Type", "Optode.Source.Type", src, "pencil", valuestring), srctypeid); + } + subitem = FIND_JSON_OBJ("Param1", "Optode.Source.Param1", src); if (subitem && cJSON_IsArray(subitem)) { @@ -2784,10 +2780,23 @@ int mcx_loadjson(cJSON* root, Config* cfg) { } } - cfg->omega = FIND_JSON_KEY("Frequency", "Optode.Source.Frequency", src, 0.f, valuedouble); - cfg->omega *= TWO_PI; - cfg->lambda = FIND_JSON_KEY("WaveLength", "Optode.Source.WaveLength", src, 0.f, valuedouble); - cfg->srcnum = FIND_JSON_KEY("SrcNum", "Optode.Source.SrcNum", src, cfg->srcnum, valueint); + if (FIND_JSON_OBJ("Frequency", "Optode.Source.Frequency", src)) { + cfg->omega = FIND_JSON_KEY("Frequency", "Optode.Source.Frequency", src, 0.f, valuedouble); + cfg->omega *= TWO_PI; + } + + if (FIND_JSON_OBJ("WaveLength", "Optode.Source.WaveLength", src)) { + cfg->lambda = FIND_JSON_KEY("WaveLength", "Optode.Source.WaveLength", src, 0.f, valuedouble); + } + + if (FIND_JSON_OBJ("SrcNum", "Optode.Source.SrcNum", src)) { + cfg->srcnum = FIND_JSON_KEY("SrcNum", "Optode.Source.SrcNum", src, cfg->srcnum, valueint); + } + + if (FIND_JSON_OBJ("ID", "Optode.Source.ID", src)) { + cfg->srcid = FIND_JSON_KEY("ID", "Optode.Source.ID", src, cfg->srcid, valueint); + } + subitem = FIND_JSON_OBJ("Pattern", "Optode.Source.Pattern", src); if (subitem) { @@ -3110,6 +3119,7 @@ int mcx_loadjson(cJSON* root, Config* cfg) { cfg->his.detnum = cfg->detnum; cfg->his.srcnum = cfg->srcnum; cfg->his.savedetflag = cfg->savedetflag; + cfg->his.totalsource = ((cfg->srcid <= 0) ? cfg->extrasrclen + 1 : 1); //cfg->his.colcount=2+(cfg->medianum-1)*(2+(cfg->ismomentum>0))+(cfg->issaveexit>0)*6; /*column count=maxmedia+2*/ return 0; } @@ -3259,7 +3269,7 @@ void mcx_savejdata(char* filename, Config* cfg) { dims[2] = cfg->srcparam2.w; cJSON_AddItemToObject(sub, "Pattern", tmp = cJSON_CreateObject()); - int ret = mcx_jdataencode(cfg->srcpattern, 2 + (cfg->srcnum > 1), dims + (cfg->srcnum == 1), "single", dims[0] * dims[1] * dims[2], cfg->zipid, tmp, 0, cfg); + int ret = mcx_jdataencode(cfg->srcpattern, 2 + (cfg->srcnum > 1), dims + (cfg->srcnum == 1), "single", dims[0] * dims[1] * dims[2], cfg->zipid, tmp, 0, 0, cfg); if (ret) { MCX_ERROR(ret, "data compression or base64 encoding failed"); @@ -3322,7 +3332,7 @@ void mcx_savejdata(char* filename, Config* cfg) { obj = cJSON_CreateObject(); ret = mcx_jdataencode(buf, 3, &cfg->dim.x, (cfg->mediabyte == 1 ? "uint8" : (cfg->mediabyte == 2 ? "uint16" : "uint32")), - outputbytes / datalen, cfg->zipid, obj, 0, cfg); + outputbytes / datalen, cfg->zipid, obj, 0, 0, cfg); if (buf) { free(buf); @@ -4259,9 +4269,9 @@ void mcx_dumpmask(Config* cfg) { mcx_convertcol2row((uint**)(&buf), (uint3*)dims); if (cfg->outputformat == ofJNifti) { - mcx_savejnii((float*)buf, 3, dims, voxelsize, fname, 0, cfg); + mcx_savejnii((float*)buf, 3, dims, voxelsize, fname, 0, 0, cfg); } else { - mcx_savebnii((float*)buf, 3, dims, voxelsize, fname, 0, cfg); + mcx_savebnii((float*)buf, 3, dims, voxelsize, fname, 0, 0, cfg); } free(buf); @@ -4375,10 +4385,12 @@ int mcx_jdatadecode(void** vol, int* ndim, uint* dims, int maxdim, char** type, * @param[in] type: a string of JData data types, such as "uint8" "float32", "int32" etc * @param[in] byte: number of byte per voxel * @param[in] zipid: zip method: 0:zlib,1:gzip,2:base64,3:lzma,4:lzip,5:lz4,6:lz4hc - * @param[in] obj: a pre-created cJSON object to store the output JData fields + * @param[in] obj: a pre-created cJSON or UBJ object to store the output JData fields + * @param[in] isubj: 1 if obj is a binary JSON (UBJ) object, 0 if obj is a cJSON object + * @param[in] cfg: mcx config struct */ -int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int zipid, void* obj, int isubj, Config* cfg) { +int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int zipid, void* obj, int isubj, int iscol, Config* cfg) { uint datalen = 1; size_t compressedbytes, totalbytes; uchar* compressed = NULL, *buf = NULL; @@ -4407,6 +4419,11 @@ int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int UBJ_WRITE_KEY(item, "_ArrayType_", string, type); ubjw_write_key(item, "_ArraySize_"); UBJ_WRITE_ARRAY(item, uint32, ndim, dims); + + if (iscol) { + UBJ_WRITE_KEY(item, "_ArrayOrder_", string, "c"); + } + UBJ_WRITE_KEY(item, "_ArrayZipType_", string, zipformat[zipid]); UBJ_WRITE_KEY(item, "_ArrayZipSize_", uint32, datalen); ubjw_write_key(item, "_ArrayZipData_"); @@ -4423,6 +4440,11 @@ int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int if (!ret) { cJSON_AddStringToObject((cJSON*)obj, "_ArrayType_", type); cJSON_AddItemToObject((cJSON*)obj, "_ArraySize_", cJSON_CreateIntArray((int*)dims, ndim)); + + if (iscol) { + cJSON_AddStringToObject((cJSON*)obj, "_ArrayOrder_", "c"); + } + cJSON_AddStringToObject((cJSON*)obj, "_ArrayZipType_", zipformat[zipid]); cJSON_AddNumberToObject((cJSON*)obj, "_ArrayZipSize_", datalen); cJSON_AddStringToObject((cJSON*)obj, "_ArrayZipData_", (char*)buf); @@ -4989,6 +5011,8 @@ void mcx_parsecmd(int argc, char* argv[], Config* cfg) { int isatomic = 1; i = mcx_readarg(argc, argv, i, &(isatomic), "char"); cfg->sradius = (isatomic) ? -2.f : 0.f; + } else if (strcmp(argv[i] + 2, "srcid") == 0) { + i = mcx_readarg(argc, argv, i, &(cfg->srcid), "int"); } else if (strcmp(argv[i] + 2, "internalsrc") == 0) { i = mcx_readarg(argc, argv, i, &(cfg->internalsrc), "int"); } else { @@ -5299,7 +5323,7 @@ void mcx_printheader(Config* cfg) { # MCX proudly developed human-readable JSON-based data formats for easy reuse,#\n\ # Please consider using JSON (" S_BLUE "https://neurojson.org/" S_MAGENTA ") for your research data #\n\ ###############################################################################\n\ -$Rev:: $" S_GREEN MCX_VERSION S_MAGENTA "$Date:: $ by $Author:: $\n\ +$Rev:: $" S_GREEN MCX_VERSION S_MAGENTA " $Date:: $ by $Author:: $\n\ ###############################################################################\n" S_RESET); } @@ -5497,6 +5521,8 @@ where possible parameters include (the first value in [*|*] is the default)\n\ --gscatter [1e9|int] after a photon completes the specified number of\n\ scattering events, mcx then ignores anisotropy g\n\ and only performs isotropic scattering for speed\n\ + --srcid [0|-1,0,1,2,..] -1 simulate multi-source separately;0 all sources\n\ + together; a positive integer runs a single source\n\ --internalsrc [0|1] set to 1 to skip entry search to speedup launch\n\ --maxvoidstep [1000|int] maximum distance (in voxel unit) of a photon that\n\ can travel before entering the domain, if \n\ diff --git a/src/mcx_utils.h b/src/mcx_utils.h index 8ade1903..46cce1d2 100644 --- a/src/mcx_utils.h +++ b/src/mcx_utils.h @@ -121,8 +121,9 @@ typedef struct MCXHistoryHeader { float normalizer; /**< what is the normalization factor */ int respin; /**< if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset */ unsigned int srcnum; /**< number of sources for simultaneous pattern sources */ - unsigned int savedetflag; /**< number of sources for simultaneous pattern sources */ - int reserved[2]; /**< reserved fields for future extension */ + unsigned int savedetflag; /**< bit-mask for storage of different types of detected photon data */ + unsigned int totalsource; /**< total source number when multiple sources are defined */ + int reserved[1]; /**< reserved fields for future extension */ } History; /** @@ -271,6 +272,7 @@ typedef struct MCXConfig { float* invcdf; /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle */ unsigned int nangle; /**< number of samples for inverse-cdf of launch angle, will be added by 2 to include -1 and 1 on the two ends */ float* angleinvcdf; /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle of photon launch */ + int srcid; /**< flag to control the simulation of multiple sources */ unsigned int extrasrclen; /**< length of additional sources */ ExtraSrc* srcdata; /**< buffer to store multiple source input data */ } Config; @@ -317,10 +319,10 @@ void mcx_flush(Config* cfg); int mcx_run_from_json(char* jsonstr); float mcx_updatemua(unsigned int mediaid, Config* cfg); void mcx_savejdata(char* filename, Config* cfg); -int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int zipid, void* obj, int isubj, Config* cfg); +int mcx_jdataencode(void* vol, int ndim, uint* dims, char* type, int byte, int zipid, void* obj, int isubj, int iscol, Config* cfg); int mcx_jdatadecode(void** vol, int* ndim, uint* dims, int maxdim, char** type, cJSON* obj, Config* cfg); -void mcx_savejnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, Config* cfg); -void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, Config* cfg); +void mcx_savejnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, int iscol, Config* cfg); +void mcx_savebnii(float* vol, int ndim, uint* dims, float* voxelsize, char* name, int isfloat, int iscol, Config* cfg); void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* cfg); int mcx_svmc_bgvoxel(int vol); void mcx_loadseedjdat(char* filename, Config* cfg); diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp index 44dda611..8f76879f 100644 --- a/src/mcxlab.cpp +++ b/src/mcxlab.cpp @@ -261,6 +261,10 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { fieldlen *= 2; } + if (cfg.extrasrclen && cfg.srcid == -1) { + fieldlen *= (cfg.extrasrclen + 1); + } + cfg.exportfield = (float*)calloc(fieldlen, sizeof(float)); } @@ -392,6 +396,10 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { fielddim[5] = 2; } + if (cfg.extrasrclen && cfg.srcid == -1) { + fielddim[5] *= (cfg.extrasrclen + 1); + } + fieldlen = fielddim[0] * fielddim[1] * fielddim[2] * fielddim[3] * fielddim[4] * fielddim[5]; if (cfg.issaveref) { @@ -561,6 +569,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf GET_ONE_FIELD(cfg, maxjumpdebug) GET_ONE_FIELD(cfg, gscatter) GET_ONE_FIELD(cfg, srcnum) + GET_ONE_FIELD(cfg, srcid) GET_ONE_FIELD(cfg, omega) GET_ONE_FIELD(cfg, issave2pt) GET_ONE_FIELD(cfg, lambda) diff --git a/src/pmcx.cpp b/src/pmcx.cpp index 96174a15..24025220 100644 --- a/src/pmcx.cpp +++ b/src/pmcx.cpp @@ -427,6 +427,7 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) { GET_SCALAR_FIELD(user_cfg, mcx_config, maxjumpdebug, py::int_); GET_SCALAR_FIELD(user_cfg, mcx_config, gscatter, py::int_); GET_SCALAR_FIELD(user_cfg, mcx_config, srcnum, py::int_); + GET_SCALAR_FIELD(user_cfg, mcx_config, srcid, py::int_); GET_SCALAR_FIELD(user_cfg, mcx_config, omega, py::float_); GET_SCALAR_FIELD(user_cfg, mcx_config, lambda, py::float_); GET_VEC3_FIELD(user_cfg, mcx_config, steps, float); @@ -1073,6 +1074,10 @@ py::dict pmcx_interface(const py::dict& user_cfg) { field_len *= 2; } + if (mcx_config.extrasrclen && mcx_config.srcid == -1) { + field_len *= (mcx_config.extrasrclen + 1); + } + mcx_config.exportfield = (float*) calloc(field_len, sizeof(float)); } @@ -1198,6 +1203,10 @@ py::dict pmcx_interface(const py::dict& user_cfg) { field_dim[5] = 2; } + if (mcx_config.extrasrclen && mcx_config.srcid == -1) { + field_dim[5] *= (mcx_config.extrasrclen + 1); + } + field_len = field_dim[0] * field_dim[1] * field_dim[2] * field_dim[3] * field_dim[4] * field_dim[5]; std::vector array_dims;