Skip to content

Commit

Permalink
remove all file IO functions from mex and oct files
Browse files Browse the repository at this point in the history
  • Loading branch information
fangq committed Sep 12, 2023
1 parent 4ca899f commit 130570b
Show file tree
Hide file tree
Showing 7 changed files with 154 additions and 83 deletions.
10 changes: 9 additions & 1 deletion commons/Makefile_common.mk
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@ AROUTPUT += -o
MAKE ?= make

ZMATLIB :=lib/libzmat.a
USERARFLAGS?=$(ZMATLIB)

LIBOPENCLDIR ?= /usr/local/cuda/lib64
LIBOPENCL ?=-lOpenCL
Expand All @@ -46,6 +45,14 @@ CUDA_STATIC=--cudart static -Xcompiler "-static-libgcc -static-libstdc++"
ECHO := echo
MKDIR := mkdir

ifneq (,$(filter $(MAKECMDGOALS),mex oct mexomp octomp cudamex cudaoct))
ZMATLIB=
else
FILES+=mmc cjson/cJSON ubj/ubjw
endif

USERARFLAGS?=$(ZMATLIB)

MEXLINKLIBS=-L"\$$MATLABROOT/extern/lib/\$$ARCH" -L"\$$MATLABROOT/bin/\$$ARCH" -lmx -lmex $(ZMATLIB)

ARCH = $(shell uname -m)
Expand Down Expand Up @@ -183,6 +190,7 @@ web: EXTRALIB=-s SIMD=1 -s WASM=1 -s EXTRA_EXPORTED_RUNTIME_METHODS='["cwrap"]'

mex oct mexomp octomp: EXTRALIB=
mex oct mexomp octomp: CCFLAGS+=$(DLLFLAG) -DMCX_CONTAINER
mex oct mexomp octomp: CUCCOPT+=-DMCX_CONTAINER
mex oct mexomp octomp: CPPFLAGS+=-g $(DLLFLAG) -DMCX_CONTAINER
mex oct mexomp octomp: BINDIR=../mmclab
mex mexomp: AR=$(MKMEX)
Expand Down
2 changes: 1 addition & 1 deletion src/Makefile
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
ROOTDIR = ..
BINARY=mmc

FILES=mmc_rand_xorshift128p mmc_mesh mmc_raytrace mmc_utils mmc_tictoc mmc cjson/cJSON mmc_host mmc_highorder mmc_cl_utils mmc_cl_host ubj/ubjw
FILES=mmc_rand_xorshift128p mmc_mesh mmc_raytrace mmc_utils mmc_tictoc mmc_host mmc_highorder mmc_cl_utils mmc_cl_host
CLPROGRAM=mmc_core

ifeq ($(findstring Darwin,$(PLATFORM)), Darwin)
Expand Down
4 changes: 4 additions & 0 deletions src/mmc_cl_host.c
Original file line number Diff line number Diff line change
Expand Up @@ -742,6 +742,8 @@ is more than what your have specified (%d), please use the -H option to specify
mesh_normalize(mesh, cfg, cfg->energyabs, cfg->energytot, 0);
}

#ifndef MCX_CONTAINER

if (cfg->issave2pt && cfg->parentid == mpStandalone) {
MMC_FPRINTF(cfg->flog, "saving data to file ...\t");
mesh_saveweight(mesh, cfg, 0);
Expand All @@ -761,6 +763,8 @@ is more than what your have specified (%d), please use the -H option to specify
mesh_saveweight(mesh, cfg, 1);
}

#endif

// total energy here equals total simulated photons+unfinished photons for all threads
MMC_FPRINTF(cfg->flog, "simulated %zu photons (%zu) with %d devices (ray-tet %.0f)\nMCX simulation speed: %.2f photon/ms\n",
cfg->nphoton, cfg->nphoton, workdev, reporter.raytet, (double)cfg->nphoton / toc);
Expand Down
4 changes: 4 additions & 0 deletions src/mmc_cu_host.cu
Original file line number Diff line number Diff line change
Expand Up @@ -769,6 +769,8 @@ void mmc_run_simulation(mcconfig* cfg, tetmesh* mesh, raytracer* tracer, GPUInfo
mesh_normalize(mesh, cfg, cfg->energyabs, cfg->energytot, 0);
}

#ifndef MCX_CONTAINER

if (cfg->issave2pt && cfg->parentid == mpStandalone) {
MMC_FPRINTF(cfg->flog, "saving data to file ...\t");
mesh_saveweight(mesh, cfg, 0);
Expand All @@ -791,6 +793,8 @@ void mmc_run_simulation(mcconfig* cfg, tetmesh* mesh, raytracer* tracer, GPUInfo
mesh_saveweight(mesh, cfg, 1);
}

#endif

// total energy here equals total simulated photons+unfinished photons for
// all threads
MMC_FPRINTF(cfg->flog,
Expand Down
23 changes: 21 additions & 2 deletions src/mmc_host.c
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,9 @@ necessary data structures. This include the command line options (cfg),
tetrahedral mesh (mesh) and the ray-tracer precomputed data (tracer).
*******************************************************************************/

#ifndef MCX_CONTAINER


/**
* \brief Initialize simulation configuration structure cfg using command line options
*
Expand Down Expand Up @@ -89,6 +92,8 @@ int mmc_init_from_json(mcconfig* cfg, tetmesh* mesh, raytracer* tracer, char* jc
return 0;
}

#endif

/**
* \brief Rest simulation related data structures
*
Expand Down Expand Up @@ -338,7 +343,9 @@ int mmc_run_mp(mcconfig* cfg, tetmesh* mesh, raytracer* tracer) {
cfg->his.normalizer = sum_normalizer / cfg->srcnum; // average normalizer value for all simulated sources
}

if (cfg->issave2pt) {
#ifndef MCX_CONTAINER

if (cfg->issave2pt && cfg->parentid == mpStandalone) {
switch (cfg->outputtype) {
case otFlux:
MMCDEBUG(cfg, dlTime, (cfg->flog, "saving flux ..."));
Expand All @@ -356,17 +363,25 @@ int mmc_run_mp(mcconfig* cfg, tetmesh* mesh, raytracer* tracer) {
mesh_saveweight(mesh, cfg, 0);
}

if (cfg->issavedet) {
#endif

if (cfg->issavedet && cfg->parentid == mpStandalone) {
MMCDEBUG(cfg, dlTime, (cfg->flog, "saving detected photons ..."));

#ifndef MCX_CONTAINER

if (cfg->issaveexit) {
mesh_savedetphoton(master.partialpath, master.photonseed, master.bufpos, (sizeof(RandType)*RAND_BUF_LEN), cfg);
}

#endif

if (cfg->issaveexit == 2) {
float* detimage = (float*)calloc(cfg->detparam1.w * cfg->detparam2.w * cfg->maxgate, sizeof(float));
mesh_getdetimage(detimage, master.partialpath, master.bufpos, cfg, mesh);
#ifndef MCX_CONTAINER
mesh_savedetimage(detimage, cfg);
#endif
free(detimage);
}

Expand All @@ -379,11 +394,15 @@ int mmc_run_mp(mcconfig* cfg, tetmesh* mesh, raytracer* tracer) {
}
}

#ifndef MCX_CONTAINER

if (cfg->issaveref) {
MMCDEBUG(cfg, dlTime, (cfg->flog, "saving surface diffuse reflectance ..."));
mesh_saveweight(mesh, cfg, 1);
}

#endif

MMCDEBUG(cfg, dlTime, (cfg->flog, "\tdone\t%d\n", GetTimeMillis() - t0));
visitor_clear(&master);

Expand Down
176 changes: 97 additions & 79 deletions src/mmc_mesh.c
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,8 @@ void mesh_init(tetmesh* mesh) {
mesh->nmax.w = 1.f;
}

#ifndef MCX_CONTAINER

/**
* @brief Loading user-specified mesh data
*
Expand All @@ -167,6 +169,8 @@ void mesh_init_from_cfg(tetmesh* mesh, mcconfig* cfg) {
}
}

#endif

/**
* @brief Error-handling in mesh operations
*
Expand Down Expand Up @@ -203,6 +207,90 @@ void mesh_filenames(const char* format, char* foutput, mcconfig* cfg) {
}
}


void mesh_createdualmesh(tetmesh* mesh, mcconfig* cfg) {
int i;
mesh->nmin.x = VERY_BIG;
mesh->nmin.y = VERY_BIG;
mesh->nmin.z = VERY_BIG;
mesh->nmax.x = -VERY_BIG;
mesh->nmax.y = -VERY_BIG;
mesh->nmax.z = -VERY_BIG;

for (i = 0; i < mesh->nn; i++) {
mesh->nmin.x = MIN(mesh->node[i].x, mesh->nmin.x);
mesh->nmin.y = MIN(mesh->node[i].y, mesh->nmin.y);
mesh->nmin.z = MIN(mesh->node[i].z, mesh->nmin.z);
mesh->nmax.x = MAX(mesh->node[i].x, mesh->nmax.x);
mesh->nmax.y = MAX(mesh->node[i].y, mesh->nmax.y);
mesh->nmax.z = MAX(mesh->node[i].z, mesh->nmax.z);
}

mesh->nmin.x -= EPS;
mesh->nmin.y -= EPS;
mesh->nmin.z -= EPS;
mesh->nmax.x += EPS;
mesh->nmax.y += EPS;
mesh->nmax.z += EPS;

cfg->dim.x = (int)((mesh->nmax.x - mesh->nmin.x) / cfg->steps.x) + 1;
cfg->dim.y = (int)((mesh->nmax.y - mesh->nmin.y) / cfg->steps.y) + 1;
cfg->dim.z = (int)((mesh->nmax.z - mesh->nmin.z) / cfg->steps.z) + 1;

cfg->crop0.x = cfg->dim.x;
cfg->crop0.y = cfg->dim.y * cfg->dim.x;
cfg->crop0.z = cfg->dim.y * cfg->dim.x * cfg->dim.z;
}

/**
* @brief Identify wide-field source and detector-related elements (type=-1 for source, type=-2 for det)
*
* @param[in] mesh: the mesh object
* @param[in] cfg: the simulation configuration structure
*/

void mesh_srcdetelem(tetmesh* mesh, mcconfig* cfg) {
int i;

mesh->srcelemlen = 0;
mesh->detelemlen = 0;

for (i = 0; i < mesh->ne; i++) {
if (mesh->type[i] == -1) { /*number of elements in the initial candidate list*/
mesh->srcelemlen++;
cfg->e0 = (cfg->e0 == 0) ? i + 1 : cfg->e0;
}

if (mesh->type[i] == -2) { /*number of elements in the initial candidate list*/
mesh->detelemlen++;
cfg->isextdet = 1;
cfg->detnum = 0; // when detecting wide-field detectors, suppress point detectors
}
}

/*Record the index of inital elements to initiate source search*/
/*Then change the type of initial elements back to 0 to continue propogation*/
if (mesh->srcelemlen > 0 || mesh->detelemlen > 0) {
int is = 0, id = 0;
mesh->srcelem = (int*)calloc(mesh->srcelemlen, sizeof(int));
mesh->detelem = (int*)calloc(mesh->detelemlen, sizeof(int));

for (i = 0; i < mesh->ne; i++) {
if (mesh->type[i] < 0) {
if (mesh->type[i] == -1) {
mesh->srcelem[is++] = i + 1;
mesh->type[i] = 0;
} else if (mesh->type[i] == -2) { /*keep -2, will be replaced to medianum+1 in loadmedia*/
mesh->detelem[id++] = i + 1;
}
}
}
}
}


#ifndef MCX_CONTAINER

/**
* @brief Load node file and initialize the related mesh properties
*
Expand Down Expand Up @@ -241,40 +329,6 @@ void mesh_loadnode(tetmesh* mesh, mcconfig* cfg) {
}
}

void mesh_createdualmesh(tetmesh* mesh, mcconfig* cfg) {
int i;
mesh->nmin.x = VERY_BIG;
mesh->nmin.y = VERY_BIG;
mesh->nmin.z = VERY_BIG;
mesh->nmax.x = -VERY_BIG;
mesh->nmax.y = -VERY_BIG;
mesh->nmax.z = -VERY_BIG;

for (i = 0; i < mesh->nn; i++) {
mesh->nmin.x = MIN(mesh->node[i].x, mesh->nmin.x);
mesh->nmin.y = MIN(mesh->node[i].y, mesh->nmin.y);
mesh->nmin.z = MIN(mesh->node[i].z, mesh->nmin.z);
mesh->nmax.x = MAX(mesh->node[i].x, mesh->nmax.x);
mesh->nmax.y = MAX(mesh->node[i].y, mesh->nmax.y);
mesh->nmax.z = MAX(mesh->node[i].z, mesh->nmax.z);
}

mesh->nmin.x -= EPS;
mesh->nmin.y -= EPS;
mesh->nmin.z -= EPS;
mesh->nmax.x += EPS;
mesh->nmax.y += EPS;
mesh->nmax.z += EPS;

cfg->dim.x = (int)((mesh->nmax.x - mesh->nmin.x) / cfg->steps.x) + 1;
cfg->dim.y = (int)((mesh->nmax.y - mesh->nmin.y) / cfg->steps.y) + 1;
cfg->dim.z = (int)((mesh->nmax.z - mesh->nmin.z) / cfg->steps.z) + 1;

cfg->crop0.x = cfg->dim.x;
cfg->crop0.y = cfg->dim.y * cfg->dim.x;
cfg->crop0.z = cfg->dim.y * cfg->dim.x * cfg->dim.z;
}

/**
* @brief Load optical property file and initialize the related mesh properties
*
Expand Down Expand Up @@ -452,51 +506,6 @@ void mesh_loadelem(tetmesh* mesh, mcconfig* cfg) {
mesh_srcdetelem(mesh, cfg);
}

/**
* @brief Identify wide-field source and detector-related elements (type=-1 for source, type=-2 for det)
*
* @param[in] mesh: the mesh object
* @param[in] cfg: the simulation configuration structure
*/

void mesh_srcdetelem(tetmesh* mesh, mcconfig* cfg) {
int i;

mesh->srcelemlen = 0;
mesh->detelemlen = 0;

for (i = 0; i < mesh->ne; i++) {
if (mesh->type[i] == -1) { /*number of elements in the initial candidate list*/
mesh->srcelemlen++;
cfg->e0 = (cfg->e0 == 0) ? i + 1 : cfg->e0;
}

if (mesh->type[i] == -2) { /*number of elements in the initial candidate list*/
mesh->detelemlen++;
cfg->isextdet = 1;
cfg->detnum = 0; // when detecting wide-field detectors, suppress point detectors
}
}

/*Record the index of inital elements to initiate source search*/
/*Then change the type of initial elements back to 0 to continue propogation*/
if (mesh->srcelemlen > 0 || mesh->detelemlen > 0) {
int is = 0, id = 0;
mesh->srcelem = (int*)calloc(mesh->srcelemlen, sizeof(int));
mesh->detelem = (int*)calloc(mesh->detelemlen, sizeof(int));

for (i = 0; i < mesh->ne; i++) {
if (mesh->type[i] < 0) {
if (mesh->type[i] == -1) {
mesh->srcelem[is++] = i + 1;
mesh->type[i] = 0;
} else if (mesh->type[i] == -2) { /*keep -2, will be replaced to medianum+1 in loadmedia*/
mesh->detelem[id++] = i + 1;
}
}
}
}
}

/**
* @brief Load tet element volume file and initialize the related mesh properties
Expand Down Expand Up @@ -677,6 +686,8 @@ void mesh_loadseedfile(tetmesh* mesh, mcconfig* cfg) {
fclose(fp);
}

#endif

/**
* @brief Clearing the mesh data structure
*
Expand Down Expand Up @@ -1214,6 +1225,8 @@ float mc_next_scatter(float g, float3* dir, RandType* ran, RandType* ran0, mccon
return nextslen;
}

#ifndef MCX_CONTAINER

/**
* @brief Save the fluence output to a file
*
Expand Down Expand Up @@ -1340,6 +1353,7 @@ void mesh_savedetphoton(float* ppath, void* seeds, int count, int seedbyte, mcco
fclose(fp);
}

#endif

/**
* @brief Save binned detected photon data over an area-detector as time-resolved 2D images
Expand Down Expand Up @@ -1403,6 +1417,8 @@ void mesh_getdetimage(float* detmap, float* ppath, int count, mcconfig* cfg, tet
}
}

#ifndef MCX_CONTAINER

/**
* @brief Save binned detected photon data over an area-detector
*
Expand Down Expand Up @@ -1431,6 +1447,8 @@ void mesh_savedetimage(float* detmap, mcconfig* cfg) {
fclose(fp);
}

#endif

/**
* @brief Recompute the detected photon weight from the partial-pathlengths
*
Expand Down
Loading

0 comments on commit 130570b

Please sign in to comment.