diff --git a/.github/workflows/ubuntu-basic.yml b/.github/workflows/ubuntu-basic.yml index 831e3bbaf..c5836ea08 100644 --- a/.github/workflows/ubuntu-basic.yml +++ b/.github/workflows/ubuntu-basic.yml @@ -55,11 +55,9 @@ jobs: --with-scalapack-dir="${ISSM_DIR}/externalpackages/petsc/install" \ --with-mumps-dir="${ISSM_DIR}/externalpackages/petsc/install" \ --with-triangle-dir="${ISSM_DIR}/externalpackages/triangle/install" \ - --with-semic-dir="${ISSM_DIR}/externalpackages/semic/install" \ - --with-m1qn3-dir="${ISSM_DIR}/externalpackages/m1qn3/install" + --with-semic-dir="${ISSM_DIR}/externalpackages/semic/install" ext_pkg_build_command: | cd $ISSM_DIR/externalpackages/triangle && ./install-linux.sh - cd $ISSM_DIR/externalpackages/m1qn3 && ./install-linux.sh cd $ISSM_DIR/externalpackages/petsc && ./install-3.22-linux.sh cd $ISSM_DIR/externalpackages/semic && ./install.sh test_cases: '["101:399", "401:899"]' diff --git a/.github/workflows/ubuntu-codipack.yml b/.github/workflows/ubuntu-codipack.yml index 0a1fdae22..e8995a6c5 100644 --- a/.github/workflows/ubuntu-codipack.yml +++ b/.github/workflows/ubuntu-codipack.yml @@ -59,12 +59,10 @@ jobs: --with-mumps-dir="${ISSM_DIR}/externalpackages/petsc/install" \ --with-triangle-dir="${ISSM_DIR}/externalpackages/triangle/install" \ --with-gsl-dir="${ISSM_DIR}/externalpackages/gsl/install" \ - --with-m1qn3-dir="${ISSM_DIR}/externalpackages/m1qn3/install" \ --with-medipack-dir="${ISSM_DIR}/externalpackages/medipack/install" \ --with-codipack-dir="${ISSM_DIR}/externalpackages/codipack/install" ext_pkg_build_command: | cd $ISSM_DIR/externalpackages/triangle && ./install-linux.sh - cd $ISSM_DIR/externalpackages/m1qn3 && ./install-linux.sh cd $ISSM_DIR/externalpackages/petsc && ./install-3.22-linux.sh cd $ISSM_DIR/externalpackages/gsl && ./install.sh cd $ISSM_DIR/externalpackages/codipack && ./install.sh diff --git a/.github/workflows/ubuntu-python.yml b/.github/workflows/ubuntu-python.yml index e68758489..fc3ecb564 100644 --- a/.github/workflows/ubuntu-python.yml +++ b/.github/workflows/ubuntu-python.yml @@ -59,11 +59,9 @@ jobs: --with-scalapack-dir="${ISSM_DIR}/externalpackages/petsc/install" \ --with-mumps-dir="${ISSM_DIR}/externalpackages/petsc/install" \ --with-triangle-dir="${ISSM_DIR}/externalpackages/triangle/install" \ - --with-semic-dir="${ISSM_DIR}/externalpackages/semic/install" \ - --with-m1qn3-dir="${ISSM_DIR}/externalpackages/m1qn3/install" + --with-semic-dir="${ISSM_DIR}/externalpackages/semic/install" ext_pkg_build_command: | cd $ISSM_DIR/externalpackages/triangle && ./install-linux.sh - cd $ISSM_DIR/externalpackages/m1qn3 && ./install-linux.sh cd $ISSM_DIR/externalpackages/petsc && ./install-3.22-linux.sh cd $ISSM_DIR/externalpackages/semic && ./install.sh test_cases: '["101:399", "401:899"]' diff --git a/src/c/classes/Elements/Element.cpp b/src/c/classes/Elements/Element.cpp index 77bc516fb..bccad9e4a 100644 --- a/src/c/classes/Elements/Element.cpp +++ b/src/c/classes/Elements/Element.cpp @@ -22,10 +22,6 @@ #include "../Inputs/ControlInput.h" #include "../Inputs/ArrayInput.h" #include "../Inputs/IntArrayInput.h" -#ifdef _HAVE_PyBind11_ - #include - namespace py=pybind11; -#endif /*}}}*/ #define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/ @@ -4209,132 +4205,6 @@ void Element::PositiveDegreeDayGCM(){/*{{{*/ delete gauss; } /*}}}*/ -#if _HAVE_PyBind11_ -void Element::SmbEmulator(IssmDouble timeinputs, EmulatorParam* emulator){/*{{{*/ - - int numvertices = this->GetNumberOfVertices(); - - IssmDouble* smb = xNew(numvertices); - IssmDouble* elev = xNew(numvertices); - IssmDouble* lon = xNew(numvertices); - IssmDouble* lat = xNew(numvertices); - IssmDouble* al = xNew(numvertices); - IssmDouble* st = xNew(numvertices); - IssmDouble* tt = xNew(numvertices); - IssmDouble* swd = xNew(numvertices); - IssmDouble* lwd = xNew(numvertices); - IssmDouble* swu = xNew(numvertices); - IssmDouble* lwu = xNew(numvertices); - IssmDouble* shf = xNew(numvertices); - IssmDouble* lhf = xNew(numvertices); - - Input* elev_input = this->GetInput(SurfaceEnum); _assert_(elev_input); - Input* al_input = this->GetInput(MariaAlEnum,timeinputs); _assert_(al_input); - Input* st_input = this->GetInput(MariaStEnum,timeinputs); _assert_(st_input); - Input* tt_input = this->GetInput(MariaTtEnum,timeinputs); _assert_(tt_input); - Input* swd_input = this->GetInput(MariaSwdEnum,timeinputs); _assert_(swd_input); - Input* lwd_input = this->GetInput(MariaLwdEnum,timeinputs); _assert_(lwd_input); - Input* swu_input = this->GetInput(MariaSwuEnum,timeinputs); _assert_(swu_input); - Input* lwu_input = this->GetInput(MariaLwuEnum,timeinputs); _assert_(lwu_input); - Input* shf_input = this->GetInput(MariaShfEnum,timeinputs); _assert_(shf_input); - Input* lhf_input = this->GetInput(MariaLhfEnum,timeinputs); _assert_(lhf_input); - - for(int iv=0;ivvertices[iv]->GetLongitude(); - lat[iv] = this->vertices[iv]->GetLatitude(); - } - this->GetInputListOnVertices(elev,elev_input,0.); - this->GetInputListOnVertices(al,al_input,0.); - this->GetInputListOnVertices(st,st_input,0.); - this->GetInputListOnVertices(tt,tt_input,0.); - this->GetInputListOnVertices(swd,swd_input,0.); - this->GetInputListOnVertices(lwd,lwd_input,0.); - this->GetInputListOnVertices(swu,swu_input,0.); - this->GetInputListOnVertices(lwu,lwu_input,0.); - this->GetInputListOnVertices(shf,shf_input,0.); - this->GetInputListOnVertices(lhf,lhf_input,0.); - - try{ - py::gil_scoped_acquire gil; - - py::array_t lon_np(numvertices); - py::array_t lat_np(numvertices); - py::array_t elev_np(numvertices); - py::array_t al_np(numvertices); - py::array_t st_np(numvertices); - py::array_t tt_np(numvertices); - py::array_t swd_np(numvertices); - py::array_t lwd_np(numvertices); - py::array_t swu_np(numvertices); - py::array_t lwu_np(numvertices); - py::array_t shf_np(numvertices); - py::array_t lhf_np(numvertices); - - auto lon_view = lon_np.mutable_unchecked<1>(); - auto lat_view = lat_np.mutable_unchecked<1>(); - auto elev_view = elev_np.mutable_unchecked<1>(); - auto al_view = al_np.mutable_unchecked<1>(); - auto st_view = st_np.mutable_unchecked<1>(); - auto tt_view = tt_np.mutable_unchecked<1>(); - auto swd_view = swd_np.mutable_unchecked<1>(); - auto lwd_view = lwd_np.mutable_unchecked<1>(); - auto swu_view = swu_np.mutable_unchecked<1>(); - auto lwu_view = lwu_np.mutable_unchecked<1>(); - auto shf_view = shf_np.mutable_unchecked<1>(); - auto lhf_view = lhf_np.mutable_unchecked<1>(); - - for(int iv=0;iv(lon[iv]); - lat_view(iv) = static_cast(lat[iv]); - elev_view(iv) = static_cast(elev[iv]); - al_view(iv) = static_cast(al[iv]); - st_view(iv) = static_cast(st[iv]); - tt_view(iv) = static_cast(tt[iv]); - swd_view(iv) = static_cast(swd[iv]); - lwd_view(iv) = static_cast(lwd[iv]); - swu_view(iv) = static_cast(swu[iv]); - lwu_view(iv) = static_cast(lwu[iv]); - shf_view(iv) = static_cast(shf[iv]); - lhf_view(iv) = static_cast(lhf[iv]); - } - - py::object pred_obj = emulator->mod.attr("predict_smb_np")( - elev_np, lon_np, lat_np, al_np, st_np, tt_np, swd_np, lwd_np, swu_np, lwu_np, shf_np, lhf_np, - py::arg("dtype") = "float64"); - py::array_t pred(pred_obj); - auto pred_view = pred.unchecked<1>(); - _assert_(pred.shape(0)==numvertices); - - for(int iv=0;iv(pred_view(iv)); - } - } - catch(const py::error_already_set& e){ - _error_(std::string("Python SMB emulator inference failed: ") + e.what()); - } - catch(const std::exception& e){ - _error_(std::string("SMB emulator inference failed: ") + e.what()); - } - - this->AddInput(SmbMassBalanceEnum,smb,P1Enum); - - xDelete(lon); - xDelete(lat); - xDelete(elev); - xDelete(al); - xDelete(st); - xDelete(tt); - xDelete(swd); - xDelete(lwd); - xDelete(swu); - xDelete(lwu); - xDelete(shf); - xDelete(lhf); - xDelete(smb); - -} -/*}}}*/ -#endif void Element::ProjectGridDataToMesh(IssmDouble* griddata,IssmDouble* x_grid,IssmDouble* y_grid,int Nx,int Ny,int input_enum){/*{{{*/ const int NUM_VERTICES = this->GetNumberOfVertices(); diff --git a/src/c/classes/Elements/Element.h b/src/c/classes/Elements/Element.h index 7fd4344cc..2507d72d9 100644 --- a/src/c/classes/Elements/Element.h +++ b/src/c/classes/Elements/Element.h @@ -39,9 +39,6 @@ template class Vector; class ElementMatrix; class ElementVector; class BarystaticContributions; -#if _HAVE_PyBind11_ -class EmulatorParam; -#endif /*}}}*/ class Element: public Object{ @@ -185,9 +182,6 @@ class Element: public Object{ void PositiveDegreeDaySicopolis(bool isfirnwarming); void PositiveDegreeDayFast(bool isfirnwarming); void PositiveDegreeDayGCM(); - #ifdef _HAVE_PyBind11_ - void SmbEmulator(IssmDouble timeinputs, EmulatorParam* emulator); - #endif void ProjectGridDataToMesh(IssmDouble* griddata,IssmDouble* x_grid,IssmDouble* y_grid,int Nx,int Ny,int input_enum); void SmbDebrisEvatt(); void RignotMeltParameterization(); diff --git a/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp b/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp index b433fcc5c..c3e8a5b9e 100644 --- a/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp +++ b/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp @@ -10,6 +10,14 @@ #include "../../classes/Inputs/TransientInput.h" #include "../../shared/Random/random.h" +#ifdef _HAVE_PyBind11_ +#include +#include +#include +#include +namespace py = pybind11; +#endif + void SmbForcingx(FemModel* femmodel){/*{{{*/ // void SmbForcingx(smb,ni){ @@ -425,11 +433,174 @@ void SmbEmulatorx(FemModel* femmodel){/*{{{*/ IssmDouble timeinputs; femmodel->parameters->FindParam(&timeinputs,TimeEnum); + IssmDouble rho_ice; + femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); + const IssmDouble dts = 86400.; + + int numberofvertices_local = femmodel->vertices->NumberOfVerticesLocalAll(); + if(numberofvertices_local==0) return; + + std::vector lid_to_batch(numberofvertices_local,-1); + std::vector batch_lids; + std::vector elev; + std::vector lon; + std::vector lat; + std::vector al; + std::vector st; + std::vector tt; + std::vector swd; + std::vector lwd; + std::vector swu; + std::vector lwu; + std::vector shf; + std::vector lhf; + std::vector maria_elev; + std::vector maria_smb; + batch_lids.reserve(numberofvertices_local); + elev.reserve(numberofvertices_local); + lon.reserve(numberofvertices_local); + lat.reserve(numberofvertices_local); + al.reserve(numberofvertices_local); + st.reserve(numberofvertices_local); + tt.reserve(numberofvertices_local); + swd.reserve(numberofvertices_local); + lwd.reserve(numberofvertices_local); + swu.reserve(numberofvertices_local); + lwu.reserve(numberofvertices_local); + shf.reserve(numberofvertices_local); + lhf.reserve(numberofvertices_local); + maria_elev.reserve(numberofvertices_local); + maria_smb.reserve(numberofvertices_local); for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SmbEmulator(timeinputs, smbemulator); + const int numvertices = element->GetNumberOfVertices(); + + std::vector elev_e(numvertices); + std::vector al_e(numvertices); + std::vector st_e(numvertices); + std::vector tt_e(numvertices); + std::vector swd_e(numvertices); + std::vector lwd_e(numvertices); + std::vector swu_e(numvertices); + std::vector lwu_e(numvertices); + std::vector shf_e(numvertices); + std::vector lhf_e(numvertices); + std::vector maria_elev_e(numvertices); + std::vector maria_smb_e(numvertices); + + Input* elev_input = element->GetInput(SurfaceEnum); _assert_(elev_input); + Input* maria_elev_input = element->GetInput(MariaElevEnum,timeinputs); _assert_(maria_elev_input); + Input* maria_smb_input = element->GetInput(MariaSmbEnum,timeinputs); _assert_(maria_smb_input); + Input* al_input = element->GetInput(MariaAlEnum,timeinputs); _assert_(al_input); + Input* st_input = element->GetInput(MariaStEnum,timeinputs); _assert_(st_input); + Input* tt_input = element->GetInput(MariaTtEnum,timeinputs); _assert_(tt_input); + Input* swd_input = element->GetInput(MariaSwdEnum,timeinputs); _assert_(swd_input); + Input* lwd_input = element->GetInput(MariaLwdEnum,timeinputs); _assert_(lwd_input); + Input* swu_input = element->GetInput(MariaSwuEnum,timeinputs); _assert_(swu_input); + Input* lwu_input = element->GetInput(MariaLwuEnum,timeinputs); _assert_(lwu_input); + Input* shf_input = element->GetInput(MariaShfEnum,timeinputs); _assert_(shf_input); + Input* lhf_input = element->GetInput(MariaLhfEnum,timeinputs); _assert_(lhf_input); + + element->GetInputListOnVertices(elev_e.data(),elev_input,0.); + element->GetInputListOnVertices(maria_elev_e.data(),maria_elev_input,0.); + element->GetInputListOnVertices(maria_smb_e.data(),maria_smb_input,0.); + element->GetInputListOnVertices(al_e.data(),al_input,0.); + element->GetInputListOnVertices(st_e.data(),st_input,0.); + element->GetInputListOnVertices(tt_e.data(),tt_input,0.); + element->GetInputListOnVertices(swd_e.data(),swd_input,0.); + element->GetInputListOnVertices(lwd_e.data(),lwd_input,0.); + element->GetInputListOnVertices(swu_e.data(),swu_input,0.); + element->GetInputListOnVertices(lwu_e.data(),lwu_input,0.); + element->GetInputListOnVertices(shf_e.data(),shf_input,0.); + element->GetInputListOnVertices(lhf_e.data(),lhf_input,0.); + + for(int iv=0;ivvertices[iv]->Lid(); + if(lid<0 || lid>=numberofvertices_local){ + _error_("SmbEmulatorx vertex local id " << lid << " is outside local vertex range [0," << numberofvertices_local << ")"); + } + if(lid_to_batch[lid]>=0) continue; + + lid_to_batch[lid] = static_cast(batch_lids.size()); + batch_lids.push_back(lid); + elev.push_back(static_cast(elev_e[iv])); + lon.push_back(static_cast(element->vertices[iv]->GetLongitude())); + lat.push_back(static_cast(element->vertices[iv]->GetLatitude())); + al.push_back(static_cast(al_e[iv])); + st.push_back(static_cast(st_e[iv])); + tt.push_back(static_cast(tt_e[iv])); + swd.push_back(static_cast(swd_e[iv])); + lwd.push_back(static_cast(lwd_e[iv])); + swu.push_back(static_cast(swu_e[iv])); + lwu.push_back(static_cast(lwu_e[iv])); + shf.push_back(static_cast(shf_e[iv])); + lhf.push_back(static_cast(lhf_e[iv])); + maria_elev.push_back(static_cast(maria_elev_e[iv])); + maria_smb.push_back(static_cast(maria_smb_e[iv])); + } + } + + int batch_size = static_cast(batch_lids.size()); + if(batch_size==0) return; + + std::vector smb_local(numberofvertices_local,0.); + + try{ + py::gil_scoped_acquire gil; + + py::array_t elev_np(batch_size,elev.data()); + py::array_t lon_np(batch_size,lon.data()); + py::array_t lat_np(batch_size,lat.data()); + py::array_t al_np(batch_size,al.data()); + py::array_t st_np(batch_size,st.data()); + py::array_t tt_np(batch_size,tt.data()); + py::array_t swd_np(batch_size,swd.data()); + py::array_t lwd_np(batch_size,lwd.data()); + py::array_t swu_np(batch_size,swu.data()); + py::array_t lwu_np(batch_size,lwu.data()); + py::array_t shf_np(batch_size,shf.data()); + py::array_t lhf_np(batch_size,lhf.data()); + + py::object pred_obj = smbemulator->mod.attr("predict_smb_np")( + elev_np, lon_np, lat_np, al_np, st_np, tt_np, swd_np, lwd_np, swu_np, lwu_np, shf_np, lhf_np, + py::arg("dtype") = "float64"); + py::array_t pred(pred_obj); + py::buffer_info pred_info = pred.request(); + + if(pred_info.ndim!=1){ + _error_("SmbEmulatorx output is not a 1D array"); + } + if(pred_info.shape[0]!=batch_size){ + std::ostringstream message; + message << "SmbEmulatorx unexpected output length: " << pred_info.shape[0] << ", expected " << batch_size; + _error_(message.str()); + } + + const double* pred_data = static_cast(pred_info.ptr); + for(int i=0;i(pred_data[i])/(rho_ice*dts); + } + if(IssmComm::GetRank()==0){ + int nprint = batch_size<3 ? batch_size : 3; + _printf0_(" SMB emulator diagnostic at time " << timeinputs << ":\n"); + for(int i=0;i=6) executable='issm_dakota.exe'; end + dakota_version_str = IssmConfig('_DAKOTA_VERSION_'); version = str2num(dakota_version_str(1:3)); + if version >= 6, executable = 'issm_dakota.exe'; end end if isoceancoupling - executable='issm_ocean.exe'; + executable = 'issm_ocean.exe'; end if ~ispc() - % Check that executable exists at the right path - if ~exist([cluster.codepath '/' executable],'file') - error(['File ' cluster.codepath '/' executable ' does not exist']); + % Verify the executable exists + exepath = [cluster.codepath '/' executable]; + if ~exist(exepath, 'file') + error('File %s does not exist', exepath); end - % Process codepath and prepend empty spaces with \ to avoid errors in queuing script - codepath=strrep(cluster.codepath,' ','\ '); + % Escape spaces in codepath for the shell script + codepath = strrep(cluster.codepath, ' ', '\ '); execpath = [cluster.executionpath '/' dirname]; - % Write queuing script - fid=fopen(filename, 'w'); - fprintf(fid,'#!%s\n',cluster.shell); + % Build the mpi prefix once (empty string when MPI is not available) + mpiprefix = ''; + if IssmConfig('_HAVE_MPI_') + mpiprefix = sprintf('mpiexec -np %i ', cluster.np); + end + + % Build the core command string if isvalgrind - %Add --gen-suppressions=all to get suppression lines - %fprintf(fid,'LD_PRELOAD=%s \\\n',cluster.valgrindlib); it could be deleted if ismac - if IssmConfig('_HAVE_MPI_') - fprintf(fid,'mpiexec -np %i %s --leak-check=full --leak-check=full --show-leak-kinds=all --error-limit=no --dsymutil=yes --suppressions=%s %s/%s %s %s %s 2> %s.errlog > %s.outlog ',... - cluster.np,cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution,execpath, modelname,modelname,modelname); - else - fprintf(fid,'%s --leak-check=full --dsymutil=yes --error-limit=no --leak-check=full --show-leak-kinds=all --suppressions=%s %s/%s %s %s %s 2> %s.errlog > %s.outlog',... - cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution, execpath, modelname,modelname,modelname); - end + vgflags = '--leak-check=full --show-leak-kinds=all --error-limit=no --dsymutil=yes'; else - if IssmConfig('_HAVE_MPI_') - fprintf(fid,'mpiexec -np %i %s --leak-check=full --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog > %s.outlog',... - cluster.np,cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution, execpath,modelname,modelname,modelname); - else - fprintf(fid,'%s --leak-check=full --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog > %s.outlog',... - cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution, execpath, modelname,modelname,modelname); - end + vgflags = '--leak-check=full --error-limit=no'; end + cmd = sprintf('%s%s %s --suppressions=%s %s/%s %s %s %s 2> %s.errlog > %s.outlog', ... + mpiprefix, cluster.valgrind, vgflags, cluster.valgrindsup, ... + codepath, executable, solution, execpath, modelname, modelname, modelname); + elseif isgprof - fprintf(fid,'\n gprof %s/issm.exe gmon.out > %s.performance',cluster.codepath,modelname); + cmd = sprintf('gprof %s/issm.exe gmon.out > %s.performance', cluster.codepath, modelname); + + elseif cluster.interactive + cmd = sprintf('%s%s/%s %s %s %s', mpiprefix, codepath, executable, solution, execpath, modelname); + else - if cluster.interactive - if IssmConfig('_HAVE_MPI_') - fprintf(fid,'mpiexec -np %i %s/%s %s %s %s\n',cluster.np,cluster.codepath,executable,solution, execpath, modelname); - else - fprintf(fid,'%s/%s %s %s %s',cluster.codepath,executable,solution, execpath, modelname); - end + % Non-interactive: redirect output and run in background + if IssmConfig('_HAVE_MPI_') + cmd = sprintf('%s%s/%s %s %s %s 2> %s/%s.errlog > %s/%s.outlog &', ... + mpiprefix, codepath, executable, solution, execpath, modelname, execpath, modelname, execpath, modelname); else - if IssmConfig('_HAVE_MPI_') - fprintf(fid,'mpiexec -np %i %s/%s %s %s %s 2> %s/%s.errlog > %s/%s.outlog &',cluster.np,cluster.codepath,executable,solution,execpath,modelname,execpath,modelname,execpath,modelname); - else - fprintf(fid,'%s/%s %s %s %s 2> %s.errlog > %s.outlog &',cluster.codepath,executable,solution, execpath,modelname,modelname,modelname); - end + cmd = sprintf('%s/%s %s %s %s 2> %s.errlog > %s.outlog &', ... + codepath, executable, solution, execpath, modelname, modelname, modelname); end end - if ~io_gather, %concatenate the output files: - fprintf(fid,'\ncat %s.outbin.* > %s.outbin',modelname,modelname); + + % Write the queue script + fid = fopen(filename, 'w'); + fprintf(fid, '#!%s\n', cluster.shell); + fprintf(fid, '%s', cmd); + if ~io_gather + % Concatenate distributed output files into one + fprintf(fid, '\ncat %s.outbin.* > %s.outbin', modelname, modelname); end fclose(fid); else % Windows - batfilename=[filename(1:end-6) '.bat']; - fid=fopen(batfilename,'w'); - fprintf(fid,'@echo off\n'); - execdir=[cluster.executionpath '\' dirname]; - if cluster.np>1 - fprintf(fid,'"C:\\Program Files\\Microsoft MPI\\Bin\\mpiexec.exe" -n %i "%s\\%s" %s "%s" %s',cluster.np,cluster.codepath,executable,solution,execdir,modelname); + batfilename = [filename(1:end-6) '.bat']; + execdir = [cluster.executionpath '\' dirname]; + fid = fopen(batfilename, 'w'); + fprintf(fid, '@echo off\n'); + if cluster.np > 1 + fprintf(fid, '"C:\\Program Files\\Microsoft MPI\\Bin\\mpiexec.exe" -n %i "%s\\%s" %s "%s" %s', ... + cluster.np, cluster.codepath, executable, solution, execdir, modelname); else - fprintf(fid,'"%s\\%s" %s "%s" %s',cluster.codepath,executable,solution,execdir,modelname); + fprintf(fid, '"%s\\%s" %s "%s" %s', cluster.codepath, executable, solution, execdir, modelname); end fclose(fid); + end end %}}} diff --git a/src/m/classes/clusters/generic.py b/src/m/classes/clusters/generic.py index 179d96bf8..eaf06819e 100644 --- a/src/m/classes/clusters/generic.py +++ b/src/m/classes/clusters/generic.py @@ -89,7 +89,7 @@ def checkconsistency(self, md, solution, analyses): # {{{ def BuildQueueScript(self, md, filename): # {{{ - # Get variables from md + # Unpack fields used below dirname = md.private.runtimename modelname = md.miscellaneous.name solution = md.private.solution @@ -99,62 +99,63 @@ def BuildQueueScript(self, md, filename): # {{{ isdakota = md.qmu.isdakota isoceancoupling = md.transient.isoceancoupling - # Which executable are we calling? - executable = 'issm.exe' # default + # Determine which executable to call + executable = 'issm.exe' if isdakota: - version = IssmConfig('_DAKOTA_VERSION_') - version = float(version[0]) - if version >= 6: executable = 'issm_dakota.exe' + if float(IssmConfig('_DAKOTA_VERSION_')[0]) >= 6: + executable = 'issm_dakota.exe' if isoceancoupling: executable = 'issm_ocean.exe' - # Write queuing script if not ispc(): - fid = open(filename, 'w') - fid.write('#!/bin/sh\n') - if not isvalgrind: - if self.interactive: - if IssmConfig('_HAVE_MPI_')[0]: - fid.write('mpiexec -np {} {}/{} {} {}/{} {}'.format(self.np, self.codepath, executable, solution, self.executionpath, dirname, modelname)) - else: - fid.write('{}/{} {} {}/{} {}'.format(self.codepath, executable, solution, self.executionpath, dirname, modelname)) - else: - if IssmConfig('_HAVE_MPI_')[0]: - fid.write('mpiexec -np {} {}/{} {} {}/{} {} 2> {}.errlog > {}.outlog'. - format(self.np, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) - else: - fid.write('{}/{} {} {}/{} {} 2> {}.errlog > {}.outlog '. - format(self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) + + # Verify the executable exists + exepath = '{}/{}'.format(self.codepath, executable) + if not os.path.isfile(exepath): + raise RuntimeError('File {} does not exist'.format(exepath)) + + execpath = '{}/{}'.format(self.executionpath, dirname) + mpiprefix = 'mpiexec -np {} '.format(self.np) if IssmConfig('_HAVE_MPI_')[0] else '' + + # Build the core command string + if isvalgrind: + supstring = ' '.join('--suppressions=' + s for s in self.valgrindsup) + cmd = '{}{} --leak-check=full {} {}/{} {} {} {} 2> {}.errlog > {}.outlog'.format( + mpiprefix, self.valgrind, supstring, + self.codepath, executable, solution, execpath, modelname, + modelname, modelname) + elif isgprof: - fid.write('\n gprof {}/{} gmon.out > {}.performance'.format(self.codepath, executable, modelname)) - else: - #Add --gen -suppressions = all to get suppression lines - #fid.write('LD_PRELOAD={} \\\n'.format(self.valgrindlib)) it could be deleted - supstring = '' - for supfile in self.valgrindsup: - supstring += ' --suppressions=' + supfile - - if IssmConfig('_HAVE_MPI_')[0]: - fid.write('mpiexec -np {} {} --leak-check=full {} {}/{} {} {}/{} {} 2> {}.errlog > {}.outlog '. - format(self.np, self.valgrind, supstring, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) - else: - fid.write('{} --leak-check=full {} {}/{} {} {}/{} {} 2> {}.errlog > {}.outlog '. - format(self.valgrind, supstring, self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) + cmd = 'gprof {}/{} gmon.out > {}.performance'.format(self.codepath, executable, modelname) - if not io_gather: #concatenate the output files: - fid.write('\ncat {}.outbin .*>{}.outbin'.format(modelname, modelname)) - fid.close() + elif self.interactive: + cmd = '{}{}/{} {} {} {}'.format(mpiprefix, self.codepath, executable, solution, execpath, modelname) + + else: + # Non-interactive: redirect output and run in background + cmd = '{}{}/{} {} {} {} 2> {}.errlog > {}.outlog &'.format( + mpiprefix, self.codepath, executable, solution, execpath, modelname, + modelname, modelname) + + # Write the queue script + with open(filename, 'w') as fid: + fid.write('#!{}\n'.format(self.shell)) + fid.write(cmd) + if not io_gather: + # Concatenate distributed output files into one + fid.write('\ncat {}.outbin.* > {}.outbin'.format(modelname, modelname)) else: # Windows + batfilename = filename[:-6] + '.bat' - fid = open(batfilename, 'w') - fid.write('@echo off\n') - if self.interactive: - fid.write('"{}/{}" {} "{}/{}" {} '.format(self.codepath, executable, solution, self.executionpath, dirname, modelname)) - else: - fid.write('"{}/{}" {} "{}/{}" {} 2>{}.errlog>{}.outlog'. - format(self.codepath, executable, solution, self.executionpath, dirname, modelname, modelname, modelname)) - fid.close() + execdir = '{}/{}'.format(self.executionpath, dirname) + with open(batfilename, 'w') as fid: + fid.write('@echo off\n') + if self.interactive: + fid.write('"{}/{}" {} "{}" {}'.format(self.codepath, executable, solution, execdir, modelname)) + else: + fid.write('"{}/{}" {} "{}" {} 2> {}.errlog > {}.outlog'.format( + self.codepath, executable, solution, execdir, modelname, modelname, modelname)) # }}} def BuildKrigingQueueScript(self, md, filename): # {{{ diff --git a/src/m/classes/clusters/unity.m b/src/m/classes/clusters/unity.m new file mode 100644 index 000000000..c33f30ac3 --- /dev/null +++ b/src/m/classes/clusters/unity.m @@ -0,0 +1,129 @@ +%UNITY (Massachusetts Green High Performance Computing Center) cluster class definition +% +% Usage: +% cluster=unity(); +% cluster=unity('np',3); +% cluster=unity('np',3,'login','username'); + +classdef unity + properties (SetAccess=public) + % {{{ + name = 'unity' + login = ''; + numnodes = 1; + cpuspernode = 16; + codepath = ''; + executionpath = ''; + time = 10; %in hours + memory = 32; %in Gb + email = ''; + end + %}}} + methods + function cluster=unity(varargin) % {{{ + + %initialize cluster using default settings if provided + if (exist('unity_settings')==2), unity_settings; end + + %use provided options to change fields + cluster=AssignObjectFields(pairoptions(varargin{:}),cluster); + end + %}}} + function disp(cluster) % {{{ + % display the object + disp(sprintf('class ''%s'' object ''%s'' = ',class(cluster),inputname(1))); + disp(sprintf(' name: %s',cluster.name)); + disp(sprintf(' login: %s',cluster.login)); + disp(sprintf(' numnodes: %i',cluster.numnodes)); + disp(sprintf(' cpuspernode: %i',cluster.cpuspernode)); + disp(sprintf(' time: %i hours',cluster.time)); + disp(sprintf(' memory: %i Gb',cluster.memory)); + disp(sprintf(' email: %s (receive notifications if END,FAIL)',cluster.email)); + disp(sprintf(' codepath: %s',cluster.codepath)); + disp(sprintf(' executionpath: %s',cluster.executionpath)); + end + %}}} + function numprocs=nprocs(cluster) % {{{ + %compute number of processors + numprocs=cluster.numnodes*cluster.cpuspernode; + end + %}}} + function md = checkconsistency(cluster,md,solution,analyses) % {{{ + %Miscellaneous + if isempty(cluster.login), md = checkmessage(md,'login empty'); end + if isempty(cluster.codepath), md = checkmessage(md,'codepath empty'); end + if isempty(cluster.executionpath), md = checkmessage(md,'executionpath empty'); end + end + %}}} + function BuildQueueScript(cluster, md, filename) % {{{ + + %Get variables from md + dirname = md.private.runtimename; + modelname = md.miscellaneous.name; + solution = md.private.solution; + io_gather = md.settings.io_gather; + isdakota = md.qmu.isdakota; + isoceancoupling = md.transient.isoceancoupling; + + %write queuing script + fid=fopen(filename, 'w'); + fprintf(fid,'#!/bin/bash -l\n'); + fprintf(fid,'#SBATCH --job-name=%s\n',modelname); + fprintf(fid,'#SBATCH -p cpu # Partition\n'); + fprintf(fid,'#SBATCH -o %s.outlog \n',modelname); + fprintf(fid,'#SBATCH -e %s.errlog \n',modelname); + fprintf(fid,'#SBATCH --nodes=%i\n',cluster.numnodes); + fprintf(fid,'#SBATCH --ntasks-per-node=%i\n',cluster.cpuspernode); + fprintf(fid,'#SBATCH --time=%s\n',eraseBetween(datestr(cluster.time/24,'dd-HH:MM:SS'),1,1)); %walltime is in d-HH:MM:SS format. cluster.time is in hour + fprintf(fid,'#SBATCH --mem=%iG\n',cluster.memory); + if ~isempty(cluster.email) + fprintf(fid,'#SBATCH --mail-type=END,FAIL\n'); + fprintf(fid,'#SBATCH --mail-user=%s\n', cluster.email); + end + fprintf(fid,'\n'); + fprintf(fid,'module load intel-oneapi-compilers/2024.1.0 intel-oneapi-mpi/2021.12.1 petsc/3.22.1\n'); + fprintf(fid,'mpiexec -n %i %s/issm.exe %s %s %s\n',cluster.nprocs(), cluster.codepath,solution,[cluster.executionpath '/' dirname],modelname); + if ~io_gather, %concatenate the output files: + fprintf(fid,'cat %s.outbin.* > %s.outbin',modelname,modelname); + end + fclose(fid); + end %}}} + function UploadQueueJob(cluster,modelname,dirname,filelist) % {{{ + + %compress the files into one zip. + %filelist contains full paths; tar with -C so only basenames are stored in the archive + root=[issmdir() '/execution/' dirname]; + compressstring=['tar -C ' root ' -zcf ' dirname '.tar.gz']; + for i=1:numel(filelist) + if ~exist(filelist{i},'file') + error(['File ' filelist{i} ' not found']); + end + [~,fname,fext]=fileparts(filelist{i}); + compressstring=[compressstring ' ' fname fext]; + end + system(compressstring); + + %upload input files + issmscpout(cluster.name,cluster.executionpath,cluster.login,0,{[dirname '.tar.gz']}); + + end %}}} + function LaunchQueueJob(cluster,modelname,dirname,filelist,restart,batch) % {{{ + + %Execute Queue job + if ~isempty(restart) + launchcommand=['cd ' cluster.executionpath ' && cd ' dirname ' && hostname && sbatch ' modelname '.queue ']; + else + launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' dirname ' && mkdir ' dirname ... + ' && cd ' dirname ' && mv ../' dirname '.tar.gz ./ && tar -zxf ' dirname '.tar.gz && sbatch ' modelname '.queue ']; + end + issmssh(cluster.name,cluster.login,0,launchcommand); + end %}}} + function Download(cluster,dirname,filelist) % {{{ + + %copy files from cluster to current directory + directory=[cluster.executionpath '/' dirname '/']; + issmscpin(cluster.name,cluster.login,0,directory,filelist, 2); %use {} and not \{\} + + end %}}} + end +end diff --git a/src/m/classes/linearbasalforcings.py b/src/m/classes/linearbasalforcings.py index 59970924c..30d62c39b 100644 --- a/src/m/classes/linearbasalforcings.py +++ b/src/m/classes/linearbasalforcings.py @@ -16,7 +16,6 @@ class linearbasalforcings(object): def __init__(self, *args): # {{{ nargs = len(args) if nargs == 0: - print('empty init') self.deepwater_melting_rate = 0. self.deepwater_elevation = 0. self.upperwater_melting_rate = 0. diff --git a/src/m/classes/m1qn3inversion.py b/src/m/classes/m1qn3inversion.py index 0a03475f2..43cc7d437 100644 --- a/src/m/classes/m1qn3inversion.py +++ b/src/m/classes/m1qn3inversion.py @@ -17,7 +17,6 @@ class m1qn3inversion(object): def __init__(self, *args): # {{{ if not len(args): - print('empty init') self.iscontrol = 0 self.incomplete_adjoint = 0 self.control_parameters = np.nan diff --git a/src/m/netcdf/OLD/export_netCDF.m b/src/m/netcdf/OLD/export_netCDF.m new file mode 100644 index 000000000..ff94a69cb --- /dev/null +++ b/src/m/netcdf/OLD/export_netCDF.m @@ -0,0 +1,510 @@ +function export_netCDF(md,filename) +%verbosity of the code, 0 is no messages, 5 is chatty + verbose = 5; + if exist(filename), + delete(filename) + % disp(sprintf('File %s allready exist', filename)); + % prompt = 'Give a new name or "delete" to replace: '; + % newname = input(prompt,'s'); + % if strcmp(newname,'delete') + % delete(filename) + % else + % disp(sprintf('New file name is %s ', newname)); + % filename=newname + % end + end + %open file and write description + mode = netcdf.getConstant('NC_NETCDF4'); + mode = bitor(mode,netcdf.getConstant('NC_NOCLOBBER')); %NOCLOBBER to avoid overwrite + ncid = netcdf.create(filename,mode); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description',['Results for run ' md.miscellaneous.name]); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'history',['Created ' datestr(now)]); + + %gather geometry and timestepping as dimensions + if isempty(fieldnames(md.results)), + %results as no field so no time is present + Duration = 0; + else + resfields = fieldnames(md.results); + Duration = size(eval(['md.results. ' resfields{1} ]),2); + end + if Duration>0, + StepNum = Duration; + else + StepNum=1; + end + DimSize(1).index=netcdf.defDim(ncid,'Time',StepNum); %time is the first dimension + [DimSize(1).name,DimSize(1).value]=netcdf.inqDim(ncid,DimSize(1).index); + DimValue(1)=DimSize(1).value; + DimSize(2).index=netcdf.defDim(ncid,'UnLim',netcdf.getConstant('NC_UNLIMITED')); % we add an unlimited dimension if needed + [DimSize(2).name,DimSize(2).value]=netcdf.inqDim(ncid,DimSize(2).index); + DimValue(2)=DimSize(2).value; + % adding mesh related dimensions + dimlist=[2 40 md.mesh.numberofelements md.mesh.numberofvertices size(md.mesh.elements,2)]; + dimnames=["DictDummy" "StringLength" "EltNum" "VertNum" "VertPerElt"]; + if isprop(md.mesh, 'edges'), + dimlist(end+1)=md.mesh.numberofedges; + dimnames(end+1)="EdgeNum"; + else + dimlist(end+1)=0; + dimnames(end+1)="EdgeNum"; + end + if verbose > 0, + disp('===Creating dimensions ==='); + end + %define netcdf dimensions + for i=1:length(dimlist) + % do not add the dimension if it exists already + if sum(dimlist(i) == DimValue) == 0 + DimSize(i+2).index=netcdf.defDim(ncid,dimnames(i),dimlist(i)); + [DimSize(i+2).name,DimSize(i+2).value]=netcdf.inqDim(ncid,DimSize(i+2).index); + DimValue(i+2)=DimSize(i+2).value; + end + end + issmclasses = fieldnames(md)'; + typelist={'half', 'single','double','int8','int16'... + ,'int32','int64','uint8','uint16','uint32'... + ,'uint64','logical','char','string'}; %all malab types that are 0D + + for cl=1:length(issmclasses), + if isempty(md.(issmclasses{cl})), + disp(sprintf("md.%s is empty and will be left as default",issmclasses{cl})); + else + subclasses=fieldnames(md.(issmclasses{cl}))'; + for sc=1:length(subclasses), + if sum(strcmp(class(md.(issmclasses{cl}).(subclasses{sc})), typelist)) == 0, + issmclasses = [issmclasses class(md.(issmclasses{cl}).(subclasses{sc}))]; + end + end + end + end + %get all model classes and create respective groups + groups=fieldnames(md); + if verbose > 0, + disp('===Creating and populating groups==='); + end + for i=1:length(groups), + if verbose >1, + disp(sprintf('===Now treating %s===',groups{i})); + end + if strcmp(groups{i}, 'qmu'), + disp('qmu is skipped until it is more stable'); + continue + end + groupID=netcdf.defGrp(ncid,groups{i}); + %In each group gather the fields of the class + if isempty(md.(groups{i})), + disp(sprintf("WARNING: md.%s is empty, we skip it.",groups{i})) + continue + end + fields=fieldnames(md.(groups{i})); + if isempty(fields), + disp(sprintf("WARNING: md.%s as no fields, we skip it.",groups{i})) + continue + end + %looping on fields in each group + for j=1:length(fields), + Var=md.(groups{i}).(fields{j}); + %treatment for lists + if isa(Var,'cell') + Stdlist=false; %first assume it is not a standard list + if length(Var) == 0 + Stdlist=true; %It is empty and so standard (for us) + else + for k=1:length(typelist) + if isa(Var{1},typelist{k}) + Stdlist=true; %if the list is of a known type (to matlab) if not it is probably some exotic ISSM stuff + end + end + end + %print the issm class as a classtype attribute + klass = class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if(Stdlist) % this is a standard or empty list just proceed + if verbose > 4, + disp(sprintf("=££=creating var for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring)) + end + if ~isempty(Var) && isa(Var{1}, 'char'), % we have a char array, pad it to a given length + Var=char(Var)'; + end + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + else % this is a list of fields, specific treatment needed (perhaps) + if verbose > 4, + disp(sprintf("=??=we have a list of fields for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring)); + end + if strcmp(groups{i}, 'outputdefinition'), + listsize=length(Var); + for k=1:listsize, + subgroupname=md.(groups{i}).(fields{j}){k}.definitionstring; + subgroupID=netcdf.defGrp(groupID,subgroupname); + klass=class(md.(groups{i}).(fields{j}){k}); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(md.(groups{i}).(fields{j}){k}); + for l=1:length(subfields) + if verbose > 4, + disp(sprintf("=--=creating var for %s.%s[%i].%s",groups{i}, fields{j}, k, subfields{l})); + end + Var = md.(groups{i}).(fields{j}){k}.(subfields{l}); + if sum(numel(Var) == size(Var)) == 0, %this is a 2D array or more (and not a vector with dimension 2 = 1) + Var = Var'; + end + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{l},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + end + end + else + disp(sprintf("WARNING: unknown treatment for md.%s",groups{i})); + end + end + elseif sum(strcmp(class(Var), typelist))==1, %this is a standard matlab class with no subgrouping + if verbose > 4, + disp(sprintf("====creating var for %s.%s", groups{i}, fields{j})) + end + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if sum(numel(Var) == size(Var)) == 0, %this is a 2D array or more (and not a vector with dimension 2 = 1) + Var = Var'; + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + + elseif isa(Var,'struct') % structures need special treatment + if strcmp(groups{i}, 'results'), + klasstring=strcat(groups{i} ,'.', groups{i}); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + Listsize= length(md.(groups{i}).(fields{j})); + subgroupname=fields{j}; + subgroupID=netcdf.defGrp(groupID,subgroupname); + klasstring='results.solutionstep'; + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(md.(groups{i}).(fields{j})); + if isempty(subfields), + disp(sprintf("WARNING: md.%s.%s as no subfields, we skip it.",groups{i}, fields{j})); + continue + end + for k=1:length(subfields), + if ~ismember(subfields{k}, {'errlog', 'outlog', 'SolutionType'}) + StackedVar=restable(); + for l=1:Listsize, + Var = md.(groups{i}).(fields{j})(l).(subfields{k}); + if length(Var) == 0, + %Some variables only have data on the first step + break + end + lastindex=l; + StackedVar=StackedVar.update(Var); + end + if verbose > 4, + disp(sprintf("=$$=creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k})); + disp(sprintf("last index on the list is %i",lastindex)); + end + StackedVar=StackedVar.finalize(lastindex); + [DimSize,DimValue,varid]=CreateVar(ncid,StackedVar,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(StackedVar,subgroupID,varid); + end + elseif ismember(subfields{k}, {'SolutionType'}) + %We just add solution type once as an attribute + Var = md.(groups{i}).(fields{j})(1).(subfields{k}); + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + + end + end + elseif strcmp(groups{i}, 'toolkits'), + klasstring=strcat(groups{i} ,'.', groups{i}); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if verbose > 4, + disp(sprintf("=}{=creating var for %s.%s",groups{i}, fields{j})); + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + elseif isempty(fieldnames(md.(groups{i}).(fields{j}))) % this is an empty struct, jus treat it as normal + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if verbose > 4, + disp(sprintf("=[]=creating var for %s.%s",groups{i}, fields{j})); + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + else + disp(sprintf("WARNING, md.%s.%s is not treated as it does not fall in one of the existing cases with class '%s'.",groups{i}, fields{j}, class(md.(groups{i}).(fields{j})))) + end + elseif sum(strcmp(class(Var), issmclasses)) == 1, % that is an issm class + if strcmp(class(Var), 'solution'), + if verbose > 4, + disp(sprintf("=$$=creating var for %s.%s",groups{i}, fields{j})) + disp("NEED treatment") + end + elseif strcmp(class(Var), 'dict'), %we have potential for a dict in py not to sure what it translates to here. + if verbose > 4, + disp(sprintf("=WW=creating var for %s.%s",groups{i}, fields{j})) + disp("NEED Treatment") + end + + else + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subgroupID=netcdf.defGrp(groupID,fields{j}); + klass=class(md.(groups{i}).(fields{j})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(Var); + for k=1:length(subfields), + if sum(strcmp(subfields{k},["outlog" "errlog"])) == 0, + if verbose > 4, + disp(sprintf("+==+creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k})) + end + Var=md.(groups{i}).(fields{j}).(subfields{k}); + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + end + end + + end + else + disp(sprintf("WARNING, md.%s.%s is not treated as it does not fall in one of the existing cases with class '%s'.",groups{i}, fields{j}, class(Var))) + end + end + end + netcdf.close(ncid); +end + +function [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,field,DimSize,DimValue) +% Grab dimensions + varsize=size(Var); + varlength=length(Var); + % treating scalar string or bool as atribute + if isa(Var,'logical'), + if Var, + LogicString='True'; + else, + LogicString='False'; + end + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,LogicString); + varid=[]; + + elseif isa(Var,'char'), + if strcmp(field,'name'), % it looks like netCDF does not like attributes that are called "name" + field = 'varname'; + end + if size(Var,1) <= 1 %that is a single string or empty + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,Var); + varid=[]; + else % that is a character array + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + + elseif isa(Var,'double'), %dealing with arrays + if all(mod(Var, 1) == 0, 'all') %those are actually integers, + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_INT64',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + else + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + elseif isa(Var,'cell'), + % cells can be a range of things, what are we dealing with here + if isempty(Var), + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,'emptycell'); + varid=[]; + else + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + if isa(Var{1}, 'double'), + varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + else + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + end + elseif isa(Var,'struct'), + if isempty(fieldnames(Var)), + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,'emptystruct'); + varid=[]; + else + %Start by getting the structure fields and size + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + else + disp(sprintf('no support for class %s of field %s',class(Var),field)); + varid=[]; + end + return +end + + +function FillVar(Var,groupID,varid) +% Grab dimensions + varsize=size(Var); + varlength=length(Var); + % treating scalar string or bool as atribute + if isa(Var,'double'), %dealing with arrays + if all(mod(Var, 1) == 0, 'all') %those are actually integers, + Var = int64(Var); + end + if length(Var)==0, + netcdf.putVar(groupID,varid,NaN); + else + netcdf.putVar(groupID,varid,Var); + end + elseif isa(Var,'char'), % at this point this should be a character array + netcdf.putVar(groupID,varid,Var); + elseif isa(Var,'cell'), % there can be a number of things in a cell array + for i=1:length(Var), + if isa(Var{i},'char') %for characters we limit the size to 40 for now + if length(Var)>1, + count=[min(length(Var{i}),40), 1]; + startpoint=[0 i-1]; + else + count=min(length(Var{i}),40); + startpoint=0; + end + + if length(Var{i})>40, + netcdf.putVar(groupID,varid,startpoint,count,Var{i}(1:40)); + disp(sprintf('some variable have been truncated')); + else + netcdf.putVar(groupID,varid,startpoint,count,Var{i}); + end + elseif isa(Var{i},'double') + startpoint=[i-1]; + count=[1 length(Var{i}) ndims(Var{i})]; + for j=1:ndims(Var{i}), + startpoint=[startpoint 0]; + end + netcdf.putVar(groupID,varid,startpoint,count,Var{i}); + else + disp(sprintf("WARNING: cell of class %s is not supported.",class(Var{i}))) + end + end + elseif isa(Var,'struct'), + %Start by getting the structure fields and size + locfields=fieldnames(Var); + for i=1:length(locfields), + for j=1:2, + if j==1, + CharVar=locfields{i}'; + disp(size(CharVar)) + if length(CharVar)==0 + CharVar='emptystruct'; + end + startpoint=[0,0,i-1] + else + if isa(Var.(locfields{i}),'char'), + CharVar=Var.(locfields{i})'; + else + CharVar=num2str(Var.(locfields{i}))'; + end + if length(CharVar)==0 + CharVar='emptystruct'; + end + startpoint=[0,1,i-1] + end + + extent=[min(length(CharVar),40), 1, 1] + if length(CharVar)>40, + netcdf.putVar(groupID,varid,startpoint,extent,CharVar(1:40)); + disp(sprintf('some variable have been truncated')); + else + netcdf.putVar(groupID,varid,startpoint,extent,CharVar); + end + end + end + else + disp(sprintf('no support for class %s',class(Var))); + end + return +end + +function [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue) + dims=[]; + celldims=[]; + dim=ndims(Var); + if isa(Var,'struct'), + varsize=length(fieldnames(Var)); + else + varsize=size(Var); + if isa(Var, 'cell') + %we add the dimension of the cells themselves, + %that will most probably fail if cells have different sizes + for i=1:dim, + newdim=size(Var{i}); + if ~ismember(newdim, celldims), + celldims=[celldims newdim]; + end + end + end + end + varsize=[varsize celldims]; + alldim=length(varsize); + if dim>0, + for i=1:alldim, + if size(Var, i)>1 || i>dim || isa(Var, 'struct'), %we skip dimensions with zero lenght but want to add dimensions from cells + indsize=find(varsize(i)==DimValue); + if length(indsize)>0 + dims=[dims DimSize(indsize).index]; + else + indsize=length(DimSize)+1; + DimSize(indsize).index=netcdf.defDim(ncid,['DimNum' num2str(indsize)],varsize(i)); + [DimSize(indsize).name,DimSize(indsize).value]=netcdf.inqDim(ncid,DimSize(indsize).index); + DimValue(indsize)=DimSize(indsize).value; + dims=[dims DimSize(indsize).index]; + end + end + end + end + if isa(Var, 'cell') && isa(Var{1}, 'char'), + %if we have an cell variable with strings we need to add a stringlength + dims=[dims DimSize(4).index]; + end + % struct also need an extra dimension 2, but only if non empty + if isa(Var,'struct'), + dims=[DimSize(4).index DimSize(3).index, dims]; + end +end diff --git a/src/m/netcdf/export_netCDF.py b/src/m/netcdf/OLD/export_netCDF.py similarity index 100% rename from src/m/netcdf/export_netCDF.py rename to src/m/netcdf/OLD/export_netCDF.py diff --git a/src/m/netcdf/OLD/restable.m b/src/m/netcdf/OLD/restable.m new file mode 100644 index 000000000..8fa8ec988 --- /dev/null +++ b/src/m/netcdf/OLD/restable.m @@ -0,0 +1,60 @@ +classdef restable + properties(SetAccess=public) + data=[]; + sizes=[]; + end + methods + function self = update(self, stepvar) + if length(stepvar) == 1 + %if we have a scalar we just add it to the end + %we save the size of the current step for further treatment + self.sizes=[self.sizes;1]; + self.data=[self.data;stepvar]; + else + % if it is an array we add the values one by one + %we save the size of the current step for further treatment + self.sizes=[self.sizes;fliplr(size(stepvar))]; + %we need to transpose to follow the indexing + flatdata=reshape(stepvar', 1, []); + self.data=[self.data,flatdata]; + end + end + function outdat = finalize(self, rows) + if length(self.data)>rows, + if size(self.sizes, 1)==1, + %just one step, data don't need treatment + outdat=self.data; + else, + %we have more scalars than steps, so we have an array + maxsize=[]; + for d=1:size(self.sizes,2) + maxsize=[maxsize,max(self.sizes(:,d))]; + end + findim=[maxsize, rows]; + %first check if all steps are the same size + SameSize = sum(self.sizes - self.sizes(1, :))==0; + if SameSize, + %same size for all steps, just reshape + outdat=reshape(self.data, findim); + else, + %different sizes at each steps, first create a table big enough for the biggest step + startpoint=1; + datadim=ndims(self.data); + outdat=nan(findim); + for r=1:rows + curlen = prod(self.sizes(r, :)); + outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+ ... + curlen-1),self.sizes(r,:)); + startpoint = startpoint+curlen; + end + + end + end, + else, + %as much scalars as steps (or less) so just one value per step + outdat=self.data; + + end + end + end +end \ No newline at end of file diff --git a/src/m/netcdf/README.txt b/src/m/netcdf/README.txt index 606c5a962..928c7c76a 100644 --- a/src/m/netcdf/README.txt +++ b/src/m/netcdf/README.txt @@ -1,87 +1,114 @@ -The write_netCDF and read_netCDF modules provide a convenient way to save and restore the state of a model class instance -in binary format via NetCDF4. This allows users to store the class state on disk and retrieve it later, facilitating seamless -transitions between Python and MATLAB environments. - -To save a model, call either write_netCDF.py or write_netCDF.m depending on whether your class is in matlab or python. -To read a saved model, call either read_netCDF.py or read_netCDF.m depending on what language you prefer to use the model in. -If you would like to log the names and locations of variables being stored, add the argument verbose = True (verbose = true for matlab). - -Usage Instructions: - - Python: - - Saving a model: - from write_netCDF import write_netCDF - - md = bamg(model(), foo.csv, .01) - - write_netCDF(md, 'adress_to_save/../filename.nc') - - - Reading a model: - from read_netCDF import read_netCDF - - md = read_netCDF('adress_to_file/../filename.nc') - - Verbose examples: - write_netCDF(md, adress_to_save/../filename.nc, verbose = True) - md = read_netCDF(adress_to_file/../filename.nc, verbose = True) - - MATLAB: - - Saving a model: - - write_netCDF(md, adress_to_save/../filename.nc); - - - Reading a model: - - md = read_netCDF(adress_to_file/../filename.nc); - - Verbose examples: - write_netCDF(md, adress_to_save/../filename.nc, verbose = true); - - or: - - write_netCDF(md, adress_to_save/../filename.nc, verbose); - md = read_netCDF(adress_to_file/../filename.nc, verbose = true); - -Dependencies: - Python: - - NumPy - - NetCDF4 / NetCDF4.Dataset - - The model() class - - results.solution / results.solutionstep / results.resultsdakota - - inversion.inversion / inversion.m1qn3inversion / inversion.taoinversion - - MATLAB: - - The model() class - - inversion.inversion / inversion.m1qn3inversion / inversion.taoinversion - - -Additional Information: - -There are currently datatypes that both write_netCDF and read_netCDF modules may not be able to handle. These datatypes might -include lists with multiple datatypes (ie, ['number', 1, 'letter', a, 'color', 'blue']), lists of dicts ect. - -To add functionality for these additional cases, one must simply create a function to handle the case and call it using a -conditional case within the create_var() function. To read the data from the NetCDF4 file, add the case to the -copy_variable_data_to_new_model() function in read_netCDF so that the data can be added to a new model() instance. - -Known issues: - -Unlike Python, MATLAB doesn't utilize subclasses in its model class. This leads to a loss of certain subclass instances. -For instance, the results.solutionstep() class poses a known issue. In MATLAB, there's no direct equivalent. The fields in -'md.results' in MATLAB might correspond to instances of resultsdakota(), solution(), or solutionstep() in Python, but -because those classes don't exist in MATLAB, there is no way for python to know which instance it needs. - -The current workaround, while not theoretically sound, involves searching for the class name string in MATLAB's 'results' -field names. For instance, 'md.results.TransientSolution' is recorded as a solution() class instance. However, problems arise -in cases like 'md.results.StressbalanceSolution', where the code notes a solution() instance, while in Python, it should be a -solutionstep() instance. - -So far, there have been no recorded problems swapping a solutionstep() instance for a solution() instance. - -Potential solutions are: - - - Restructure both Python and MATLAB solve frameworks. In Python, when creating an md.results. instance, - embed 'solutionstep' in the class instance name. - >> This solution is very involved, and would include the tedious modification of >5 files in total - - Create a hash table linking solutions with their corresponding 'md.results.' for reference when saving models to - the netCDF file. +write_netCDF / read_netCDF – Cross-language ISSM model serialisation +====================================================================== + +These four files let you save an ISSM model (md) to NetCDF4 and reload it +in either MATLAB or Python, interchangeably. + + write_netCDF.m – save from MATLAB + read_netCDF.m – load into MATLAB + write_netCDF.py – save from Python + read_netCDF.py – load into Python + +Version: 2.0 + + +USAGE +----- + +Python: + from write_netCDF import write_netCDF + from read_netCDF import read_netCDF + + write_netCDF(md, 'model.nc') + write_netCDF(md, 'model.nc', verbose=True) + + md2 = read_netCDF('model.nc') + md2 = read_netCDF('model.nc', verbose=True) + +MATLAB: + write_netCDF(md, 'model.nc') + write_netCDF(md, 'model.nc', 'verbose', true) + + md2 = read_netCDF('model.nc') + md2 = read_netCDF('model.nc', 'verbose', true) + +Cross-language round-trip: + %% MATLAB → Python + write_netCDF(md, 'model.nc'); % MATLAB + md2 = read_netCDF('model.nc') # Python + + ## Python → MATLAB + write_netCDF(md, 'model.nc') # Python + md2 = read_netCDF('model.nc'); % MATLAB + + +DEPENDENCIES +------------ +Python: numpy, netCDF4, and the ISSM model() class (plus any sub-class + modules that your model uses, e.g. m1qn3inversion, SMBpdd, etc.) +MATLAB: built-in netcdf package (no extra toolboxes required), plus + the ISSM model() class and any sub-class .m files. + + +FILE FORMAT (ISSM-NetCDF-2.0) +------------------------------ +Global attributes: + Conventions = 'ISSM-NetCDF-2.0' + history = creation timestamp + +Layout: + /mesh/ ← group per top-level md field + classtype ← NC_GLOBAL attribute: Python/MATLAB class name + x ← variable + y ← variable + … + /geometry/ + … + /results/ + /TransientSolution/ ← sub-group (classtype='struct') + nsteps ← NC_GLOBAL attribute: number of steps + /step_1/ + Vel ← variable + … + /step_2/ + … + /inversion/ + classtype = 'm1qn3inversion' ← enables correct class reconstruction + … + +Key design decisions: + - Every group carries a 'classtype' global attribute so readers can + reconstruct the exact Python/MATLAB class without heuristics. + - Struct arrays (results.TransientSolution) are stored as step_1 … step_n + sub-groups, one per time step. + - The file is overwritten silently if it already exists (no interactive + prompts – safe for use in scripts and Jupyter notebooks). + - No module-level global state in Python (fully re-entrant; calling + read_netCDF twice returns two independent model objects). + - No MATLAB 'persistent' variable bugs (the old code skipped writing + inversion/smb/friction/hydrology class names on the second call). + + +KNOWN LIMITATIONS +----------------- +- md.qmu is not yet supported (complex OrderedDict-based structure). +- Cell arrays whose elements are themselves struct arrays are not yet + supported. +- MATLAB classes that have no Python equivalent (e.g. SMBgemb) can be + saved from MATLAB and round-tripped back to MATLAB, but the Python + reader will fall back to a plain dict for those fields. +- Very large arrays (> a few GB) work fine thanks to NetCDF4 chunking, + but you may want to increase the zlib compression level in write_netCDF.py + for highly compressible fields. + + +OTHER FILES IN THIS DIRECTORY +------------------------------ +export_netCDF.m / export_netCDF.py + An older, alternative serialisation that stacks results variables by + time dimension rather than using per-step sub-groups. It writes a + different file format and does not have a matching read_netCDF. + Kept for backward compatibility. + +restable.m + Helper used by export_netCDF.m. diff --git a/src/m/netcdf/read_netCDF.m b/src/m/netcdf/read_netCDF.m index 81cca0655..2d85a3d49 100644 --- a/src/m/netcdf/read_netCDF.m +++ b/src/m/netcdf/read_netCDF.m @@ -1,526 +1,478 @@ -%{ -Given a NetCDF4 file, this set of functions will perform the following: - 1. Enter each group of the file. - 2. For each variable in each group, update an empty model with the variable's data - 3. Enter nested groups and repeat - - -If the model you saved has subclass instances that are not in the standard model() class -you can: - 1. Copy lines 30-35, set the "results" string to the name of the subclass instance, - 2. Copy and modify the make_results_subclasses() function to create the new subclass - instances you need. -From there, the rest of this script will automatically create the new subclass -instance in the model you're writing to and store the data from the netcdf file there. -%} - - -function model_copy = read_netCDF(filename, varargin) - if nargin > 1 - verbose = true; +function md = read_netCDF(filename, varargin) +% READ_NETCDF - Load an ISSM model from a NetCDF4 file written by write_netCDF.m or write_netCDF.py. +% +% Usage: +% md = read_netCDF(filename) +% md = read_netCDF(filename, 'verbose', true) +% +% Inputs: +% filename - path to the .nc file +% +% Optional name-value pair: +% 'verbose' - true/false (default false) +% +% Returns: +% md - an ISSM model() object populated from the file + + % --- parse options ------------------------------------------------ + if nargin > 1 && islogical(varargin{1}) + verbose = varargin{1}; + elseif nargin > 1 + p = inputParser(); + p.addParameter('verbose', false, @islogical); + p.parse(varargin{:}); + verbose = p.Results.verbose; else verbose = false; end - + if verbose - fprintf('NetCDF42C v1.1.14\n'); + fprintf('read_netCDF v2.0 (MATLAB)\n'); end - % make a model framework to fill that is in the scope of this file - model_copy = model(); - - % Check if path exists - if exist(filename, 'file') - if verbose - fprintf('Opening %s for reading\n', filename); - end - - % Open the given netCDF4 file - NCData = netcdf.open(filename, 'NOWRITE'); - % Remove masks from netCDF data for easy conversion: NOT WORKING - %netcdf.setMask(NCData, 'NC_NOFILL'); - - % see if results is in there, if it is we have to instantiate some classes - try - results_group_id = netcdf.inqNcid(NCData, "results"); - model_copy = make_results_subclasses(model_copy, NCData, verbose); - catch - end % 'results' group doesn't exist - % see if inversion is in there, if it is we may have to instantiate some classes - try - inversion_group_id = netcdf.inqNcid(NCData, "inversion"); - model_copy = check_inversion_class(model_copy, NCData, verbose); - catch - end % 'inversion' group doesn't exist - - % loop over first layer of groups in netcdf file - for group = netcdf.inqGrps(NCData) - group_id = netcdf.inqNcid(NCData, netcdf.inqGrpName(group)); - %disp(netcdf.inqGrpNameFull(group_id)) - % hand off first level to recursive search - model_copy = walk_nested_groups(group_id, model_copy, NCData, verbose); - end - - % Close the netCDF file - netcdf.close(NCData); - if verbose - disp('Model Successfully Copied') - end - else - fprintf('File %s does not exist.\n', filename); + if ~exist(filename, 'file') + error('read_netCDF: file not found: %s', filename); end -end + % Start with a fully-initialised model + md = model(); -function model_copy = make_results_subclasses(model_copy, NCData, verbose) - resultsGroup = netcdf.inqNcid(NCData, "results"); - variables = netcdf.inqVarIDs(resultsGroup); - for name = variables - class_instance = netcdf.inqVar(resultsGroup, name); - class_instance_names_raw = netcdf.getVar(resultsGroup, name, 'char').'; - class_instance_names = cellstr(class_instance_names_raw); - for index = 1:numel(class_instance_names) - class_instance_name = class_instance_names{index}; - model_copy.results = setfield(model_copy.results, class_instance_name, struct()); - end - %model_copy.results = setfield(model_copy.results, class_instance, class_instance_name); + ncid = netcdf.open(filename, 'NOWRITE'); + try + md = read_all_groups(ncid, md, verbose); + catch ME + netcdf.close(ncid); + rethrow(ME); end - model_copy = model_copy; + netcdf.close(ncid); + if verbose - disp('Successfully recreated results structs:') - for fieldname = string(fieldnames(model_copy.results)) - disp(fieldname) - end + disp('Model successfully loaded from NetCDF4.'); end end +function md = read_all_groups(ncid, md, verbose) % {{{ +% ===================================================================== +% read_all_groups – iterate over every top-level group +% ===================================================================== + top_groups = netcdf.inqGrps(ncid); + for gi = 1:numel(top_groups) + gid = top_groups(gi); + gname = netcdf.inqGrpName(gid); -function model_copy = check_inversion_class(model_copy, NCData, verbose) - % get the name of the inversion class: either inversion or m1qn3inversion or taoinversion - inversionGroup = netcdf.inqNcid(NCData, "inversion"); - varid = netcdf.inqVarID(inversionGroup, 'inversion_class_name'); - inversion_class = convertCharsToStrings(netcdf.getVar(inversionGroup, varid,'char')); - if strcmp(inversion_class, 'm1qn3inversion') - model_copy.inversion = m1qn3inversion(); - if verbose - disp('Successfully created inversion class instance: m1qn3inversion') - end - elseif strcmp(inversion_class, 'taoinversion') - model_copy.inversion = taoinversion(); if verbose - disp('Successfully created inversion class instance: taoinversion') + fprintf(' Reading group: %s\n', gname); end - else - if verbose - disp('No inversion class was found') - end - end - model_copy = model_copy; -end + % Instantiate the right subclass based on the stored classtype attribute + md = instantiate_subclass(md, gname, gid, verbose); -function model_copy = walk_nested_groups(group_location_in_file, model_copy, NCData, verbose) - % we search the current group level for variables by getting this struct - variables = netcdf.inqVarIDs(group_location_in_file); - - % from the variables struct get the info related to the variables - for variable = variables - [varname, xtype, dimids, numatts] = netcdf.inqVar(group_location_in_file, variable); - - % keep an eye out for nested structs: - if strcmp(varname, 'this_is_a_nested') - is_object = true; - model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose); - elseif strcmp(varname, 'name_of_cell_array') - is_object = true; - model_copy = copy_cell_array_of_objects(variables, group_location_in_file, model_copy, NCData, verbose); - elseif strcmp(varname, 'solution') - % band-aid pass.. - else - if logical(exist('is_object', 'var')) - % already handled - else - model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose); - end - end + % Now populate all fields + md = read_group_into_model(gid, md, gname, verbose); + end +end % }}} +function md = instantiate_subclass(md, gname, gid, verbose) % {{{ +% ===================================================================== +% instantiate_subclass – create the right class for polymorphic fields +% ===================================================================== + % Read the classtype attribute from this group (written by write_netCDF) + ct = get_group_classtype(gid); + if isempty(ct) + return % no classtype stored, keep existing model default end - % try to find groups in current level, if it doesn't work it's because there is nothing there - %try - % if it's a nested struct the function copy_nested_struct has already been called - if logical(exist('is_object', 'var')) - % do nothing - else - % search for nested groups in the current level to feed back to this function - groups = netcdf.inqGrps(group_location_in_file); - if not(isempty(groups)) - for group = groups - group_id = netcdf.inqNcid(group_location_in_file, netcdf.inqGrpName(group)); - %disp(netcdf.inqGrpNameFull(group_id)) - model_copy = walk_nested_groups(group, model_copy, NCData, verbose); + % Only act for fields whose type can vary at runtime + switch gname + case 'inversion' + switch ct + case 'm1qn3inversion' + md.inversion = m1qn3inversion(); + case 'taoinversion' + md.inversion = taoinversion(); + % 'inversion' is already the default end - end + case 'smb' + md = instantiate_smb(md, ct, verbose); + case 'friction' + md = instantiate_friction(md, ct, verbose); + case 'hydrology' + md = instantiate_hydrology(md, ct, verbose); + case 'mesh' + switch ct + case 'mesh3dprisms' + md.mesh = mesh3dprisms(); + % mesh2d is the default + end + % All other groups keep their default instance; we just overwrite fields + end + if verbose && ~isempty(ct) + fprintf(' classtype = %s\n', ct); end - %catch % no nested groups here - %end end +% }}} +function md = read_group_into_model(gid, md, gname, verbose) % {{{ +% ===================================================================== +% read_group_into_model – populate md. from a NetCDF group +% ===================================================================== + % Skip empty/missing model fields gracefully + try + target = md.(gname); + catch + if verbose + fprintf(' [SKIP] md.%s does not exist in model()\n', gname); + end + return + end + % Read variables at this level + var_ids = netcdf.inqVarIDs(gid); + for vi = 1:numel(var_ids) + varid = var_ids(vi); + vname = netcdf.inqVar(gid, varid); + data = read_variable(gid, varid, verbose); + md.(gname) = set_field(md.(gname), vname, data, verbose); + end -% to read cell arrays with objects: -function model_copy = copy_cell_array_of_objects(variables, group_location_in_file, model_copy, NCData, verbose); - %{ - The structure in netcdf for groups with the name_of_cell_array variable is like: - - group: 2x6_cell_array_of_objects { - name_of_cell_array = - - group: Row_1_of_2 { - group: Col_1_of_6 { - ... other groups can be here that refer to objects - } // group Col_6_of_6 - } // group Row_1_of_2 - - group: Row_2_of_2 { - group: Col_1_of_6 { - ... other groups can be here that refer to objects - } // group Col_6_of_6 - } // group Row_2_of_2 - } // group 2x6_cell_array_of_objects - - We have to navigate this structure to extract all the data and recreate the - original structure when the model was saved - %} - - % get the name_of_cell_array, rows and cols vars - name_of_cell_array_varID = netcdf.inqVarID(group_location_in_file, 'name_of_cell_array'); - rows_varID = netcdf.inqVarID(group_location_in_file, 'rows'); - cols_varID = netcdf.inqVarID(group_location_in_file, 'cols'); - - name_of_cell_array = netcdf.getVar(group_location_in_file, name_of_cell_array_varID).'; % transpose - rows = netcdf.getVar(group_location_in_file, rows_varID); - cols = netcdf.getVar(group_location_in_file, cols_varID); - - % now we work backwards: make the cell array, fill it in, and assign it to the model - - % make the cell array - cell_array_placeholder = cell(rows, cols); - - % get subgroups which are elements of the cell array - subgroups = netcdf.inqGrps(group_location_in_file); % numerical cell array with ID's of subgroups - - % enter each subgroup, get the data, assign it to the corresponding index of cell array - if rows > 1 - % we go over rows - % set index for cell array rows - row_idx = 1; - for row = subgroups - % now columns - columns = netcdf.inqGrps(group_location_in_file); - - % set index for cell array cols - col_idx = 1; - for column = columns - % now variables - current_column_varids = netcdf.inqVarIDs(column); - - % if 'class_is_a' or 'this_is_a_nested' variables is present at this level we have to handle them accordingly - try - class_is_aID = netcdf.inqVarID(column, 'class_is_a'); - col_data = deserialize_class(column, NCData, verbose); - is_object = true; - catch - end - - try - this_is_a_nestedID = netcdf.inqVarID(column, 'this_is_a_nested'); - % functionality not supported - disp('Error: Cell Arrays of structs not yet supported!') - % copy_nested_struct(column, model_copy, NCData, verbose) - is_object = true; - catch - end + % Recurse into subgroups + sub_groups = netcdf.inqGrps(gid); + for si = 1:numel(sub_groups) + sub_gid = sub_groups(si); + sub_name = netcdf.inqGrpName(sub_gid); + ct = get_group_classtype(sub_gid); - if logical(exist('is_object', 'var')) - % already taken care of - else - % store the variables as normal -- to be added later - disp('Error: Cell Arrays of mixed objects not yet supported!') - for var = current_column_varids - % not supported - end - end + if strcmp(ct, 'struct') + % This is a struct array (e.g. results.TransientSolution) + n = get_group_int_att(sub_gid, 'nsteps', 1); + s = read_struct_array(sub_gid, n, verbose); + md.(gname) = set_field(md.(gname), sub_name, s, verbose); - cell_array_placeholder{row_idx, col_idx} = col_data; - col_idx = col_idx + 1; - end - row_idx = row_idx + 1; - end - else - % set index for cell array - col_idx = 1; - for column = subgroups - % now variables - current_column_varids = netcdf.inqVarIDs(column); + elseif strcmp(ct, 'cell_of_objects') + % Cell array of ISSM class instances + c = read_cell_of_objects(sub_gid, verbose); + md.(gname) = set_field(md.(gname), sub_name, c, verbose); - % if 'class_is_a' or 'this_is_a_nested' variables is present at this level we have to handle them accordingly + else + % Regular sub-class: read its variables into the existing sub-object try - classID = netcdf.inqVarID(column, 'class_is_a'); - col_data = deserialize_class(classID, column, NCData, verbose); - is_object = true; - catch ME - rethrow(ME) + sub_obj = md.(gname).(sub_name); + catch + % Field doesn't exist in the default model – create a struct + sub_obj = struct(); end - + + sub_obj = read_subgroup_into_obj(sub_gid, sub_obj, verbose); + md.(gname) = set_field(md.(gname), sub_name, sub_obj, verbose); + end + end +end +% }}} +function obj = read_subgroup_into_obj(gid, obj, verbose) % {{{ +% ===================================================================== +% read_subgroup_into_obj – generic recursive group reader +% ===================================================================== + % Variables + var_ids = netcdf.inqVarIDs(gid); + for vi = 1:numel(var_ids) + varid = var_ids(vi); + vname = netcdf.inqVar(gid, varid); + data = read_variable(gid, varid, verbose); + obj = set_field(obj, vname, data, verbose); + end + + % Child subgroups + sub_groups = netcdf.inqGrps(gid); + for si = 1:numel(sub_groups) + sub_gid = sub_groups(si); + sub_name = netcdf.inqGrpName(sub_gid); + ct = get_group_classtype(sub_gid); + + if strcmp(ct, 'struct') + n = get_group_int_att(sub_gid, 'nsteps', 1); + s = read_struct_array(sub_gid, n, verbose); + obj = set_field(obj, sub_name, s, verbose); + elseif strcmp(ct, 'cell_of_objects') + c = read_cell_of_objects(sub_gid, verbose); + obj = set_field(obj, sub_name, c, verbose); + else + % Try to get existing sub-object, fall back to empty struct try - this_is_a_nestedID = netcdf.inqVarID(column, 'this_is_a_nested'); - % functionality not supported - disp('Error: Cell Arrays of structs not yet supported!') - % col_data = copy_nested_struct(column, model_copy, NCData, verbose); - is_object = true; + sub_obj = obj.(sub_name); catch + sub_obj = struct(); end - if logical(exist('is_object', 'var')) - % already taken care of - else - % store the variables as normal -- to be added later - disp('Error: Cell Arrays of mixed objects not yet supported!') - for var = current_column_varids - % col_data = not supported - end - end - - cell_array_placeholder{col_idx} = col_data; - col_idx = col_idx + 1; + sub_obj = read_subgroup_into_obj(sub_gid, sub_obj, verbose); + obj = set_field(obj, sub_name, sub_obj, verbose); + end + end +end % }}} +function s = read_struct_array(gid, n, verbose) % {{{ +% ===================================================================== +% read_struct_array – reconstruct a 1×n struct array from step_k groups +% ===================================================================== + % Gather all step subgroups (step_1, step_2, …) + step_groups = netcdf.inqGrps(gid); + if isempty(step_groups) + s = struct(); + return + end - end + % Collect all field names from the first step + first_gid = step_groups(1); + first_varids = netcdf.inqVarIDs(first_gid); + fnames = {}; + for vi = 1:numel(first_varids) + fnames{end+1} = netcdf.inqVar(first_gid, first_varids(vi)); %#ok end - - % Like in copy_nested_struct, we can only handle things 1 layer deep. - % assign cell array to model - address_to_attr_list = split(netcdf.inqGrpNameFull(group_location_in_file), '/'); - address_to_attr = address_to_attr_list{2}; - if isprop(model_copy.(address_to_attr), name_of_cell_array); - model_copy.(address_to_attr).(name_of_cell_array) = cell_array_placeholder; + % Pre-allocate struct array + if ~isempty(fnames) + args = [fnames; repmat({[]}, 1, numel(fnames))]; + s = repmat(struct(args{:}), 1, numel(step_groups)); else - model_copy = addprop(model_copy.(address_to_attr), name_of_cell_array, cell_array_placeholder); + s = repmat(struct(), 1, numel(step_groups)); end + for ki = 1:numel(step_groups) + step_gid = step_groups(ki); + var_ids = netcdf.inqVarIDs(step_gid); + for vi = 1:numel(var_ids) + varid = var_ids(vi); + vname = netcdf.inqVar(step_gid, varid); + data = read_variable(step_gid, varid, verbose); + s(ki).(vname) = data; + end + % Recurse into any sub-subgroups (rare but possible) + sub_grps = netcdf.inqGrps(step_gid); + for si = 1:numel(sub_grps) + sg_name = netcdf.inqGrpName(sub_grps(si)); + sub_obj = read_subgroup_into_obj(sub_grps(si), struct(), verbose); + s(ki).(sg_name) = sub_obj; + end + end if verbose - fprintf("Successfully loaded cell array %s to %s\n", name_of_cell_array,address_to_attr_list{2}) + fprintf(' [struct] loaded 1x%d struct array\n', numel(step_groups)); end -end - - - - -function output = deserialize_class(classID, group, NCData, verbose) - %{ - This function will recreate a class - %} - - % get the name of the class - name = netcdf.getVar(group, classID).'; +end % }}} +function c = read_cell_of_objects(gid, verbose) % {{{ +% ===================================================================== +% read_cell_of_objects – reconstruct a cell array of ISSM objects +% ===================================================================== + nrows = get_group_int_att(gid, 'nrows', 1); + ncols = get_group_int_att(gid, 'ncols', 1); + c = cell(nrows, ncols); + + item_groups = netcdf.inqGrps(gid); + for gi = 1:numel(item_groups) + ig = item_groups(gi); + name = netcdf.inqGrpName(ig); % e.g. 'item_1_3' + ct = get_group_classtype(ig); + + % Parse row/col from name + tok = regexp(name, '^item_(\d+)_(\d+)$', 'tokens'); + if isempty(tok) + continue + end + r = str2double(tok{1}{1}); + col_idx = str2double(tok{1}{2}); - % instantiate it - class_instance = eval([name, '()']); + % Instantiate the object + try + obj = eval([ct '()']); + catch + obj = struct(); + end + obj = read_subgroup_into_obj(ig, obj, verbose); + c{r, col_idx} = obj; + end + if verbose + fprintf(' [cell] loaded %dx%d cell of objects\n', nrows, ncols); + end +end % }}} +function data = read_variable(gid, varid, verbose) % {{{ +% ===================================================================== +% read_variable – read one NetCDF variable and convert to MATLAB type +% ===================================================================== + [~, xtype, dimids, ~] = netcdf.inqVar(gid, varid); + + % Get type_is attribute if present + type_is = ''; + try + type_is = netcdf.getAtt(gid, varid, 'type_is'); + catch + end - % get and assign properties - subgroups = netcdf.inqGrps(group); % numerical cell array with ID's of subgroups + raw = netcdf.getVar(gid, varid); - if numel(subgroups) == 1 - % get properties - varIDs = netcdf.inqVarIDs(subgroups); - for varID = varIDs - % var metadata - [varname, xtype, dimids, numatts] = netcdf.inqVar(subgroups, varID); - % data - data = netcdf.getVar(subgroups, varID); + % Empty / NaN sentinel + if strcmp(type_is, 'empty') + data = []; + return + end - % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim - if all(size(data)~=1) || xtype == 2 - data = data.'; - end + % Bool + if strcmp(type_is, 'bool') + data = logical(raw); + return + end - % some classes have permissions... so we skip those - try - % if property already exists, assign new value - if isprop(class_instance, varname) - class_instance.(varname) = data; - else - addprop(class_instance, varname, data); - end - catch - end + % String + if strcmp(type_is, 'string') || xtype == 2 % NC_CHAR = 2 + if strcmp(type_is, 'cell_of_strings') + % char matrix → cell array of strings (trim trailing spaces) + data = cellstr(raw.'); + else + data = char(raw.'); end - else - % not supported + return end - output = class_instance; -end + % Cell of strings (attribute check) + if strcmp(type_is, 'cell_of_strings') + data = cellstr(raw.'); + return + end -function model_copy = copy_nested_struct(group_location_in_file, model_copy, NCData, verbose) - %{ - A common multidimensional struct array is the 1xn md.results.TransientSolution struct. - The process to recreate is as follows: - 1. Get the name of the struct from group name - 2. Get the fieldnames from the subgroups - 3. Recreate the struct with fieldnames - 4. Populate the fields with their respective values - %} - - % step 1 - name_of_struct = netcdf.inqGrpName(group_location_in_file); - - % step 2 - subgroups = netcdf.inqGrps(group_location_in_file); % numerical cell array with ID's of subgroups - % get single subgroup's data - single_subgroup_ID = subgroups(1); - subgroup_varids = netcdf.inqVarIDs(single_subgroup_ID); - fieldnames = {}; - for variable = subgroup_varids - [varname, xtype, dimids, numatts] = netcdf.inqVar(single_subgroup_ID, variable); - fieldnames{end+1} = varname; - end - - % step 3 - address_in_model_raw = split(netcdf.inqGrpNameFull(group_location_in_file), '/'); - address_in_model = address_in_model_raw{2}; - - % we cannot assign a variable to represent this object as MATLAB treats all variables as copies - % and not pointers to the same memory address - % this means that if address_in_model has more than 1 layer, we need to modify the code. For now, - % we just hope this will do. An example of a no-solution would be model().abc.def.ghi.field whereas we're only assuming model().abc.field now - - model_copy.(address_in_model).(name_of_struct) = struct(); - % for every fieldname in the subgroup, create an empty field - for fieldname = string(fieldnames) - model_copy.(address_in_model).(name_of_struct).(fieldname) = {}; - end - - % use repmat to make the struct array multidimensional along the fields axis - number_of_dimensions = numel(subgroups); - model_copy.(address_in_model).(name_of_struct) = repmat(model_copy.(address_in_model).(name_of_struct), 1, number_of_dimensions); - - % step 4 - % for every layer of the multidimensional struct array, populate the fields - for current_layer = 1:number_of_dimensions - % choose subgroup - current_layer_subgroup_ID = subgroups(current_layer); - % get all vars - current_layer_subgroup_varids = netcdf.inqVarIDs(current_layer_subgroup_ID); - % get individual vars and set fields at layer current_layer - for varid = current_layer_subgroup_varids - [varname, xtype, dimids, numatts] = netcdf.inqVar(current_layer_subgroup_ID, varid); - data = netcdf.getVar(current_layer_subgroup_ID, varid); - - % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim - if all(size(data)~=1) || xtype == 2 - data = data.'; + % Numeric: NetCDF is row-major, MATLAB is column-major → transpose 2-D arrays + data = double(raw); + if ndims(data) == 2 && ~any(size(data) == 1) + data = data.'; + elseif iscolumn(data) + % keep as column vector (MATLAB convention) + end +end % }}} +function obj = set_field(obj, fname, data, verbose) % {{{ +% ===================================================================== +% set_field – safely set a field on either a class or struct +% ===================================================================== + try + if isstruct(obj) + obj.(fname) = data; + elseif isobject(obj) + if isprop(obj, fname) || isfield(obj, fname) + obj.(fname) = data; + else + % Property not in the default class – try dynamic property + try + obj.(fname) = data; + catch + if verbose + fprintf(' [SKIP] cannot set property %s on %s\n', fname, class(obj)); + end + end end - - % set the field - model_copy.(address_in_model).(name_of_struct)(current_layer).(varname) = data; - %address_to_struct_in_model = setfield(address_to_struct_in_model(current_layer), varname, data) end - model_copy.(address_in_model).(name_of_struct)(current_layer); + catch ME if verbose - fprintf("Successfully loaded layer %s to multidimension struct array\n", num2str(current_layer)) + fprintf(' [WARN] set_field failed for %s: %s\n', fname, ME.message); end end - model_copy = model_copy; - if verbose - fprintf('Successfully recreated multidimensional structure array %s in md.%s\n', name_of_struct, address_in_model) - end end - - - - -%{ -Since there are two types of objects that MATLAB uses (classes and structs), we have to check -which object we're working with before we can set any fields/attributes of it. After this is completed, -we can write the data to that location in the model. -%} - -function model_copy = copy_variable_data_to_new_model(group_location_in_file, varname, xtype, model_copy, NCData, verbose) - %disp(varname) - % this is an inversion band-aid - if strcmp(varname, 'inversion_class_name') || strcmp(varname, 'name_of_struct') || strcmp(varname, 'solution') - % we don't need this - else - % putting try/catch here so that any errors generated while copying data are logged and not lost by the try/catch in walk_nested_groups function - try - %disp(netcdf.inqGrpNameFull(group_location_in_file)) - %disp(class(netcdf.inqGrpNameFull(group_location_in_file))) - address_to_attr = strrep(netcdf.inqGrpNameFull(group_location_in_file), '/', '.'); - varid = netcdf.inqVarID(group_location_in_file, varname); - data = netcdf.getVar(group_location_in_file, varid); - - - % if we have an empty string - if xtype == 2 && isempty(all(data)) - data = cell(char()); - % if we have an empty cell-char array - elseif numel(data) == 1 && xtype == 3 && data == -32767 - data = cell(char()); - elseif isempty(all(data)) - data = [] +% }}} +function md = instantiate_smb(md, ct, verbose) % {{{ +% ===================================================================== +% Subclass instantiation helpers +% ===================================================================== + switch ct + case 'SMBforcing', md.smb = SMBforcing(); + case 'SMBpdd', md.smb = SMBpdd(); + case 'SMBd18opdd', md.smb = SMBd18opdd(); + case 'SMBgradients', md.smb = SMBgradients(); + case 'SMBcomponents', md.smb = SMBcomponents(); + case 'SMBmeltcomponents', md.smb = SMBmeltcomponents(); + case 'SMBgradientscomponents' + if exist('SMBgradientscomponents','class') + md.smb = SMBgradientscomponents(); end - % band-aid for some cell-char-arrays: - if xtype == 2 && strcmp(data, 'default') - data = {'default'}; + case 'SMBgradientsela' + if exist('SMBgradientsela','class') + md.smb = SMBgradientsela(); end - - % netcdf uses Row Major Order but MATLAB uses Column Major Order so we need to transpose all arrays w/ more than 1 dim - if all(size(data)~=1) || xtype == 2 - data = data.'; + case 'SMBhenning' + if exist('SMBhenning','class') + md.smb = SMBhenning(); end - - % if we have a list of strings - if xtype == 2 - try - if strcmp(netcdf.getAtt(group_location_in_file, varid, "type_is"), 'cell_array_of_strings') - data = cellstr(data); - end - catch - % no attr found so we pass - end + case 'SMBgemb' + if exist('SMBgemb','class') + md.smb = SMBgemb(); end - - % the issm c compiler does not work with int64 datatypes, so we need to convert those to int16 - % reference this (very hard to find) link for netcdf4 datatypes: https://docs.unidata.ucar.edu/netcdf-c/current/netcdf_8h_source.html - %xtype - if xtype == 10 - arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = ' , 'double(data);']; - eval(arg_to_eval); - %disp('Loaded int64 as int16') - else - arg_to_eval = ['model_copy', address_to_attr, '.', varname, ' = data;']; - eval(arg_to_eval); + case 'SMBpddSicopolis' + if exist('SMBpddSicopolis','class') + md.smb = SMBpddSicopolis(); end - + case 'SMBsemic' + if exist('SMBsemic','class') + md.smb = SMBsemic(); + end + case 'SMBdebrisEvatt' + if exist('SMBdebrisEvatt','class') + md.smb = SMBdebrisEvatt(); + end + otherwise if verbose - full_addy = netcdf.inqGrpNameFull(group_location_in_file); - %disp(xtype) - %class(data) - fprintf('Successfully loaded %s to %s\n', varname, full_addy); + fprintf(' [WARN] unknown SMB class: %s\n', ct); end - - catch ME %ME is an MException struct - % Some error occurred if you get here. - fprintf(1,'There was an error with %s! \n', varname) - errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ME.stack.name, ME.stack.line, ME.message); - fprintf(1, '%s\n', errorMessage); - uiwait(warndlg(errorMessage)); - %line = ME.stack.line - %fprintf(1,'There was an error with %s! \n', varname) - %fprintf('The message was:\n%s\n',ME.message); - %fprintf(1,'The identifier was:\n%s\n',ME.identifier); - - % more error handling... - end end - model_copy = model_copy; +end % }}} +function md = instantiate_friction(md, ct, verbose) % {{{ + switch ct + case 'friction', md.friction = friction(); + case 'frictioncoulomb', md.friction = frictioncoulomb(); + case 'frictioncoulomb2', md.friction = frictioncoulomb2(); + case 'frictionhydro', md.friction = frictionhydro(); + case 'frictionjosh', md.friction = frictionjosh(); + case 'frictionpism', md.friction = frictionpism(); + case 'frictionregcoulomb', md.friction = frictionregcoulomb(); + case 'frictionregcoulomb2', md.friction = frictionregcoulomb2(); + case 'frictionschoof', md.friction = frictionschoof(); + case 'frictionshakti', md.friction = frictionshakti(); + case 'frictiontsai', md.friction = frictiontsai(); + case 'frictionwaterlayer', md.friction = frictionwaterlayer(); + case 'frictionweertman', md.friction = frictionweertman(); + case 'frictionweertmantemp', md.friction = frictionweertmantemp(); + otherwise + if verbose + fprintf(' [WARN] unknown friction class: %s\n', ct); + end + end +end +% }}} +function md = instantiate_hydrology(md, ct, verbose) % {{{ + switch ct + case 'hydrologyshreve', md.hydrology = hydrologyshreve(); + case 'hydrologydc', md.hydrology = hydrologydc(); + case 'hydrologyglads', md.hydrology = hydrologyglads(); + case 'hydrologypism', md.hydrology = hydrologypism(); + case 'hydrologyshakti', md.hydrology = hydrologyshakti(); + case 'hydrologytws' + if exist('hydrologytws','class') + md.hydrology = hydrologytws(); + end + case 'hydrologyarmapw' + if exist('hydrologyarmapw','class') + md.hydrology = hydrologyarmapw(); + end + otherwise + if verbose + fprintf(' [WARN] unknown hydrology class: %s\n', ct); + end + end +end % }}} +function ct = get_group_classtype(gid) % {{{ +% ===================================================================== +% Utility helpers +% ===================================================================== + ct = ''; + try + ct = netcdf.getAtt(gid, netcdf.getConstant('NC_GLOBAL'), 'classtype'); + catch + end +end +% }}} +function val = get_group_int_att(gid, att_name, default_val) % {{{ + try + val = double(netcdf.getAtt(gid, netcdf.getConstant('NC_GLOBAL'), att_name)); + catch + val = default_val; + end end +% }}} diff --git a/src/m/netcdf/read_netCDF.py b/src/m/netcdf/read_netCDF.py index f2d5001ca..a59f56c45 100644 --- a/src/m/netcdf/read_netCDF.py +++ b/src/m/netcdf/read_netCDF.py @@ -1,679 +1,501 @@ -# imports -from netCDF4 import Dataset -import numpy as np -import numpy.ma as ma -from os import path, remove -from model import * -import re -from results import * -from m1qn3inversion import m1qn3inversion -from taoinversion import taoinversion -from collections import OrderedDict -import sys -from massfluxatgate import massfluxatgate +""" +read_netCDF – Load an ISSM model from a NetCDF4 file written by + write_netCDF.py (Python) or write_netCDF.m (MATLAB). +Usage +----- + from read_netCDF import read_netCDF + md = read_netCDF('model.nc') + md = read_netCDF('model.nc', verbose=True) +Every call returns an independent model() instance (fully re-entrant, +no module-level global state). +""" -''' -Given a NetCDF4 file, this set of functions will perform the following: - 1. Enter each group of the file. - 2. For each variable in each group, update an empty model with the variable's data - 3. Enter nested groups and repeat -''' +import sys +from collections import OrderedDict +import numpy as np +from netCDF4 import Dataset +from model import model +from results import results, solution, solutionstep +from toolkits import toolkits -# make a model framework to fill that is in the scope of this file -model_copy = model() -# make a blank requested output for hydrology -hydro_rout_list = ['default'] +# --------------------------------------------------------------------------- +# Public entry point +# --------------------------------------------------------------------------- -def read_netCDF(filename, verbose = False): +def read_netCDF(filename: str, verbose: bool = False): + """Read *filename* and return a populated ISSM model object.""" if verbose: - print('NetCDF42C v1.2.0') + print('read_netCDF v2.0 (Python)') - ''' - filename = path and name to save file under - verbose = T/F = show or muted log statements. Naturally muted - ''' + from os.path import exists + if not exists(filename): + raise FileNotFoundError(f'read_netCDF: file not found: {filename}') - # this is a precaution so that data is not lost + md = model() + + nc = Dataset(filename, 'r') + nc.set_auto_mask(False) try: - # check if path exists - if path.exists(filename): - if verbose: - print('Opening {} for reading'.format(filename)) - else: pass - - # open the given netCDF4 file - NCData = Dataset(filename, 'r') - # remove masks from numpy arrays for easy conversion - NCData.set_auto_mask(False) - else: - return 'The file you entered does not exist or cannot be found in the current directory' - - # in order to handle some subclasses in the results class, we have to utilize this band-aid - # there will likely be more band-aids added unless a class name library is created with all class names that might be added to a md - try: - # if results has meaningful data, save the name of the subclass and class instance - NCData.groups['results'] - make_results_subclasses(NCData, verbose) - except: - pass - - # similarly, we need to check and see if we have an m1qn3inversion class instance - try: - NCData.groups['inversion'] - check_inversion_class(NCData, verbose) - except: - pass - - # check the smb class used - try: - NCData.groups['smb'] - check_smb_class(NCData, verbose) - except: - pass - - # check the friction class used - try: - NCData.groups['friction'] - check_friction_class(NCData, verbose) - except: - pass - - # check the hydrology class used + md = _read_all_groups(nc, md, verbose) + finally: + nc.close() + + if verbose: + print('Model successfully loaded from NetCDF4.') + return md + + +# --------------------------------------------------------------------------- +# Top-level group walk +# --------------------------------------------------------------------------- + +def _read_all_groups(nc, md, verbose: bool): + for gname, grp in nc.groups.items(): + if verbose: + print(f' Reading group: {gname}') + + ct = _get_classtype(grp) + + # results group: walk solution sub-groups and build results()/solution() objects + if gname == 'results': + md.results = _read_results_group(grp, verbose) + continue + + # toolkits group: each analysis sub-group holds an OrderedDict of solver options + if gname == 'toolkits': + md.toolkits = _read_toolkits_group(grp, verbose) + continue + + # Instantiate the correct sub-class for polymorphic model fields + md = _instantiate_subclass(md, gname, ct, verbose) + + # Populate all fields try: - NCData.groups['hydrology'] - hydro_rout_list = check_hydrology_class(NCData, verbose) - except: - pass - - - - # walk through each group looking for subgroups and variables - for group in NCData.groups.keys(): - if 'debris' in group: - pass + target = getattr(md, gname) + except AttributeError: + if verbose: + print(f' [SKIP] md.{gname} does not exist in model()') + continue + + target = _read_group_into_obj(grp, target, verbose) + setattr(md, gname, target) + + return md + + +# --------------------------------------------------------------------------- +# Results group reader – builds results() / solution() / solutionstep() +# --------------------------------------------------------------------------- + +def _read_results_group(grp, verbose: bool): + """Reconstruct md.results as a proper results()/solution()/solutionstep() + hierarchy, matching what loadresultsfromdisk produces.""" + res = results() + + for sol_name, sol_grp in grp.groups.items(): + ct = _get_classtype(sol_grp) + + if ct == 'struct': + # Each step_k sub-group → one solutionstep + steps = [] + for sg_name in sorted(sol_grp.groups.keys(), + key=lambda n: int(n.split('_')[-1]) + if n.startswith('step_') else 0): + sg = sol_grp.groups[sg_name] + step = solutionstep() + for vname, var in sg.variables.items(): + setattr(step, vname, _read_variable(sg, vname, var, verbose)) + # Recurse into any nested groups within a step (rare) + for sub_name, sub_sg in sg.groups.items(): + setattr(step, sub_name, _read_group_into_obj(sub_sg, {}, verbose)) + steps.append(step) + + if not steps: + # Empty struct – store as empty solution + setattr(res, sol_name, solution([])) + elif len(steps) == 1 and sol_name != 'TransientSolution': + # Single-step non-transient: store the solutionstep directly, + # wrapped in a solution so attribute access still works + sol_obj = solution(steps) + setattr(res, sol_name, sol_obj) else: - # have to send a custom name to this function: filename.groups['group'] - name = "NCData.groups['" + str(group) + "']" - walk_nested_groups(name, NCData, verbose) - model_copy.hydrology.requested_outputs = hydro_rout_list - - if verbose: - print("Model Successfully Loaded.") - - NCData.close() - - return model_copy + setattr(res, sol_name, solution(steps)) - # just in case something unexpected happens - except Exception as e: - if 'NCData' in locals(): - NCData.close() - raise e - -def make_results_subclasses(NCData, verbose = False): - ''' - There are 3 possible subclasses: solution, solutionstep, resultsdakota. - In the NetCDF file these are saved as a list of strings. Ie, say there are 2 - instances of solution under results, StressbalanceSolution and TransientSolution. - In the NetCDF file we would see solution = "StressbalanceSolution", "TransientSolution" - To deconstruct this, we need to iteratively assign md.results.StressbalanceSolution = solution() - and md.results.TransientSolution = solution() and whatever else. - ''' - # start with the subclasses - for subclass in NCData.groups['results'].variables.keys(): - class_instance = subclass + '()' - - # now handle the instances - for instance in NCData.groups['results'].variables[subclass][:]: - # this is an ndarray of numpy bytes_ that we have to convert to strings - class_instance_name = instance.tobytes().decode('utf-8').strip() - # from here we can make new subclasses named as they were in the model saved - setattr(model_copy.results, class_instance_name, eval(class_instance)) if verbose: - print(f'Successfully created results subclass instance {class_instance} named {class_instance_name}.') + print(f' [results] {sol_name}: {len(steps)}-step solution') + else: + # Unexpected sub-group inside results – fall back to dict + if verbose: + print(f' [results] {sol_name}: unrecognised classtype "{ct}", stored as dict') + setattr(res, sol_name, _read_group_into_obj(sol_grp, {}, verbose)) + return res -def check_inversion_class(NCData, verbose = False): - # get the name of the inversion class: either inversion or m1qn3inversion or taoinversion - inversion_class_is = NCData.groups['inversion'].variables['inversion_class_name'][:][...].tobytes().decode() - if inversion_class_is == 'm1qn3inversion': - # if it is m1qn3inversion we need to instantiate that class since it's not native to model() - model_copy.inversion = m1qn3inversion(model_copy.inversion) - if verbose: - print('Conversion successful') - elif inversion_class_is == 'taoinversion': - # if it is taoinversion we need to instantiate that class since it's not native to model() - model_copy.inversion = taoinverion() - if verbose: - print('Conversion successful') - else: pass - -def check_smb_class(NCData, verbose = False): - # get the name of the smb class, either: SMBforcing, SMB, SMBcomponents, SMBd18opdd, - # SMBdebrisEvatt, SMBgemb, SMBgradients, SMBgradientscomponents, SMBgradientsela, - # SMBhenning, SMBmeltcomponents, SMBpdd, SMBpddSicopolis, or SMBsemic - # - # Note: SMB, SMBdebrisEvatt, SMBgemb, SMBgradientscomponents, SMBgradientsela, SMBhenning, - # SMBpddSicopolis, and SMBsemic do not have python versions as of yet - smb_class_is = NCData.groups['smb'].variables['smb_class_name'][:][...].tobytes().decode() - if smb_class_is == 'SMBforcing': - # if it is SMBforcing there is no real need to do anything since that is the default, - # but I have done it anyway to match everything else - model_copy.smb = SMBforcing() - if verbose: - print('Conversion successful') - elif smb_class_is == 'SMBcomponents': - # if it is SMBcomponents load the SMBcomponents module - model_copy.smb = SMBcomponents() - if verbose: - print('Conversion successful') - elif smb_class_is == 'SMBd18opdd': - # if it is SMBd18opdd load the SMBd18opdd module - model_copy.smb = SMBd18opdd() - if verbose: - print('Conversion successful') - elif smb_class_is == 'SMBgradients': - # if it is SMBgradients load the SMBgradients module - model_copy.smb = SMBgradients() - if verbose: - print('Conversion successful') - elif smb_class_is == 'SMBmeltcomponents': - # you guessed it! if it is the SMBmeltcomponents module, load SMBmeltcomponents - model_copy.smb = SMBmeltcomponents() - if verbose: - print('Conversion successful') - elif smb_class_is == 'SMBpdd': - # one more for luck: if it is SMBpdd, load SMBpdd - model_copy.smb = SMBpdd() - if verbose: - print('Conversion successful') - else: - print('Conversion unsuccessful, SMB requested is not yet in Python!') - pass - -def check_friction_class(NCData, verbose = False): - # get the name of the friction class, either: friction(default), frictioncoulomb, frictioncoulomb2, - # frictionhydro, frictionjosh, frictionpism, frictionregcoulomb, frictionregcoulomb2,frictionschoof, - # frictionshakti, frictiontemp, frictiontsai, frictionwaterlayer, frictionweertman, or frictionweertmantemp - # - # Note: frictionregcoulomb, frictionregcoulomb, frictiontemp, frictiontsai, and frictionweertmantemp either - # are not present or tripped an error - friction_class_is = NCData.groups['friction'].variables['friction_class_name'][:][...].tobytes().decode() - if friction_class_is == 'friction': - # if it is friction, load friction (the default) - model_copy.friction = friction() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictioncoulomb': - # if it is frictioncoulomb, import the module, load frictioncoulomb - from frictioncoulomb import frictioncoulomb - model_copy.friction = frictioncoulomb() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictioncoulomb2': - # if it is frictioncoulomb2, import the module, load frictioncoulomb2 - from frictioncoulomb2 import frictioncoulomb2 - model_copy.friction = frictioncoulomb2() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionhydro': - # if it is frictionhydro, import the module, load frictionhydro - from frictionhydro import frictionhydro - model_copy.friction = frictionhydro() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionjosh': - # if it is frictionjosh, import the module, load frictionjosh - from frictionjosh import frictionjosh - model_copy.friction = frictionjosh() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionpism': - # if it is frictionpism, import the module, load frictionpism - from frictionpism import frictionpism - model_copy.friction = frictionpism() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionschoof': - # if it is frictionschoof, import the module, load frictionschoof - from frictionschoof import frictionschoof - model_copy.friction = frictionschoof() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionshakti': - #if it is frictionshakti, import the module, load frictionshakti - from frictionshakti import frictionshakti - model_copy.friction = frictionshakti() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionwaterlayer': - #if it is frictionwaterlayer, import the module, load frictionwaterlayer - from frictionwaterlayer import frictionwaterlayer - model_copy.friction = frictionwaterlayer() - if verbose: - print('Conversion successful') - elif friction_class_is == 'frictionweertman': - #if it is frictionweertman, import the module, load frictionweertman - from frictionweertman import frictionweertman - model_copy.friction = frictionweertman() - if verbose: - print('Conversion successful') - else: - print('Conversion unsuccessful, friction requested is not present and/or working yet in Python!') - pass - -def check_hydrology_class(NCData, verbose = False): - # get the name of the hydrology class, either: hydrologyshreve (default), hydrology, hydrologyarmapw, - # hydrologydc, hydrologyglads, hydrologypism, hydrologyshakti, hydrologtws - # - # Note: hydrology, hydrologyarmapw, hydrologytws do not have a python implementation - hydrology_class_is = NCData.groups['hydrology'].variables['hydrology_class_name'][:][...].tobytes().decode() - hydro_rout = NCData.groups['hydrology'].variables['requested_outputs'][:][...].tobytes().decode() - hydro_rout_list = hydro_rout.split() - model_copy.hydrology.requested_outputs = hydro_rout_list - if hydrology_class_is == 'hydrologyshreve': - #if it is hydrologyshreve, load hydrologyshreve - model_copy.hydrology = hydrologyshreve() - if verbose: - print('Conversion successful') - elif hydrology_class_is == 'hydrologydc': - #if it is hydrologydc, load hydrologydc - model_copy.hydrology = hydrologydc() - if verbose: - print('Conversion successful') - elif hydrology_class_is == 'hydrologyglads': - #if it is hydrologyglads, load hydrologyglads - model_copy.hydrology = hydrologyglads() - if verbose: - print('Conversion successful') - elif hydrology_class_is == 'hydrologypism': - #if it is hydrologypism, load hydrologypism - model_copy.hydrology = hydrologypism() - if verbose: - print('Conversion successful') - elif hydrology_class_is == 'hydrologyshakti': - #if it is hydrologyshakti, load hydrologyshakti - model_copy.hydrology = hydrologyshakti() - if verbose: - print('Conversion successful') - else: - print('Conversion unsuccessful, hydrology requested is not present and/or working yet in Python!') - pass - return hydro_rout_list - -def walk_nested_groups(group_location_in_file, NCData, verbose = False): - # first, we enter the group by: filename.groups['group_name'] - # second we search the current level for variables: filename.groups['group_name'].variables.keys() - # at this step we check for multidimensional structure arrays/ arrays of objects and filter them out - # third we get nested group keys by: filename.groups['group_name'].groups.keys() - # if a nested groups exist, repeat all - - for variable in eval(group_location_in_file + '.variables.keys()'): - if 'is_object' not in locals(): - if variable == 'this_is_a_nested' and 'results' in group_location_in_file and 'qmu' not in group_location_in_file: - # have to do some string deconstruction to get the name of the class instance/last group from 'NetCDF.groups['group1'].groups['group1.1']' - pattern = r"\['(.*?)'\]" - matches = re.findall(pattern, group_location_in_file) - name_of_struct = matches[-1] #eval(group_location_in_file + ".variables['solution']") - deserialize_nested_results_struct(group_location_in_file, name_of_struct, NCData) - is_object = True - - elif variable == 'name_of_cell_array': - # reconstruct an array of elements - deserialize_array_of_objects(group_location_in_file, model_copy, NCData, verbose) - is_object = True - - elif variable == 'this_is_a_nested' and 'qmu' in group_location_in_file: - if verbose: - print('encountered qmu structure that is not yet supported.') - else: pass - - is_object = True - - else: - location_of_variable_in_file = group_location_in_file + ".variables['" + str(variable) + "']" - # group_location_in_file is like filename.groups['group1'].groups['group1.1'].groups['group1.1.1'] - # Define the regex pattern to match the groups within brackets - pattern = r"\['(.*?)'\]" - # Use regex to find all matches and return something like 'group1.group1.1.group1.1.1 ...' where the last value is the name of the variable - matches = re.findall(pattern, location_of_variable_in_file) - variable_name = matches[-1] - location_of_variable_in_model = '.'.join(matches[:-1]) - deserialize_data(location_of_variable_in_file, location_of_variable_in_model, variable_name, NCData, verbose=verbose) - - # if one of the variables above was an object, further subclasses will be taken care of when reconstructing it - if 'is_object' in locals(): - pass - else: - for nested_group in eval(group_location_in_file + '.groups.keys()'): - new_nested_group = group_location_in_file + ".groups['" + str(nested_group) + "']" - walk_nested_groups(new_nested_group, NCData, verbose=verbose) - - - -''' - MATLAB has Multidimensional Structure Arrays in 2 known classes: results and qmu. - The python classes results.py and qmu.py emulate this MATLAB object in their own - unique ways. The functions in this script will assign data to either of these - classes such that the final structure is compatible with its parent class. -''' - -def deserialize_nested_results_struct(group_location_in_file, name_of_struct, NCData, verbose = False): - ''' - A common multidimensional array is the 1xn md.results.TransientSolution object. - - The way that this object emulates the MATLAB mutli-dim. struct. array is with - the solution().steps attr. which is a list of solutionstep() instances - The process to recreate is as follows: - 1. Get instance of solution() with solution variable (the instance is made in make_results_subclasses) - 2. For each subgroup, create a solutionstep() class instance - 2a. Populate the instance with the key:value pairs - 2b. Append the instance to the solution().steps list - ''' - # step 1 - class_instance_name = name_of_struct - - # for some reason steps is not already a list - setattr(model_copy.results.__dict__[class_instance_name], 'steps', list()) - - steps = model_copy.results.__dict__[class_instance_name].steps - - # step 2 - layer = 1 - for subgroup in eval(group_location_in_file + ".groups.keys()"): - solutionstep_instance = solutionstep() - # step 2a - subgroup_location_in_file = group_location_in_file + ".groups['" + subgroup + "']" - for key in eval(subgroup_location_in_file + ".variables.keys()"): - value = eval(subgroup_location_in_file + ".variables['" + str(key) + "'][:]") - setattr(solutionstep_instance, key, value) - # step 2b - steps.append(solutionstep_instance) + +# --------------------------------------------------------------------------- +# Toolkits group reader – rebuilds toolkits() with OrderedDict per analysis +# --------------------------------------------------------------------------- + +def _read_toolkits_group(grp, verbose: bool): + """Reconstruct md.toolkits as a toolkits() instance whose analysis + attributes are OrderedDicts of solver options (matching mumpsoptions() etc.)""" + tk = toolkits.__new__(toolkits) # bypass __init__ / setdefaultparameters + tk.DefaultAnalysis = None + tk.RecoveryAnalysis = None + + for aname, agrp in grp.groups.items(): + ct = _get_classtype(agrp) + opts = OrderedDict() + + if ct == 'struct': + # Variables are stored one level deeper in step_1 + step_grps = agrp.groups + src = step_grps.get('step_1', agrp) if step_grps else agrp + else: + src = agrp + + for vname, var in src.variables.items(): + opts[vname] = _read_variable(src, vname, var, verbose) + + setattr(tk, aname, opts) if verbose: - print('Succesfully loaded layer ' + str(layer) + ' to results.' + str(class_instance_name) + ' struct.') - else: pass - layer += 1 + print(f' [toolkits] {aname}: {list(opts.keys())}') - if verbose: - print('Successfully recreated results structure ' + str(class_instance_name)) - - - -def deserialize_array_of_objects(group_location_in_file, model_copy, NCData, verbose): - ''' - The structure in netcdf for groups with the name_of_cell_array variable is like: - - group: 2x6_cell_array_of_objects { - name_of_cell_array = - - group: Row_1_of_2 { - group: Col_1_of_6 { - ... other groups can be here that refer to objects - } // group Col_6_of_6 - } // group Row_1_of_2 - - group: Row_2_of_2 { - group: Col_1_of_6 { - ... other groups can be here that refer to objects - } // group Col_6_of_6 - } // group Row_2_of_2 - } // group 2x6_cell_array_of_objects - - We have to navigate this structure to extract all the data and recreate the - original structure when the model was saved - ''' - - if verbose: - print(f"Loading array of objects.") - - # get the name_of_cell_array, rows and cols vars - name_of_cell_array_varID = eval(group_location_in_file + ".variables['name_of_cell_array']") - rows_varID = eval(group_location_in_file + ".variables['rows']") - cols_varID = eval(group_location_in_file + ".variables['cols']") - - name_of_cell_array = name_of_cell_array_varID[:][...].tobytes().decode() - rows = rows_varID[:] - cols = cols_varID[:] - - # now we work backwards: make the array, fill it in, and assign it to the model - - # make the array - array = list() - - subgroups = eval(group_location_in_file + ".groups") #.keys()") - - # enter each subgroup, get the data, assign it to the corresponding index of cell array - if rows > 1: - # we go over rows - # set index for rows - row_idx = 0 - for row in list(subgroups): - # make list for each row - current_row = list() - columns = subgroups[str(row)].groups.keys() - - # set index for columns - col_idx = 0 - - # iterate over columns - for col in list(columns): - # now get the variables - current_col_vars = columns.groups[str(col)].variables - - # check for special datastructures - if "class_is_a" in current_col_vars: - class_name = subgroups[str(col)].variables['class_is_a'][:][...].tobytes().decode() - col_data = deserialize_class_instance(class_name, columns.groups[str(col)], NCData, verbose) - is_object = True - elif "this_is_a_nested" in current_col_vars: - # functionality not yet supported - print('Error: Cell Arrays of structs not yet supported!') - is_object = True - else: - if 'is_object_' in locals(): - pass - # already taken care of - else: - # store the variables as normal -- to be added later - print('Error: Arrays of mixed objects not yet supported!') - for var in current_col_vars: - # this is where that functionality would be handled - pass - col_idx += 1 - # add the entry to our row list - current_row.append(col_data) - - # add the list of columns to the array - array.append(current_row) - row_idx += 1 + return tk - else: - # set index for columns - col_idx = 0 - - # iterate over columns - for col in list(subgroups): - # now get the variables - current_col_vars = subgroups[str(col)].variables - - # check for special datastructures - if "class_is_a" in current_col_vars: - class_name = subgroups[str(col)].variables['class_is_a'][:][...].tobytes().decode() - col_data = deserialize_class_instance(class_name, subgroups[str(col)], NCData, verbose) - is_object = True - elif "this_is_a_nested" in current_col_vars: - # functionality not yet supported - print('Error: Cell Arrays of structs not yet supported!') - is_object = True - else: - if 'is_object_' in locals(): - pass - # already taken care of - else: - # store the variables as normal -- to be added later - print('Error: Arrays of mixed objects not yet supported!') - for var in current_col_vars: - # this is where that functionality would be handled - pass - col_idx += 1 - # add the list of columns to the array - array.append(col_data) - - # finally, add the attribute to the model - pattern = r"\['(.*?)'\]" - matches = re.findall(pattern, group_location_in_file) - variable_name = matches[0] - setattr(model_copy.__dict__[variable_name], name_of_cell_array, array) - if verbose: - print(f"Successfully loaded array of objects: {name_of_cell_array} to {variable_name}") +# --------------------------------------------------------------------------- +# Generic recursive group reader +# --------------------------------------------------------------------------- +def _read_group_into_obj(grp, obj, verbose: bool): + """Populate *obj* from the variables and sub-groups in *grp*.""" + # Read variables at this level + for vname, var in grp.variables.items(): + data = _read_variable(grp, vname, var, verbose) + obj = _set_attr(obj, vname, data, verbose) -def deserialize_class_instance(class_name, group, NCData, verbose=False): + # Recurse into sub-groups + for sg_name, sg in grp.groups.items(): + ct = _get_classtype(sg) - if verbose: - print(f"Loading class: {class_name}") + if ct == 'struct': + # Struct array (results.TransientSolution etc.) + nsteps = int(getattr(sg, 'nsteps', len(sg.groups))) + data = _read_struct_array(sg, nsteps, verbose) - # this function requires the class module to be imported into the namespace of this file. - # we make a custom error in case the class module is not in the list of imported classes. - # most ISSM classes are imported by from import - class ModuleError(Exception): - pass - - if class_name not in sys.modules: - raise ModuleError(str('Model requires the following class to be imported from a module: ' + class_name + ". Please add the import to read_netCDF.py in order to continue.")) - - # Instantiate the class - class_instance = eval(class_name + "()") - - # Get and assign properties - subgroups = list(group.groups.keys()) - - if len(subgroups) == 1: - # Get properties - subgroup = group[subgroups[0]] - varIDs = subgroup.variables.keys() - for varname in varIDs: - # Variable metadata - var = subgroup[varname] - - # Data - if 'char' in var.dimensions[0]: - data = var[:][...].tobytes().decode() - else: - data = var[:] + elif ct == 'cell_of_objects': + data = _read_cell_of_objects(sg, verbose) + + elif ct == 'dict': + data = _read_dict(sg, verbose) - # Some classes may have permissions, so we skip those + else: + # Regular sub-class or unknown: read into existing attr or new dict try: - setattr(class_instance, varname, data) - except: - pass - else: - # Not supported - pass + sub_obj = getattr(obj, sg_name) + except AttributeError: + sub_obj = {} + + # If the classtype tells us which class to instantiate + if ct and ct not in ('', 'struct', 'cell_of_objects', 'dict'): + sub_obj = _instantiate_class(ct, sub_obj, verbose) + + data = _read_group_into_obj(sg, sub_obj, verbose) + + obj = _set_attr(obj, sg_name, data, verbose) + + return obj + + +# --------------------------------------------------------------------------- +# Struct-array reader (results.TransientSolution etc.) +# --------------------------------------------------------------------------- + +def _read_struct_array(grp, nsteps: int, verbose: bool): + """Return a list of dicts (one per step_k sub-group).""" + steps = [] + for sg_name in sorted(grp.groups.keys(), + key=lambda n: int(n.split('_')[-1]) if n.startswith('step_') else 0): + sg = grp.groups[sg_name] + step = {} + for vname, var in sg.variables.items(): + step[vname] = _read_variable(sg, vname, var, verbose) + # Recurse into any nested groups within a step (rare) + for sub_name, sub_grp in sg.groups.items(): + step[sub_name] = _read_group_into_obj(sub_grp, {}, verbose) + steps.append(step) + + if verbose: + print(f' [struct] loaded {len(steps)}-step struct array') + return steps + + +# --------------------------------------------------------------------------- +# Cell-of-objects reader +# --------------------------------------------------------------------------- + +def _read_cell_of_objects(grp, verbose: bool): + nrows = int(getattr(grp, 'nrows', 1)) + ncols = int(getattr(grp, 'ncols', len(grp.groups))) + result = [[None] * ncols for _ in range(nrows)] + + for ig_name, ig in grp.groups.items(): + import re + m = re.match(r'^item_(\d+)_(\d+)$', ig_name) + if not m: + continue + r = int(m.group(1)) - 1 + c = int(m.group(2)) - 1 + ct = _get_classtype(ig) + obj = _instantiate_class(ct, {}, verbose) + obj = _read_group_into_obj(ig, obj, verbose) + result[r][c] = obj + + if nrows == 1: + result = result[0] # return flat list for 1-D case + if verbose: + print(f' [cell] loaded {nrows}x{ncols} cell of objects') + return result + + +# --------------------------------------------------------------------------- +# Dict reader +# --------------------------------------------------------------------------- + +def _read_dict(grp, verbose: bool) -> dict: + d = {} + for vname, var in grp.variables.items(): + d[vname] = _read_variable(grp, vname, var, verbose) + for sg_name, sg in grp.groups.items(): + d[sg_name] = _read_group_into_obj(sg, {}, verbose) + return d + - if verbose: - print(f"Successfully loaded class instance {class_name} to model") - return class_instance +# --------------------------------------------------------------------------- +# Variable reader +# --------------------------------------------------------------------------- +def _read_variable(grp, vname: str, var, verbose: bool): + """Read one NetCDF variable and convert to the appropriate Python type.""" + type_is = getattr(var, 'type_is', '') + # Empty sentinel + if type_is == 'empty': + return [] -def deserialize_data(location_of_variable_in_file, location_of_variable_in_model, variable_name, NCData, verbose = False): - # as simple as navigating to the location_of_variable_in_model and setting it equal to the location_of_variable_in_file - # NetCDF4 has a property called "_FillValue" that sometimes saves empty lists, so we have to catch those - FillValue = -9223372036854775806 + raw = var[:] + + # Bool + if type_is == 'bool': + return np.array(raw, dtype=bool) + + # String / cell-of-strings + if type_is == 'string' or (hasattr(var, 'dimensions') and + any('char' in d for d in var.dimensions)): + return _decode_string(raw) + + if type_is == 'cell_of_strings': + return _decode_cell_strings(raw) + + # Numpy array / scalar + data = np.array(raw) + + # Squeeze trailing size-1 dims + data = np.squeeze(data) + + # Single-element array → Python scalar + if data.ndim == 0: + v = data.item() + # Return int if value is a whole number (guard against nan/inf first) + if isinstance(v, float) and np.isfinite(v) and v == int(v): + return int(v) + return v + + return data + + +def _decode_string(raw) -> str: + """Convert raw NetCDF char data to a Python str.""" try: - # results band-aid... - if str(location_of_variable_in_model + '.' + variable_name) in ['results.solutionstep', 'results.solution', 'results.resultsdakota']: - pass - # qmu band-aid - elif 'qmu.statistics.method' in str(location_of_variable_in_model + '.' + variable_name): - pass - # handle any strings: - elif 'char' in eval(location_of_variable_in_file + '.dimensions[0]'): - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:][...].tobytes().decode()')) - # handle ndarrays + lists - elif len(eval(location_of_variable_in_file + '[:]'))>1: - # check for bool - try: # there is only one datatype assigned the attribute 'units' and that is bool, so anything else will go right to except - if eval(location_of_variable_in_file + '.units') == 'bool': - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, np.array(eval(location_of_variable_in_file + '[:]'), dtype = bool)) - else: - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) - except: - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) - # catch everything else - else: - # check for FillValue. use try/except because try block will only work on datatypes like int64, float, single element lists/arrays ect and not nd-arrays/n-lists etc - try: - # this try block will only work on single ints/floats/doubles and will skip to the except block for all other cases - var_to_save = eval(location_of_variable_in_file + '[:][0]') # note the [0] on the end - if FillValue == var_to_save: - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, []) - else: - if var_to_save.is_integer(): - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, int(var_to_save)) - else: - # we have to convert numpy datatypes to native python types with .item() - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, var_to_save.item()) - except: - setattr(eval('model_copy.' + location_of_variable_in_model), variable_name, eval(location_of_variable_in_file + '[:]')) - except AttributeError: - deserialize_dict(location_of_variable_in_file, location_of_variable_in_model, NCData, verbose=verbose) + return raw.tobytes().decode('utf-8').strip('\x00').strip() + except Exception: + try: + return ''.join(c.decode('utf-8') if isinstance(c, bytes) else c + for c in np.ndarray.flatten(raw)) + except Exception: + return str(raw) + + +def _decode_cell_strings(raw) -> list: + """Convert a 2-D char array to a list of strings.""" + if raw.ndim == 1: + return [_decode_string(raw)] + result = [] + for row in raw: + s = _decode_string(row) + result.append(s) + return result + + +# --------------------------------------------------------------------------- +# Subclass instantiation +# --------------------------------------------------------------------------- + +def _instantiate_subclass(md, gname: str, ct: str, verbose: bool): + """For polymorphic model fields, swap in the right class instance.""" + if not ct: + return md + + try: + if gname == 'inversion': + if ct == 'm1qn3inversion': + from m1qn3inversion import m1qn3inversion + md.inversion = m1qn3inversion() + elif ct == 'taoinversion': + from taoinversion import taoinversion + md.inversion = taoinversion() + elif gname == 'smb': + md = _instantiate_smb(md, ct, verbose) + + elif gname == 'friction': + md = _instantiate_friction(md, ct, verbose) + + elif gname == 'hydrology': + md = _instantiate_hydrology(md, ct, verbose) + + elif gname == 'mesh': + if ct == 'mesh3dprisms': + from mesh3dprisms import mesh3dprisms + md.mesh = mesh3dprisms() + + # All other groups: keep default instance, just overwrite fields below + + except Exception as e: + if verbose: + print(f' [WARN] could not instantiate {ct} for {gname}: {e}') + + return md + + +def _instantiate_class(ct: str, fallback, verbose: bool): + """Try to instantiate *ct* by importing it; return *fallback* on failure.""" + if not ct or ct in ('struct', 'cell_of_objects', 'dict', ''): + return fallback + try: + mod = __import__(ct, fromlist=[ct]) + cls = getattr(mod, ct) + return cls() + except Exception: + pass + # Try to find ct in already-imported modules + for mod_name, mod in list(sys.modules.items()): + if mod is not None and hasattr(mod, ct): + cls = getattr(mod, ct) + if isinstance(cls, type): + try: + return cls() + except Exception: + pass if verbose: - print('Successfully loaded ' + location_of_variable_in_model + '.' + variable_name + ' into model.') - - - -def deserialize_dict(location_of_variable_in_file, location_of_variable_in_model, NCData, verbose = False): - FillValue = -9223372036854775806 - - # the key will be the last item in the location - key = ''.join(location_of_variable_in_model.split('.')[-1]) - - # update the location to point to the dict instead of the dict key - location_of_variable_in_model = '.'.join(location_of_variable_in_model.split('.')[:-1]) - - # verify we're working with a dict: - if isinstance(eval('model_copy.' + location_of_variable_in_model), OrderedDict): - dict_object = eval('model_copy.' + location_of_variable_in_model) - - # handle any strings: - if 'char' in eval(location_of_variable_in_file + '.dimensions[0]'): - data = eval(location_of_variable_in_file + '[:][...].tobytes().decode()') - dict_object.update({key: data}) - - # handle ndarrays + lists - elif len(eval(location_of_variable_in_file + '[:]'))>1: - # check for bool - try: # there is only one datatype assigned the attribute 'units' and that is bool, so anything else will go right to except - if eval(location_of_variable_in_file + '.units') == 'bool': - data = np.array(eval(location_of_variable_in_file + '[:]'), dtype = bool) - dict_object.update({key: data}) - else: - data = eval(location_of_variable_in_file + '[:]') - dict_object.update({key: data}) - except: - data = eval(location_of_variable_in_file + '[:]') - dict_object.update({key: data}) - # catch everything else - else: - # check for FillValue. use try/except because try block will only work on datatypes like int64, float, single element lists/arrays ect and not nd-arrays/n-lists etc - try: - # this try block will only work on single ints/floats/doubles and will skip to the except block for all other cases - if FillValue == eval(location_of_variable_in_file + '[:][0]'): - dict_object.update({key: []}) - else: - # we have to convert numpy datatypes to native python types with .item() - var_to_save = eval(location_of_variable_in_file + '[:][0]') # note the [0] on the end - dict_object.update({key: var_to_save.item()}) - except: - data = eval(location_of_variable_in_file + '[:]') - dict_object.update({key: data}) + print(f' [WARN] could not instantiate class {ct}, using fallback') + return fallback + + +def _instantiate_smb(md, ct: str, verbose: bool): + mapping = { + 'SMBforcing': ('SMBforcing', 'SMBforcing'), + 'SMBpdd': ('SMBpdd', 'SMBpdd'), + 'SMBd18opdd': ('SMBd18opdd', 'SMBd18opdd'), + 'SMBgradients': ('SMBgradients', 'SMBgradients'), + 'SMBcomponents': ('SMBcomponents', 'SMBcomponents'), + 'SMBmeltcomponents': ('SMBmeltcomponents', 'SMBmeltcomponents'), + 'SMBgradientscomponents':('SMBgradientscomponents','SMBgradientscomponents'), + 'SMBgradientsela': ('SMBgradientsela', 'SMBgradientsela'), + 'SMBhenning': ('SMBhenning', 'SMBhenning'), + 'SMBgemb': ('SMBgemb', 'SMBgemb'), + 'SMBpddSicopolis': ('SMBpddSicopolis', 'SMBpddSicopolis'), + 'SMBsemic': ('SMBsemic', 'SMBsemic'), + 'SMBdebrisEvatt': ('SMBdebrisEvatt', 'SMBdebrisEvatt'), + } + if ct in mapping: + mod_name, cls_name = mapping[ct] + try: + mod = __import__(mod_name, fromlist=[cls_name]) + md.smb = getattr(mod, cls_name)() + except Exception as e: + if verbose: + print(f' [WARN] could not load SMB class {ct}: {e}') + else: + if verbose: + print(f' [WARN] unknown SMB class: {ct}') + return md + + +def _instantiate_friction(md, ct: str, verbose: bool): + known = [ + 'friction', 'frictioncoulomb', 'frictioncoulomb2', + 'frictionhydro', 'frictionjosh', 'frictionpism', + 'frictionregcoulomb', 'frictionregcoulomb2', 'frictionschoof', + 'frictionshakti', 'frictiontsai', 'frictionwaterlayer', + 'frictionweertman', 'frictionweertmantemp', + ] + if ct in known: + try: + mod = __import__(ct, fromlist=[ct]) + md.friction = getattr(mod, ct)() + except Exception as e: + if verbose: + print(f' [WARN] could not load friction class {ct}: {e}') + else: + if verbose: + print(f' [WARN] unknown friction class: {ct}') + return md + + +def _instantiate_hydrology(md, ct: str, verbose: bool): + known = [ + 'hydrologyshreve', 'hydrologydc', 'hydrologyglads', + 'hydrologypism', 'hydrologyshakti', 'hydrologytws', 'hydrologyarmapw', + ] + if ct in known: + try: + mod = __import__(ct, fromlist=[ct]) + md.hydrology = getattr(mod, ct)() + except Exception as e: + if verbose: + print(f' [WARN] could not load hydrology class {ct}: {e}') else: - print(f"Unrecognized object was saved to NetCDF file and cannot be reconstructed: {location_of_variable_in_model}") \ No newline at end of file + if verbose: + print(f' [WARN] unknown hydrology class: {ct}') + return md + + +# --------------------------------------------------------------------------- +# Attribute setter +# --------------------------------------------------------------------------- + +def _set_attr(obj, name: str, value, verbose: bool): + """Set attribute *name* on *obj* (works for class instances and dicts).""" + try: + if isinstance(obj, dict): + obj[name] = value + else: + setattr(obj, name, value) + except Exception as e: + if verbose: + print(f' [WARN] could not set {name}: {e}') + return obj + + +# --------------------------------------------------------------------------- +# Utility +# --------------------------------------------------------------------------- + +def _get_classtype(grp) -> str: + return getattr(grp, 'classtype', '') diff --git a/src/m/netcdf/write_netCDF.m b/src/m/netcdf/write_netCDF.m index 852d25b1e..97337c2d2 100644 --- a/src/m/netcdf/write_netCDF.m +++ b/src/m/netcdf/write_netCDF.m @@ -1,1012 +1,444 @@ -%{ -Given a model, this set of functions will perform the following: - 1. Enter each nested class of the model. - 2. View each attribute of each nested class. - 3. Compare state of attribute in the model to an empty model class. - 4. If states are identical, pass. - 5. Otherwise, create nested groups named after class structure. - 6. Create variable named after class attribute and assign value to it. -%} - - -function write_netCDF(model_var, filename, varargin) - if nargin > 2 - verbose = true; +function write_netCDF(md, filename, varargin) +%WRITE_NETCDF - Save an ISSM model to a NetCDF4 file. +% +% Usage: +% write_netCDF(md, filename) +% write_netCDF(md, filename, 'verbose', true) +% +% Inputs: +% md - ISSM model object +% filename - path/name for the output .nc file (overwritten if it exists) +% +% Optional name-value pair: +% 'verbose' - true/false (default false) +% +% The file can be read back by read_netCDF.m (MATLAB) or read_netCDF.py (Python). +% +%NetCDF group layout +%------------------- +%Each first-level field of md (mesh, mask, geometry, …) becomes a NetCDF +%group. Fields within those classes are stored as variables inside the +%group. Nested subclasses get nested subgroups. +% +%Class-type encoding +%------------------- +%For fields whose MATLAB class can vary at run-time (inversion, smb, +%friction, hydrology), the concrete class name is stored as a group +%attribute 'classtype' so that the reader can reconstruct the correct +%object. Every group also gets a 'classtype' attribute with the MATLAB +%class name of the object it represents (e.g. "mesh2d", "m1qn3inversion"). +% +%Results storage +%--------------- +%md.results. is a MATLAB struct array (1×n). Each element is +%stored as a subgroup named "step_" inside the solution subgroup. Every +%variable in every step is stored as a numeric or char variable. + % --- parse options ------------------------------------------------ + p = inputParser(); + p.addParameter('verbose', false, @islogical); + % also accept positional 'verbose' flag (legacy) + if nargin > 2 && islogical(varargin{1}) + verbose = varargin{1}; else - verbose = false; + p.parse(varargin{:}); + verbose = p.Results.verbose; end + if verbose - disp('MATLAB C2NetCDF4 v1.1.14'); - end - - % model_var = class object to be saved - % filename = path and name to save file under - - % Create a NetCDF file to write to - NetCDF = make_NetCDF(filename, verbose); - - % Create an instance of an empty model class to compare model_var against - empty_model = model(); - - % Walk through the model_var class and compare subclass states to empty_model - walk_through_model(model_var, empty_model, NetCDF, verbose); - - % in order to handle some subclasses in the results class, we have to utilize this band-aid - % there will likely be more band-aids added unless a class name library is created with all class names that might be added to a model - try - % if results had meaningful data, save the name of the subclass and class instance - netcdf.inqNcid(NetCDF,'results'); - results_subclasses_bandaid(model_var, NetCDF, verbose); - % otherwise, ignore - catch - end - - netcdf.close(NetCDF); - if verbose - disp('Model successfully saved as NetCDF4'); + disp('write_netCDF v2.0 (MATLAB)'); end -end - - -function NetCDF = make_NetCDF(filename, verbose) - % matlab can't handle input in the jupyter interface, so we just yell at the user to rename - % their file if needed - % If file already exists delete / rename it + % --- create / overwrite file -------------------------------------- if exist(filename, 'file') == 2 - fprintf('File %s already exists\n', filename); - disp('Please rename your file.') - return - - % If so, inquire for a new name or to delete the existing file - %newname = input('Give a new name or input "delete" to replace: ', 's'); - - %if strcmpi(newname, 'delete') - %delete filename; - %else - %fprintf('New file name is %s\n', newname); - %filename = newname; - %end - else - % Otherwise create the file and define it globally so other functions can call it - - NetCDF = netcdf.create(filename, 'NETCDF4'); - netcdf.putAtt(NetCDF, netcdf.getConstant('NC_GLOBAL'), 'history', ['Created ', datestr(now)]); - netcdf.defDim(NetCDF, 'Unlim', netcdf.getConstant('NC_UNLIMITED')); % unlimited dimension - netcdf.defDim(NetCDF, 'float', 1); % single integer dimension - netcdf.defDim(NetCDF, 'int', 1); % single float dimension - + delete(filename); if verbose - fprintf('Successfully created %s\n', filename); + fprintf('Existing file %s deleted.\n', filename); end + end - return + ncid = netcdf.create(filename, 'NETCDF4'); + netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'history', ... + ['Created by write_netCDF.m on ' datestr(now)]); + netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'Conventions', 'ISSM-NetCDF-2.0'); + + % Global dimensions reused throughout the file + netcdf.defDim(ncid, 'scalar', 1); % for single-value variables + netcdf.defDim(ncid, 'Unlim', netcdf.getConstant('NC_UNLIMITED')); + + if verbose + fprintf('Created %s\n', filename); end -end + % --- walk all top-level fields of md ------------------------------ + top_fields = fieldnames(md); + for fi = 1:numel(top_fields) + fname = top_fields{fi}; + val = md.(fname); -%{ - Since python uses subclass instances and MATLAB uses fields, we need to guess which subclass instance python will need - given the name of the sub-field in MATLAB. We make this guess based on the name of the MATLAB subfield that will contain - the name of the python subclass instance. For example, md.results.StressbalanceSolution is an subfield in MATLAB, - but a class instance of solution(). Notice that StressbalanceSolution contains the name "Solution" in it. This is what - we will save to the netCDF file for python to pick up. -%} - -function results_subclasses_bandaid(model_var, NetCDF, verbose) - - % The results class may have nested fields within it, so we need to record the name of - % the nested field as it appears in the model that we're trying to save - quality_control = {}; - - % Access the results subclass of model_var - results_var = model_var.results; - - % get the results group id so we can write to it - groupid = netcdf.inqNcid(NetCDF,'results'); - - % Loop through each class instance in results - class_instance_names = fieldnames(results_var); - - % we save lists of instances to the netcdf - solutions = {}; - solutionsteps = {}; - resultsdakotas = {}; - - for i = 1:numel(class_instance_names) - class_instance_name = class_instance_names{i}; - % there are often mutliple instances of the same class/struct so we have to number them - % Check to see if there is a solutionstep class instance - if contains(class_instance_name, 'solutionstep',IgnoreCase=true) - quality_control{end+1} = 1; - solutionsteps{end+1} = class_instance_name; - if verbose - disp('Successfully stored class python subclass instance: solutionstep') - end - end - - % Check to see if there is a solution class instance - if contains(class_instance_name, 'solution',IgnoreCase=true) - quality_control{end+1} = 1; - solutions{end+1} = class_instance_name; - if verbose - disp('Successfully stored class python subclass instance: solution') - end - end - - % Check to see if there is a resultsdakota class instance - if contains(class_instance_name, 'resultsdakota',IgnoreCase=true) - quality_control{end+1} = 1; - resultsdakotas{end+1} = class_instance_name; + if isempty(val) if verbose - disp('Successfully stored class python subclass instance: resultsdakota') + fprintf(' md.%s is empty – skipped\n', fname); end + continue end - end - if ~isempty(solutionsteps) - write_cell_with_strings('solutionstep', solutionsteps, groupid, NetCDF, verbose) - end - if ~isempty(solutions) - write_cell_with_strings('solution', solutions, groupid, NetCDF, verbose) - end - if ~isempty(resultsdakotas) - write_cell_with_strings('resultsdakota', resultsdakotas, groupid, NetCDF, verbose) - end - - - % Check if all class instances were processed correctly - if numel(quality_control) ~= numel(class_instance_names) - disp('Error: The class instance within your model.results class is not currently supported by this application'); - else - if verbose - disp('The results class was successfully stored on disk'); - end - end -end + gid = netcdf.defGrp(ncid, fname); + % Store the concrete class name so readers can reconstruct the right object + netcdf.putAtt(gid, netcdf.getConstant('NC_GLOBAL'), 'classtype', class(val)); + write_group(val, gid, ncid, fname, verbose); + end -function walk_through_model(model_var, empty_model, NetCDF, verbose) - % Iterate over first layer of model_var attributes and assume this first layer is only classes fundamental to the model() class - % note that groups are the same as class instances/subfields in this context - groups = fieldnames(model_var); - for group = 1:numel(groups) - % now this new variable takes the form model.mesh , model.damage etc. - model_subclass = model_var.(groups{group}); - empty_model_subclass = empty_model.(groups{group}); - % Now we can recursively walk through the remaining subclasses - list_of_layers = {groups{group}}; - walk_through_subclasses(model_subclass, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); + netcdf.close(ncid); + if verbose + disp('Model successfully saved as NetCDF4.'); end end - - -function walk_through_subclasses(model_subclass, empty_model_subclass, given_list_of_layers, empty_model, NetCDF, verbose) - % Recursivley iterate over each subclass' attributes and look for more subclasses and variables with relevant data - % model_subclass is an object (ie, md.mesh.elements) - % list_of_layers is a cell array of subclasses/attributes/fields so that we can copy the structure into netcdf (ie, {'mesh', 'elements'}) - % need to check if inversion or m1qn3inversion or taoinversion class - if numel(given_list_of_layers) == 1 - if strcmp(given_list_of_layers{1}, 'inversion') - create_group(model_subclass, given_list_of_layers, NetCDF, verbose); - check_inversion_class(model_subclass, NetCDF, verbose); - - elseif strcmp(given_list_of_layers{1},'smb') - create_group(model_subclass, given_list_of_layers, NetCDF, verbose); - check_smb_class(model_subclass, NetCDF, verbose); - - elseif strcmp(given_list_of_layers{1},'friction') - create_group(model_subclass, given_list_of_layers, NetCDF, verbose); - check_friction_class(model_subclass, NetCDF, verbose); - - elseif strcmp(given_list_of_layers{1},'hydrology') - create_group(model_subclass, given_list_of_layers, NetCDF, verbose); - check_hydrology_class(model_subclass, NetCDF, verbose); - end - end - - % Use try/except since model_subclass is either a subclass/struct w/ props/fields or it's not, no unknown exceptions - try - % look for children - this is where the catch would be called - children = fieldnames(model_subclass); - - % if there are children, loop through them and see if we need to save any data - for child = 1:numel(children) - % record our current location - list_of_layers = given_list_of_layers; - current_child = children{child}; - list_of_layers{end+1} = current_child; - - % this is the value of the current location in the model (ie, md.mesh.elements) - location_of_child = model_subclass.(current_child); - - % if the empty model does not have this attribute, it's because it's new so we save it to netcdf - % there are 2 cases: the location is a struct, the location is a class - if isstruct(model_subclass) - % if the current field is a nested struct assume it has valuable data that needs to be saved - if isstruct(location_of_child) && any(size(location_of_child) > 1) - create_group(location_of_child, list_of_layers, NetCDF, verbose); - - % this would mean that the layer above the layer we're interested in is a struct, so - % we can navigate our empty model as such - elseif isfield(empty_model_subclass, current_child) - % the layer we're interested in does exist, we just need to compare states - location_of_child_in_empty_model = empty_model_subclass.(current_child); - - % if the current attribute is a numerical array assume it has valuable data that needs to be saved - if isnumeric(location_of_child) && logical(numel(location_of_child) > 1) - create_group(location_of_child, list_of_layers, NetCDF, verbose); - % if the attributes are identical we don't need to save anything - elseif (all(isnan(location_of_child)) && all(isnan(location_of_child_in_empty_model))) || isempty(setxor(location_of_child, location_of_child_in_empty_model)) - walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); - % if the attributes are not the same we need to save ours - else - % THE ORDER OF THESE LINES IS CRITICAL - walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - % this would mean that the layer we're interested in is not fundamental to the model architecture - % and thus needs to be saved to the netcdf - else - walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - % this would mean it's not a struct, and must be a class/subclass - % we now check the state of the class property - else - try - if isprop(empty_model_subclass, current_child) - % the layer we're interested in does exist, we just need to compare states - location_of_child_in_empty_model = empty_model_subclass.(current_child); - % if the current attribute is a numerical array assume it has valuable data that needs to be saved - if isnumeric(location_of_child) && logical(numel(location_of_child) > 1) - create_group(location_of_child, list_of_layers, NetCDF, verbose); - - elseif iscell(location_of_child) - % if the attributes are identical we don't need to save anything - if isempty(setxor(location_of_child, location_of_child_in_empty_model)) - % pass - else - % otherwise we need to save - walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - elseif (all(isnan(location_of_child)) && all(isnan(location_of_child_in_empty_model))) - walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); - % if the attributes are not the same we need to save ours - else - % THE ORDER OF THESE LINES IS CRITICAL - walk_through_subclasses(location_of_child, location_of_child_in_empty_model, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - else - walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - catch - walk_through_subclasses(location_of_child, empty_model_subclass, list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - end - end - end - catch ME - % If the caught error is a fieldname error, it's just saying that a variable has no fields and thus can be ignored - if strcmp(ME.identifier, 'MATLAB:fieldnames:InvalidInput') - % do nothing - % this is if we come accross instances/subfields in our model that are not fundamental to the model class (ie, taoinversion) - elseif strcmp(ME.identifier, 'MATLAB:UndefinedFunction') - walk_through_subclasses(location_of_child, empty_model_subclass, given_list_of_layers, empty_model, NetCDF, verbose); - create_group(location_of_child, list_of_layers, NetCDF, verbose); - % If it's a different error, rethrow it to MATLAB's default error handling - else - disp(ME.identifier) - disp(given_list_of_layers) - rethrow(ME); - end - end -end - - -function create_group(location_of_child, list_of_layers, NetCDF, verbose) - %disp(list_of_layers) - % location_of_child is an object - % list_of_layers is a list like {'inversion', 'StressbalanceSolution','cost_functions_coefficients'} - % first we make the group at the highest level (ie, inversion) - group_name = list_of_layers{1}; - variable_name = list_of_layers{end}; - - % if the group is already made, get it's ID instead of creating it again - try % group hasn't been made - group = netcdf.defGrp(NetCDF, group_name); - catch % group was already made - group = netcdf.inqNcid(NetCDF, group_name); - end - - % if the data is nested, create nested groups to match class structure - if numel(list_of_layers) > 2 - % the string() method is really important here since matlab apparently can't handle the infinite complexity of a string without the string method. - for name = string(list_of_layers(2:end-1)) - % the group levels may have already been made - try % group hasn't been made - group = netcdf.defGrp(group, name); - catch % group was already made - group = netcdf.inqNcid(group, name); - end - end - end - % sometimes objects are passed through twice so we account for that with this try/catch + + +% ===================================================================== +% write_group – write all fields of a class/struct into a group +% ===================================================================== +function write_group(obj, gid, root_ncid, path_str, verbose) + % Get field names (works for both classes and structs) try - % we may be dealing with an object - % first we screen for structs - if isstruct(location_of_child) % && any(size(location_of_child) > 1) -- this is being tested - % we have a struct - copy_nested_struct(variable_name, location_of_child, group, NetCDF, verbose); - - % now for cell arrays of datastructures: - elseif logical(~isstruct(location_of_child) && iscell(location_of_child) && isobject(location_of_child{1})) - copy_cell_array_of_objects(variable_name, location_of_child, group, NetCDF, verbose); - else - if ~isobject(location_of_child) && ~isstruct(location_of_child) - % we're dealing with raw data - create_var(variable_name, location_of_child, group, NetCDF, verbose); - end - end + fields = fieldnames(obj); catch - % do nothing + % obj has no fields (e.g. a raw value that slipped through) – ignore + return end -end + for fi = 1:numel(fields) + fname = fields{fi}; + val = obj.(fname); + child_path = [path_str '.' fname]; + + % Skip log fields + if ismember(fname, {'errlog','outlog'}) + continue + end + write_value(val, fname, gid, root_ncid, child_path, verbose); + end +end -function copy_cell_array_of_objects(variable_name, address_of_child, group, NetCDF, verbose) - % make subgroup to represent the array - [rows, cols] = size(address_of_child); - name_of_subgroup = [num2str(rows), 'x', num2str(cols), '_cell_array_of_objects']; - subgroup = netcdf.defGrp(group, name_of_subgroup); - % save the name of the cell array - write_string_to_netcdf('name_of_cell_array', variable_name, subgroup, NetCDF, verbose); +% ===================================================================== +% write_value – dispatch a single value to the right writer +% ===================================================================== +function write_value(val, vname, gid, root_ncid, path_str, verbose) - % save the dimensions of the cell array - create_var('rows', rows, subgroup, NetCDF, verbose); - create_var('cols', cols, subgroup, NetCDF, verbose); + % ---- empty ------------------------------------------------------- + if isempty(val) + % Store as scalar NaN so the variable exists and readers know it was empty + write_scalar_nan(vname, gid, root_ncid); + if verbose, fprintf(' [empty] %s\n', path_str); end + return + end - % if this is a multidimensional cell array, iterate over rows here and cols in copy_objects - if rows>1 - for row = 1:rows - % make a subgroup for each row - name_of_subgroup = ['Row_', num2str(row), '_of_', num2str(rows)]; - subgroup = netcdf.defGrp(group, name_of_subgroup); - copy_objects(address_of_child, subgroup, NetCDF, cols, verbose); - end - else - copy_objects(address_of_child, subgroup, NetCDF, cols, verbose); + % ---- struct (e.g. md.results.TransientSolution which is 1xN) ----- + if isstruct(val) + write_struct_value(val, vname, gid, root_ncid, path_str, verbose); + return end -end + % ---- ISSM sub-class (has its own fieldnames) ---------------------- + if isobject(val) + sub_gid = get_or_create_group(gid, vname); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'classtype', class(val)); + write_group(val, sub_gid, root_ncid, path_str, verbose); + return + end + % ---- cell array of ISSM objects (e.g. outputdefinition list) ------ + if iscell(val) && ~isempty(val) && isobject(val{1}) + write_cell_of_objects(val, vname, gid, root_ncid, path_str, verbose); + return + end -function copy_objects(address_of_child, group, NetCDF, cols, verbose) - for col = 1:cols - % make subgroup to contain each col of array - name_of_subgroup = ['Col_', num2str(col), '_of_', num2str(cols)]; - subgroup = netcdf.defGrp(group, name_of_subgroup); + % ---- cell array of strings --------------------------------------- + if iscellstr(val) || (iscell(val) && ~isempty(val) && ischar(val{1})) + write_cell_strings(val, vname, gid, root_ncid, path_str, verbose); + if verbose, fprintf(' [cellstr] %s\n', path_str); end + return + end - % get the kind of object we're working with: - if isstruct(address_of_child{col}) - % handle structs - name_raw = fields(address_of_child{col}); - variable_name = name_raw{1}; - copy_nested_struct(variable_name, address_of_child, subgroup, NetCDF, verbose); - - elseif numel(properties(address_of_child{col})) > 0 - % handle class instances - copy_class_instance(address_of_child{col}, subgroup, NetCDF, verbose); - else - disp('ERROR: Cell arrays of mixed types are not yet supported in read_netCDF!\n Deserialization will not be able to complete!') - % handle regular datastructures that are already supported - name_raw = fields(address_of_child); - variable_name = name_raw{col}; - create_var(variable_name, address_of_child, subgroup, NetCDF, verbose); + % ---- plain cell array (numbers etc.) ----------------------------- + if iscell(val) + % Try to convert to numeric matrix; if not, skip with warning + try + mat = cell2mat(val); + write_numeric(mat, vname, gid, root_ncid, path_str, verbose); + catch + if verbose + fprintf(' [SKIP] %s – cell array of mixed/unsupported types\n', path_str); + end end + return end -end - -function copy_class_instance(address_of_child, subgroup, NetCDF, verbose) - % get parent class name - name = class(address_of_child); - - % save the name of the class - write_string_to_netcdf('class_is_a', name, subgroup, NetCDF, verbose); - - % make subgroup to contain properties - name_of_subgroup = ['Properties_of_', name]; - subgroup = netcdf.defGrp(subgroup, name_of_subgroup); + % ---- logical (bool) ---------------------------------------------- + if islogical(val) + write_bool(val, vname, gid, root_ncid, path_str, verbose); + if verbose, fprintf(' [bool] %s\n', path_str); end + return + end - % get properties - props = properties(address_of_child); + % ---- char / string ----------------------------------------------- + if ischar(val) || isstring(val) + write_string(char(val), vname, gid, root_ncid, path_str, verbose); + if verbose, fprintf(' [string] %s\n', path_str); end + return + end - for property = 1:length(props) - variable_name = props{property}; - create_var(variable_name, address_of_child.(variable_name), subgroup, NetCDF, verbose); + % ---- numeric (scalar, vector, matrix, ND) ------------------------ + if isnumeric(val) + write_numeric(val, vname, gid, root_ncid, path_str, verbose); + if verbose, fprintf(' [numeric] %s %s\n', path_str, mat2str(size(val))); end + return end + % ---- fallback: skip with warning --------------------------------- + if verbose + fprintf(' [SKIP] %s – unsupported type %s\n', path_str, class(val)); + end end -function copy_nested_struct(parent_struct_name, address_of_struct, group, NetCDF, verbose) - %{ - This function takes a struct of structs and saves them to netcdf. - - It also works with single structs. - - To do this, we get the number of dimensions (substructs) of the parent struct. - Next, we iterate through each substruct and record the data. - For each substruct, we create a subgroup of the main struct. - For each variable, we create dimensions that are assigned to each subgroup uniquely. - %} - - % make a new subgroup to contain all the others: - group = netcdf.defGrp(group, parent_struct_name); - - % make sure other systems can flag the nested struct type - dimID = netcdf.defDim(group, 'struct', 6); - string_var = netcdf.defVar(group, 'this_is_a_nested', "NC_CHAR", dimID); - uint_method=uint8('struct').'; - method_ID = char(uint_method); - netcdf.putVar(group, string_var, method_ID); - - % other systems know the name of the parent struct because it's covered by the results/qmu functions above - - % 'a' will always be 1 and is not useful to us - [a, no_of_dims] = size(address_of_struct); - - for substruct = 1:no_of_dims - % we start by making subgroups with nice names like "TransientSolution_substruct_44" - name_of_subgroup = ['1x', num2str(substruct)]; - subgroup = netcdf.defGrp(group, name_of_subgroup); - - % do some housekeeping to keep track of the current layer - current_substruct = address_of_struct(substruct); - substruct_fields = fieldnames(current_substruct)'; % transpose because matlab only interates over n x 1 arrays - - % now we need to iterate over each variable of the nested struct and save it to this new subgroup - for variable_name = string(substruct_fields) - address_of_child = current_substruct.(variable_name); - create_var(variable_name, address_of_child, subgroup, NetCDF, verbose); +% ===================================================================== +% write_struct_value – handle a struct field (scalar or 1-D array) +% ===================================================================== +function write_struct_value(val, vname, gid, root_ncid, path_str, verbose) + n = numel(val); + if n == 0 + write_scalar_nan(vname, gid, root_ncid); + return + end + + sub_gid = get_or_create_group(gid, vname); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'classtype', 'struct'); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'nsteps', int32(n)); + + for k = 1:n + step_name = sprintf('step_%d', k); + step_gid = netcdf.defGrp(sub_gid, step_name); + step_fields = fieldnames(val(k)); + for fi = 1:numel(step_fields) + sfield = step_fields{fi}; + if ismember(sfield, {'errlog','outlog'}), continue; end + sval = val(k).(sfield); + write_value(sval, sfield, step_gid, root_ncid, [path_str '.' step_name '.' sfield], verbose); end end if verbose - fprintf(["Succesfully transferred nested MATLAB struct ", parent_struct_name, " to the NetCDF\n"]) + fprintf(' [struct] %s (1x%d)\n', path_str, n); end end - -% ironically inversion does not have the same problem as results as inversion subfields -% are actually subclasses and not fields -function check_inversion_class(model_var, NetCDF, verbose) - - % Define a persistent variable to ensure this function is only run once - persistent executed; - % Check if the function has already been executed - if isempty(executed) - if verbose - disp('Deconstructing Inversion class instance') - end - % Need to make sure that we have the right inversion class: inversion, m1qn3inversion, taoinversion - groupid = netcdf.inqNcid(NetCDF,'inversion'); - - if isa(model_var, 'm1qn3inversion') - write_string_to_netcdf('inversion_class_name', 'm1qn3inversion', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved inversion class instance m1qn3inversion') - end - elseif isa(model_var, 'taoinversion') - write_string_to_netcdf('inversion_class_name', 'taoinversion', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved inversion class instance taoinversion') - end - else - write_string_to_netcdf('inversion_class_name', 'inversion', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved inversion class instance inversion') +% ===================================================================== +% write_cell_of_objects – cell array whose elements are ISSM classes +% ===================================================================== +function write_cell_of_objects(val, vname, gid, root_ncid, path_str, verbose) + [rows, cols] = size(val); + sub_gid = get_or_create_group(gid, vname); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'classtype', 'cell_of_objects'); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'nrows', int32(rows)); + netcdf.putAtt(sub_gid, netcdf.getConstant('NC_GLOBAL'), 'ncols', int32(cols)); + + for r = 1:rows + for c = 1:cols + elem = val{r, c}; + elem_name = sprintf('item_%d_%d', r, c); + elem_gid = netcdf.defGrp(sub_gid, elem_name); + netcdf.putAtt(elem_gid, netcdf.getConstant('NC_GLOBAL'), 'classtype', class(elem)); + if isstruct(elem) + write_group(elem, elem_gid, root_ncid, [path_str '.' elem_name], verbose); + elseif isobject(elem) + write_group(elem, elem_gid, root_ncid, [path_str '.' elem_name], verbose); end end - % Set the persistent variable to indicate that the function has been executed - executed = true; + end + if verbose + fprintf(' [cell] %s (%dx%d objects)\n', path_str, rows, cols); end end -%Write a name for the smb class -function check_smb_class(model_var, NetCDF, verbose) -% Define a persistent variable to ensure this function is only run once -persistent executed; -% Check if the function has already been executed -if isempty(executed) - if verbose - disp('Deconstructing smb class instance') - end - % Need to make sure that we have the right smb class - groupid = netcdf.inqNcid(NetCDF,'smb'); +% ===================================================================== +% Leaf writers +% ===================================================================== - if isa(model_var, 'SMBforcing') - write_string_to_netcdf('smb_class_name', 'SMBforcing', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBforcing') - end - elseif isa(model_var, 'SMB') - write_string_to_netcdf('smb_class_name', 'SMB', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMB') - end - elseif isa(model_var, 'SMBcomponents') - write_string_to_netcdf('smb_class_name', 'SMBcomponents', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBcomponents') - end - elseif isa(model_var, 'SMBd18opdd') - write_string_to_netcdf('smb_class_name', 'SMBd18opdd', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBd18opdd') - end - elseif isa(model_var, 'SMBdebrisEvatt') - write_string_to_netcdf('smb_class_name', 'SMBdebrisEvatt', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBdebrisEvatt') - end - elseif isa(model_var, 'SMBgemb') - write_string_to_netcdf('smb_class_name', 'SMBgemb', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBgemb') - end - elseif isa(model_var, 'SMBgradients') - write_string_to_netcdf('smb_class_name', 'SMBgradients', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBgradients') - end - elseif isa(model_var, 'SMBgradientscomponents') - write_string_to_netcdf('smb_class_name', 'SMBgradientscomponents', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBgradientscomponents') - end - elseif isa(model_var, 'SMBgradientsela') - write_string_to_netcdf('smb_class_name', 'SMBgradientsela', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBgradientsela') - end - elseif isa(model_var, 'SMBhenning') - write_string_to_netcdf('smb_class_name', 'SMBhenning', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBhenning') - end - elseif isa(model_var, 'SMBmeltcomponents') - write_string_to_netcdf('smb_class_name', 'SMBmeltcomponents', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBmeltcomponents') - end - elseif isa(model_var, 'SMBpdd') - write_string_to_netcdf('smb_class_name', 'SMBpdd', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBpdd') - end - elseif isa(model_var, 'SMBpddSicopolis') - write_string_to_netcdf('smb_class_name', 'SMBpddSicopolis', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBpddSicopolis') - end - elseif isa(model_var, 'SMBsemic') - write_string_to_netcdf('smb_class_name', 'SMBsemic', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved smb class instance SMBsemic') - end +function write_numeric(val, vname, gid, root_ncid, path_str, verbose) + % MATLAB stores data in column-major order; NetCDF uses row-major. + % We transpose so Python reads the same shape as MATLAB. + if isempty(val) + write_scalar_nan(vname, gid, root_ncid); + return end - % Set the persistent variable to indicate that the function has been executed - executed = true; -end -end + data = double(val); -function check_friction_class(model_var, NetCDF, verbose) -% Define a persistent variable to ensure this function is only run once -persistent executed; -% Check if the function has already been executed -if isempty(executed) - if verbose - disp('Deconstructing friction class instance') - end - % Need to make sure that we have the right friction class - groupid = netcdf.inqNcid(NetCDF,'friction'); - if isa(model_var, 'friction') - write_string_to_netcdf('friction_class_name', 'friction', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance friction') - end - elseif isa(model_var, 'frictioncoulomb') - write_string_to_netcdf('friction_class_name', 'frictioncoulomb', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictioncoulomb') - end - elseif isa(model_var, 'frictioncoulomb2') - write_string_to_netcdf('friction_class_name', 'frictioncoulomb2', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictioncoulomb2') - end - elseif isa(model_var, 'frictionhydro') - write_string_to_netcdf('friction_class_name', 'frictionhydro', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionhydro') - end - elseif isa(model_var, 'frictionjosh') - write_string_to_netcdf('friction_class_name', 'frictionjosh', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionjosh') - end - elseif isa(model_var, 'frictionpism') - write_string_to_netcdf('friction_class_name', 'frictionpism', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionpism') - end - elseif isa(model_var, 'frictionregcoulomb') - write_string_to_netcdf('friction_class_name', 'frictionregcoulomb', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionregcoulomb') - end - elseif isa(model_var, 'frictionregcoulomb2') - write_string_to_netcdf('friction_class_name', 'frictionregcoulomb2', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionregcoulomb2') - end - elseif isa(model_var, 'frictionschoof') - write_string_to_netcdf('friction_class_name', 'frictionschoof', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionschoof') - end - elseif isa(model_var, 'frictionshakti') - write_string_to_netcdf('friction_class_name', 'frictionshakti', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionshakti') - end - elseif isa(model_var, 'frictiontsai') - write_string_to_netcdf('friction_class_name', 'frictiontsai', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictiontsai') - end - elseif isa(model_var, 'frictionwaterlayer') - write_string_to_netcdf('friction_class_name', 'frictionwaterlayer', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionwaterlayer') - end - elseif isa(model_var, 'frictionweertman') - write_string_to_netcdf('friction_class_name', 'frictionweertman', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionweertman') - end - elseif isa(model_var, 'frictionweertmantemp') - write_string_to_netcdf('friction_class_name', 'frictionweertmantemp', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved friction class instance frictionweertmantemp') - end + % Flatten trailing singleton dimensions (keep at least 2 dims for reshape) + sz = size(data); + while numel(sz) > 2 && sz(end) == 1 + sz = sz(1:end-1); end - % Set the persistent variable to indicate that the function has been executed - executed = true; -end -end + data = reshape(data, sz); -function check_hydrology_class(model_var, NetCDF, verbose) -% Define a persistent variable to ensure this function is only run once -persistent executed; -% Check if the function has already been executed -if isempty(executed) - if verbose - disp('Deconstructing hydrology class instance') + % Transpose for row-major storage (only for 2D) + if ndims(data) == 2 && ~any(sz == 1) + data = data.'; + sz = size(data); + elseif any(sz == 1) && numel(sz) == 2 + % Vector: store as 1-D + data = data(:); + sz = numel(data); end - % Need to make sure that we have the right hydrology class - groupid = netcdf.inqNcid(NetCDF,'hydrology'); - if isa(model_var, 'hydrologyshreve') - write_string_to_netcdf('hydrology_class_name', 'hydrologyshreve', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologyshreve') - end - elseif isa(model_var, 'hydrology') - write_string_to_netcdf('hydrology_class_name', 'hydrology', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrology') - end - elseif isa(model_var, 'hydrologyarmapw') - write_string_to_netcdf('hydrology_class_name', 'hydrologyarmapw', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologyarmapw') - end - elseif isa(model_var, 'hydrologydc') - write_string_to_netcdf('hydrology_class_name', 'hydrologydc', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologydc') - end - elseif isa(model_var, 'hydrologyglads') - write_string_to_netcdf('hydrology_class_name', 'hydrologyglads', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologyglads') - end - elseif isa(model_var, 'hydrologypism') - write_string_to_netcdf('hydrology_class_name', 'hydrologypism', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologypism') - end - elseif isa(model_var, 'hydrologyshakti') - write_string_to_netcdf('hydrology_class_name', 'hydrologyshakti', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologyshakti') - end - elseif isa(model_var, 'hydrologytws') - write_string_to_netcdf('hydrology_class_name', 'hydrologytws', groupid, NetCDF, verbose); - if verbose - disp('Successfully saved hydrology class instance hydrologytws') - end + dim_ids = make_dims(sz, gid, root_ncid); + + varid = netcdf.defVar(gid, vname, 'NC_DOUBLE', dim_ids); + if numel(data) > 1 + netcdf.defVarDeflate(gid, varid, true, true, 4); end - % Set the persistent variable to indicate that the function has been executed - executed = true; -end + netcdf.putVar(gid, varid, data); end -function create_var(variable_name, address_of_child, group, NetCDF, verbose) - % There are lots of different variable types that we need to handle from the model class - - % get the dimensions we'll need - intdim = netcdf.inqDimID(NetCDF,'int'); - floatdim = netcdf.inqDimID(NetCDF,'float'); - unlimdim = netcdf.inqDimID(NetCDF,'Unlim'); - - % This first conditional statement will catch numeric arrays (matrices) of any dimension and save them - if any(size(address_of_child)>1) && ~iscellstr(address_of_child) && ~ischar(address_of_child) - write_numeric_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose); - - % check if it's a string - elseif ischar(address_of_child) - write_string_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose); - - % or an empty variable - elseif isempty(address_of_child) - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", intdim); - - % or a list of strings - elseif iscellstr(address_of_child) || iscell(address_of_child) && ischar(address_of_child{1}) - write_cell_with_strings(variable_name, address_of_child, group, NetCDF, verbose) - - % or an empty list - elseif iscell(address_of_child) && isempty(address_of_child) || isa(address_of_child, 'double') && isempty(address_of_child) - variable = netcdf.defVar(group, variable_name, "NC_INT", intdim); - netcdf.putVar(group,variable, -32767); - - % or a bool - elseif islogical(address_of_child) - % netcdf4 can't handle bool types like true/false so we convert all to int 1/0 and add an attribute named units with value 'bool' - variable = netcdf.defVar(group, variable_name, 'NC_SHORT', intdim); - netcdf.putVar(group,variable,int8(address_of_child)); - % make sure other systems can flag the bool type - netcdf.putAtt(group,variable,'units','bool'); - - % or a regular list - elseif iscell(address_of_child) - disp('made list w/ unlim dim') - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); - netcdf.putVar(group,variable,address_of_child); - - % or a float - elseif isfloat(address_of_child) && numel(address_of_child) == 1 - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", floatdim); - netcdf.putVar(group,variable,address_of_child); - - % or a int - elseif mod(address_of_child,1) == 0 || isinteger(address_of_child) && numel(address_of_child) == 1 - variable = netcdf.defVar(group, variable_name, "NC_SHORT", intdim); - netcdf.putVar(group,variable,address_of_child); - - % anything else... (will likely need to add more cases; ie dict) - else + +function write_string(str, vname, gid, root_ncid, path_str, verbose) + if isempty(str) + % Zero-length string: store empty char variable + dim_name = 'char0'; try - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); - netcdf.putVar(group,variable,address_of_child); - catch ME - disp(ME.message); - disp(['Datatype given: ', class(address_of_child)]); + dimid = netcdf.defDim(gid, dim_name, 0); + catch + dimid = netcdf.inqDimID(gid, dim_name); end + netcdf.defVar(gid, vname, 'NC_CHAR', dimid); + return end - if verbose - fprintf('Successfully transferred data from %s to the NetCDF\n', variable_name); + n = numel(str); + dim_name = sprintf('char%d', n); + try + dimid = netcdf.defDim(gid, dim_name, n); + catch + dimid = netcdf.inqDimID(gid, dim_name); end + varid = netcdf.defVar(gid, vname, 'NC_CHAR', dimid); + netcdf.putVar(gid, varid, str); + netcdf.putAtt(gid, varid, 'type_is', 'string'); end -function write_cell_with_strings(variable_name, address_of_child, group, NetCDF, verbose) - %{ - Write cell array (ie {'one' 'two' 'three'}) to netcdf - %} - - if isempty(address_of_child) - % if the char array is empty, save an empty char - name_of_dimension = ['char', num2str(0)]; +function write_cell_strings(val, vname, gid, root_ncid, path_str, verbose) + if isempty(val) + dim_name = 'char0'; try - dimID = netcdf.defDim(group, name_of_dimension, 0); + dimid = netcdf.defDim(gid, dim_name, 0); catch - dimID = netcdf.inqDimID(group, name_of_dimension); + dimid = netcdf.inqDimID(gid, dim_name); end - % Now we can make a variable in this dimension: - string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); - % we leave empty now - else - % covert data to char array - method_ID = char(address_of_child); - - % make dimensions - [rows, cols] = size(method_ID); - - IDDim1 = netcdf.defDim(group,'cols',cols); - IDDim2 = netcdf.defDim(group,'rows',rows); - - % create the variable slot - IDVarId = netcdf.defVar(group,variable_name,'NC_CHAR', [IDDim1 IDDim2]); - - % save the variable - netcdf.putVar(group, IDVarId, method_ID'); %transpose - - % tell other platforms that this is a cell of strings - netcdf.putAtt(group, IDVarId, 'type_is','cell_array_of_strings'); + varid = netcdf.defVar(gid, vname, 'NC_CHAR', dimid); + netcdf.putAtt(gid, varid, 'type_is', 'cell_of_strings'); + return end -end + % Pad all strings to the same length, store as NC_CHAR rows×cols + strs = cellstr(val); + nrows = numel(strs); + max_len = max(cellfun(@numel, strs)); + if max_len == 0, max_len = 1; end + char_mat = repmat(' ', nrows, max_len); + for k = 1:nrows + s = strs{k}; + char_mat(k, 1:numel(s)) = s; + end -function write_string_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose) - % netcdf and strings don't get along.. we have to do it 'custom': + % NC dims: [cols rows] because netcdf.putVar transposes + try + dim_cols = netcdf.defDim(gid, sprintf('strcols_%s', vname), max_len); + catch + dim_cols = netcdf.inqDimID(gid, sprintf('strcols_%s', vname)); + end + try + dim_rows = netcdf.defDim(gid, sprintf('strrows_%s', vname), nrows); + catch + dim_rows = netcdf.inqDimID(gid, sprintf('strrows_%s', vname)); + end - the_string_to_save = address_of_child; + varid = netcdf.defVar(gid, vname, 'NC_CHAR', [dim_cols dim_rows]); + netcdf.putVar(gid, varid, char_mat.'); % transpose for MATLAB→NetCDF + netcdf.putAtt(gid, varid, 'type_is', 'cell_of_strings'); +end - if isempty(the_string_to_save) - % if the char array is empty, save an empty char - name_of_dimension = ['char', num2str(0)]; - try - dimID = netcdf.defDim(group, name_of_dimension, 0); - catch - dimID = netcdf.inqDimID(group, name_of_dimension); - end - % Now we can make a variable in this dimension: - string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); - % we leave empty now + +function write_bool(val, vname, gid, root_ncid, path_str, verbose) + if isscalar(val) + scalardim = netcdf.inqDimID(root_ncid, 'scalar'); + varid = netcdf.defVar(gid, vname, 'NC_SHORT', scalardim); + netcdf.putVar(gid, varid, int16(val)); else - % convert string to - uint_method=uint8(the_string_to_save).'; - method_ID = char(uint_method); - length_of_the_string = numel(method_ID); - - % Convert the string to character data using string array - %str_out = char(the_string_to_save) - - % Determine the length of the string - %length_of_the_string = numel(str_out) - - % Check if the dimension already exists, and if not, create it - name_of_dimension = ['char', num2str(length_of_the_string)]; + data = int16(val(:)); + n = numel(data); try - dimID = netcdf.defDim(group, name_of_dimension, length_of_the_string); + dimid = netcdf.defDim(gid, sprintf('booldim%d', n), n); catch - dimID = netcdf.inqDimID(group, name_of_dimension); + dimid = netcdf.inqDimID(gid, sprintf('booldim%d', n)); end - % Now we can make a variable in this dimension: - string_var = netcdf.defVar(group, variable_name, "NC_CHAR", [dimID]); - % Finally, we can write the variable (always transpose for matlab): - netcdf.putVar(group, string_var, method_ID); + varid = netcdf.defVar(gid, vname, 'NC_SHORT', dimid); + netcdf.putVar(gid, varid, data); end + netcdf.putAtt(gid, varid, 'type_is', 'bool'); +end - if verbose - disp(['Successfully transferred data from ', variable_name, ' to the NetCDF']); - end + +function write_scalar_nan(vname, gid, root_ncid) + scalardim = netcdf.inqDimID(root_ncid, 'scalar'); + varid = netcdf.defVar(gid, vname, 'NC_DOUBLE', scalardim); + netcdf.putVar(gid, varid, NaN); + netcdf.putAtt(gid, varid, 'type_is', 'empty'); end -function write_numeric_array_to_netcdf(variable_name, address_of_child, group, NetCDF, verbose) - - % get the dimensions we'll need - intdim = netcdf.inqDimID(NetCDF,'int'); - floatdim = netcdf.inqDimID(NetCDF,'float'); - unlimdim = netcdf.inqDimID(NetCDF,'Unlim'); - - typeis = class(address_of_child); - - if isa(typeis, 'logical') - % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat - % we need to transpose data before it goes into netcdf - data = address_of_child.'; - - % make the dimensions - dimensions = []; - for dimension = size(data) - dim_name = ['dim',int2str(dimension)]; - % if the dimension already exists we can't have a duplicate +% ===================================================================== +% Helpers +% ===================================================================== + +function dim_ids = make_dims(sz, gid, root_ncid) + % Create or reuse dimensions for an array of size sz. + % sz must be a row vector of positive integers. + dim_ids = zeros(1, numel(sz)); + for k = 1:numel(sz) + n = sz(k); + if n == 1 + dim_ids(k) = netcdf.inqDimID(root_ncid, 'scalar'); + else + dim_name = sprintf('dim%d', n); + % Try in the local group first, then root + try + dim_ids(k) = netcdf.defDim(gid, dim_name, n); + catch try - dimID = netcdf.defDim(group, dim_name, dimension); + dim_ids(k) = netcdf.inqDimID(gid, dim_name); catch - dimID = netcdf.inqDimID(group, dim_name); - end - % record the dimension for the variable - dimensions(end+1) = dimID; - end - - % write the variable - netcdf.putVar(group,variable,data); - - % make sure other systems can flag the bool type - netcdf.putAtt(group,variable,'units','bool'); - - % handle all other datatypes here - else - % sometimes an array has just 1 element in it, we account for those cases here: - if numel(address_of_child) == 1 - if isinteger(address_of_child) - variable = netcdf.defVar(group, variable_name, "NC_SHORT", intdim); - netcdf.putVar(group,variable,address_of_child); - elseif isa(address_of_child, 'double') || isa(address_of_child, 'float') - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", floatdim); - netcdf.putVar(group,variable,address_of_child); - else - disp('Encountered single datatype that was not float64 or int64, saving under unlimited dimension, may cause errors.') - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE", unlimdim); - netcdf.putVar(group,variable,address_of_child); - end - % this is in case of lists so that python doesn't get a (nx1) numpy array and instead gets an n-element list - elseif any(size(address_of_child)==1) - % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat - % we need to transpose data before it goes into netcdf - data = address_of_child.'; - - % make the dimensions - dimensions = []; - for dimension = size(data) - if dimension ~= 1 - dim_name = ['dim',int2str(dimension)]; - % if the dimension already exists we can't have a duplicate try - dimID = netcdf.defDim(group, dim_name, dimension); + dim_ids(k) = netcdf.defDim(root_ncid, dim_name, n); catch - dimID = netcdf.inqDimID(group, dim_name); + dim_ids(k) = netcdf.inqDimID(root_ncid, dim_name); end - % record the dimension for the variable - dimensions(end+1) = dimID; - end - end - % create the variable - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE",dimensions); - - % write the variable - netcdf.putVar(group,variable,data); - - % This catches all remaining arrays: - else - % because matlab transposes all data into and out of netcdf and because we want cross-platform-compat - % we need to transpose data before it goes into netcdf - data = address_of_child.'; - - % make the dimensions - dimensions = []; - for dimension = size(data) - dim_name = ['dim',int2str(dimension)]; - % if the dimension already exists we can't have a duplicate - try - dimID = netcdf.defDim(group, dim_name, dimension); - catch - dimID = netcdf.inqDimID(group, dim_name); end - % record the dimension for the variable - dimensions(end+1) = dimID; end - % create the variable - variable = netcdf.defVar(group, variable_name, "NC_DOUBLE",dimensions); - - % write the variable - netcdf.putVar(group,variable,data); end end end + + +function gid = get_or_create_group(parent_gid, name) + try + gid = netcdf.defGrp(parent_gid, name); + catch + gid = netcdf.inqNcid(parent_gid, name); + end +end diff --git a/src/m/netcdf/write_netCDF.py b/src/m/netcdf/write_netCDF.py index f3817e0c2..25c0d8552 100644 --- a/src/m/netcdf/write_netCDF.py +++ b/src/m/netcdf/write_netCDF.py @@ -1,595 +1,405 @@ -# imports -import netCDF4 -from netCDF4 import Dataset -import numpy as np -import numpy.ma as ma -import time +""" +write_netCDF – Save an ISSM model to a NetCDF4 file. + +Usage +----- + from write_netCDF import write_netCDF + write_netCDF(md, 'model.nc') + write_netCDF(md, 'model.nc', verbose=True) + +The file can be read back with read_netCDF.py (Python) or read_netCDF.m (MATLAB). + +NetCDF layout +------------- +Each attribute of md (mesh, mask, geometry, …) becomes a NetCDF group. +Attributes within those objects are stored as variables inside the group. +Nested sub-classes get nested sub-groups. + +Every group carries a 'classtype' global attribute with the Python class name +so that read_netCDF can reconstruct the correct object on load. + +Results storage +--------------- +md.results. is a solution() instance whose .steps list holds +solutionstep() objects. Each step is stored as a sub-group named 'step_' +inside the solution sub-group. +""" + import os -from model import * -from results import * -from m1qn3inversion import m1qn3inversion -from taoinversion import taoinversion -#import OrderedStruct +import time +import numpy as np +from netCDF4 import Dataset + + +def write_netCDF(md, filename: str, verbose: bool = False) -> None: + """Save model *md* to *filename* (NetCDF4). Overwrites silently.""" + if verbose: + print('write_netCDF v2.0 (Python)') + + # Overwrite silently + if os.path.exists(filename): + os.remove(filename) + if verbose: + print(f'Existing file {filename} deleted.') + + nc = Dataset(filename, 'w', format='NETCDF4') + try: + nc.history = 'Created by write_netCDF.py on ' + time.ctime() + nc.Conventions = 'ISSM-NetCDF-2.0' + + # Global dimensions reused throughout the file + nc.createDimension('scalar', 1) + nc.createDimension('Unlim', None) # unlimited + + # Walk every top-level attribute of the model + for gname, val in md.__dict__.items(): + if val is None: + if verbose: + print(f' md.{gname} is None – skipped') + continue + grp = nc.createGroup(gname) + grp.classtype = type(val).__name__ -''' -Given a md, this set of functions will perform the following: - 1. View each attribute of each nested class. - 2. Compare state of attribute in the model to an empty model. - 3. If states are identical, pass. (except for np arrays which will always be saved) - 4. Otherwise, create nested groups named after class structure. - 5. Create variable named after class attribute and assign value to it. -''' + _write_obj(val, grp, nc, f'md.{gname}', verbose) + except Exception: + nc.close() + raise -def write_netCDF(md, filename: str, verbose = False): + nc.close() if verbose: - print('Python C2NetCDF4 v1.2.0') - else: pass - ''' - md = model class instance to be saved - filename = path and name to save file under - verbose = T/F = show or muted log statements. Naturally muted - ''' - # this is a precaution so that data is not lost + print('Model successfully saved as NetCDF4.') + + +# ===================================================================== +# _write_obj – write all attributes of a Python object into a group +# ===================================================================== +def _write_obj(obj, grp, root_nc, path: str, verbose: bool) -> None: + """Iterate over obj.__dict__ and write each attribute into *grp*.""" try: - # Create a NCData file to write to - NCData = create_NetCDF(filename, verbose) - - # Create an instance of an empty md class to compare md_var against - empty_model = model() - - # Walk through the md class and compare subclass states to empty_model - walk_through_model(md, empty_model, NCData, verbose) - - # in order to handle some subclasses in the results class, we have to utilize this band-aid - # there will likely be more band-aids added unless a class name library is created with all class names that might be added to a md - try: - # if results has meaningful data, save the name of the subclass and class instance - NCData.groups['results'] - results_subclasses_bandaid(md, NCData, verbose) - # otherwise, ignore - except KeyError: - pass - - NCData.close() - if verbose: - print('Model successfully saved as NetCDF4') - else: pass - - # just in case something unexpected happens - except Exception as e: - if 'NCData' in locals(): - NCData.close() - raise e - - -def results_subclasses_bandaid(md, NCData, verbose = False): - # since the results class may have nested classes within it, we need to record the name of the - # nested class instance variable as it appears in the md that we're trying to save - quality_control = [] - - # we save lists of instances to the NCData - solutions = [] - solutionsteps = [] - resultsdakotas = [] - - for class_instance_name in md.results.__dict__.keys(): + items = obj.__dict__.items() + except AttributeError: + return # not an object with attributes + + for fname, val in items: + if fname in ('errlog', 'outlog'): + continue + _write_value(val, fname, grp, root_nc, f'{path}.{fname}', verbose) + + +# ===================================================================== +# _write_value – dispatch a single value to the right writer +# ===================================================================== +def _write_value(val, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + + # ---- None / empty ----------------------------------------------- + if val is None: + _write_scalar_nan(vname, grp, root_nc) if verbose: - print(class_instance_name) - # for each class instance in results, see which class its from and record that info in the NCData to recreate structure later - # check to see if there is a solutionstep class instance - if isinstance(md.results.__dict__[class_instance_name],solutionstep): - quality_control.append(1) - solutionsteps.append(class_instance_name) - - # check to see if there is a solution class instance - if isinstance(md.results.__dict__[class_instance_name],solution): - quality_control.append(1) - solutions.append(class_instance_name) - - # check to see if there is a resultsdakota class instance - if isinstance(md.results.__dict__[class_instance_name],resultsdakota): - quality_control.append(1) - resultsdakotas.append(class_instance_name) - - if solutionsteps != []: - serialize_string(variable_name=str('solutionstep'), address_of_child=solutionsteps, group=NCData.groups['results'], list=True, NCData=NCData, verbose=verbose) - - if solutions != []: - serialize_string(variable_name=str('solution'), address_of_child=solutions, group=NCData.groups['results'], list=True, NCData=NCData, verbose=verbose) - - if resultsdakotas != []: - serialize_string(variable_name=str('resultsdakota'), address_of_child=resultsdakotas, group=NCData.groups['results'], list=True, NCData=NCData, verbose=verbose) - - - if len(quality_control) != len(md.results.__dict__.keys()): - print('Error: The class instance within your md.results class is not currently supported by this application') - print(type(md.results.__dict__[class_instance_name])) - else: + print(f' [None] {path}') + return + + # ---- numpy ndarray ---------------------------------------------- + if isinstance(val, np.ndarray): + _write_numpy(val, vname, grp, root_nc, path, verbose) + return + + # ---- bool (must come before int because bool is a subclass of int) + if isinstance(val, (bool, np.bool_)): + _write_bool(np.array([int(val)], dtype=np.int16), vname, grp, root_nc) if verbose: - print('The results class was successfully stored on disk') - else: pass + print(f' [bool] {path}') + return + # ---- int / float scalars ---------------------------------------- + if isinstance(val, (int, np.integer)): + _write_scalar(int(val), vname, 'i8', grp, root_nc) + if verbose: + print(f' [int] {path}') + return -def create_NetCDF(filename: str, verbose = False): - # If file already exists delete / rename it - if os.path.exists(filename): - print('File {} allready exist'.format(filename)) - - # If so, inqure for a new name or to do delete the existing file - newname = input('Give a new name or "delete" to replace: ') + if isinstance(val, (float, np.floating)): + _write_scalar(float(val), vname, 'f8', grp, root_nc) + if verbose: + print(f' [float] {path}') + return - if newname == 'delete': - os.remove(filename) - else: - print(('New file name is {}'.format(newname))) - filename = newname - else: - # Otherwise create the file and define it globally so other functions can call it - NCData = Dataset(filename, 'w', format='NETCDF4') - NCData.history = 'Created ' + time.ctime(time.time()) - NCData.createDimension('Unlim', None) # unlimited dimension - NCData.createDimension('float', 1) # single integer dimension - NCData.createDimension('int', 1) # single float dimension - + # ---- str -------------------------------------------------------- + if isinstance(val, str): + _write_string(val, vname, grp) + if verbose: + print(f' [string] {path}') + return + + # ---- list ------------------------------------------------------- + if isinstance(val, list): + _write_list(val, vname, grp, root_nc, path, verbose) + return + + # ---- dict / OrderedDict ----------------------------------------- + if isinstance(val, dict): + _write_dict(val, vname, grp, root_nc, path, verbose) + return + + # ---- ISSM sub-class (has __dict__) ------------------------------ + if hasattr(val, '__dict__'): + sub_grp = grp.createGroup(vname) + sub_grp.classtype = type(val).__name__ + _write_obj(val, sub_grp, root_nc, path, verbose) + return + + # ---- fallback --------------------------------------------------- if verbose: - print('Successfully created ' + filename) + print(f' [SKIP] {path} – unsupported type {type(val).__name__}') - return NCData +# ===================================================================== +# _write_list – dispatch a Python list +# ===================================================================== +def _write_list(val: list, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + if len(val) == 0: + # Empty list: store as a scalar NaN sentinel + _write_scalar_nan(vname, grp, root_nc) + return -def walk_through_model(md, empty_model, NCData, verbose= False): - # Iterate over first layer of md attributes and assume this first layer is only classes - for group in md.__dict__.keys(): - address = md.__dict__[group] - empty_address = empty_model.__dict__[group] - # we need to record the layers of the md so we can save them to the NCData file - layers = [group] + first = val[0] - # Recursively walk through subclasses - walk_through_subclasses(address, empty_address, layers, NCData, empty_model, verbose) + # list of ISSM objects → cell_of_objects group + if hasattr(first, '__dict__') and not isinstance(first, dict): + _write_cell_of_objects(val, vname, grp, root_nc, path, verbose) + return + # list of str + if isinstance(first, str): + _write_cell_strings(val, vname, grp) + if verbose: + print(f' [cellstr] {path}') + return -def walk_through_subclasses(address, empty_address, layers: list, NCData, empty_model, verbose = False): - # See if we have an object with keys or a not - try: - address.__dict__.keys() - is_object = True - except: is_object = False # this is not an object with keys - - if is_object: - # enter the subclass, see if it has nested classes and/or attributes - # then compare attributes between mds and write to NCData if they differ - # if subclass found, walk through it and repeat - for child in address.__dict__.keys(): - # record the current location - current_layer = layers.copy() - current_layer.append(child) - - # navigate to child in each md - address_of_child = address.__dict__[child] - - # if the current object is a results. object and has nonzero steps attr it needs special treatment - if isinstance(address_of_child, solution) and len(address_of_child.steps) != 0: - create_group(address_of_child, current_layer, is_struct = True, is_special_list = False, NCData=NCData, verbose = verbose) - - # if the current object is a list of objects (currently only filters for lists/arrays of classes) - elif isinstance(address_of_child, list) and len(address_of_child) > 0 and hasattr(address_of_child[0], '__dict__'): - create_group(address_of_child, current_layer, is_struct = False, is_special_list = True, NCData=NCData, verbose = verbose) - - # if the variable is an array, assume it has relevant data (this is because the next line cannot evaluate "==" with an array) - elif isinstance(address_of_child, np.ndarray): - create_group(address_of_child, current_layer, is_struct = False, is_special_list = False, NCData=NCData, verbose = verbose) - - # see if the child exists in the empty md. If not, record it in the NCData - else: - try: - address_of_child_in_empty_class = empty_address.__dict__[child] - # if that line worked, we can see how the mds' attributes at this layer compare: - - # if the attributes are identical we don't need to save anything - if address_of_child == address_of_child_in_empty_class: - walk_through_subclasses(address_of_child, address_of_child_in_empty_class, current_layer, NCData, empty_model, verbose) - - # If it has been modified, record it in the NCData file - else: - create_group(address_of_child, current_layer, is_struct = False, is_special_list = False, NCData=NCData, verbose = verbose) - walk_through_subclasses(address_of_child, address_of_child_in_empty_class, current_layer, NCData, empty_model, verbose) - - except KeyError: # record in NCData and continue to walk thru md - walk_through_subclasses(address_of_child, empty_address, current_layer, NCData, empty_model, verbose) - create_group(address_of_child, current_layer, is_struct = False, is_special_list = False, NCData=NCData, verbose = verbose) - else: pass - - -def create_group(address_of_child, layers, is_struct = False, is_special_list = False, NCData=None, verbose = False): - - # Handle the first layer of the group(s) - group_name = layers[0] - - # try to make a group unless the group is already made + # list of solution-step objects (results.solution.steps) + # This is handled by _write_solution_steps when the parent is a solution(); + # if we land here generically, treat as cell_of_objects + if hasattr(first, '__dict__'): + _write_cell_of_objects(val, vname, grp, root_nc, path, verbose) + return + + # Plain numeric list → convert to numpy array try: - group = NCData.createGroup(str(group_name)) - except: - group = NCData.groups[str(group_name)] - - # need to check if inversion or m1qn3inversion class - if group_name == 'inversion': - check_inversion_class(address_of_child, NCData, verbose) - else: pass - - # if the data is nested in md, create nested groups to match class structure - if len(layers) > 2: - for name in layers[1:-1]: - try: - group = group.createGroup(str(name)) - except: - group = NCData.groups[str(name)] - else: pass - - # Lastly, handle the variable(s) - if is_struct: - parent_struct_name = layers[-1] - serialize_nested_results_struct(parent_struct_name, address_of_child, group, NCData, verbose) - - elif is_special_list: - list_name = layers[-1] - serialize_array_of_objects(list_name, address_of_child, group, NCData, verbose) - - else: - variable_name = layers[-1] - serialize_var(variable_name, address_of_child, group, NCData, verbose) - - -def singleton(func): - """ - A decorator to ensure a function is only executed once. - """ - def wrapper(*args, **kwargs): - if not wrapper.has_run: - wrapper.result = func(*args, **kwargs) - wrapper.has_run = True - return wrapper.result - wrapper.has_run = False - wrapper.result = None - return wrapper - - -@singleton -def check_inversion_class(address_of_child, NCData, verbose = False): - # need to make sure that we have the right inversion class: inversion, m1qn3inversion, taoinversion - if isinstance(address_of_child, m1qn3inversion): - serialize_string(variable_name=str('inversion_class_name'), address_of_child=str('m1qn3inversion'), group=NCData.groups['inversion'], NCData=NCData, verbose = verbose) + arr = np.array(val) + _write_numpy(arr, vname, grp, root_nc, path, verbose) + except Exception: if verbose: - print('Successfully saved inversion class instance ' + 'm1qn3inversion') - elif isinstance(address_of_child, taoinversion): - serialize_string(variable_name=str('inversion_class_name'), address_of_child=str('taoinversion'), group=NCData.groups['inversion'], NCData=NCData, verbose = verbose) - if verbose: - print('Successfully saved inversion class instance ' + 'taoinversion') - else: - serialize_string(variable_name=str('inversion_class_name'), address_of_child=str('inversion'), group=NCData.groups['inversion'], NCData=NCData, verbose = verbose) - if verbose: - print('Successfully saved inversion class instance ' + 'inversion') + print(f' [SKIP] {path} – list of mixed/unsupported types') -def serialize_nested_results_struct(parent_struct_name, address_of_struct, group, NCData, verbose = False): - ''' - This function takes a results.solution class instance and saves the solutionstep instances from .steps to the NCData. +# ===================================================================== +# Solution / solutionstep special handling +# ===================================================================== +def _write_obj_with_results(obj, grp, root_nc, path: str, verbose: bool) -> None: + """Like _write_obj but handles solution().steps lists as step_k subgroups.""" + try: + items = obj.__dict__.items() + except AttributeError: + return + + for fname, val in items: + if fname in ('errlog', 'outlog'): + continue + + # Detect a list of solutionstep-like objects + if (isinstance(val, list) and len(val) > 0 + and hasattr(val[0], '__dict__') + and not isinstance(val[0], dict)): + _write_solution_steps(val, fname, grp, root_nc, path, verbose) + else: + _write_value(val, fname, grp, root_nc, f'{path}.{fname}', verbose) - To do this, we get the number of dimensions (substructs) of the parent struct (list). - Next, we iterate through each substruct and record the data. - For each substruct, we create a subgroup of the main struct. - For each variable, we create dimensions that are assigned to each subgroup uniquely. - ''' - if verbose: - print("Beginning transfer of nested MATLAB struct to the NCData") - - # make a new subgroup to contain all the others: - group = group.createGroup(str(parent_struct_name)) - - # make sure other systems can flag the nested struct type - serialize_string('this_is_a_nested', 'struct', group, list=False, NCData=NCData, verbose = verbose) - - # other systems know the name of the parent struct because it's covered by the results/qmu functions above - no_of_dims = len(address_of_struct) - for substruct in range(0, no_of_dims): - # we start by making subgroups with nice names like "1x4" - name_of_subgroup = '1x' + str(substruct) - subgroup = group.createGroup(str(name_of_subgroup)) - - # do some housekeeping to keep track of the current layer - current_substruct = address_of_struct[substruct] - substruct_fields = current_substruct.__dict__.keys() - - # now we need to iterate over each variable of the nested struct and save it to this new subgroup - for variable in substruct_fields: - address_of_child = current_substruct.__dict__[variable] - serialize_var(variable, address_of_child, subgroup, NCData, verbose = verbose) - - if verbose: - print(f'Successfully transferred struct {parent_struct_name} to the NCData\n') +def _write_solution_steps(steps: list, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + """Write a list of solutionstep objects as step_1, step_2, … sub-groups.""" + sub_grp = grp.createGroup(vname) + sub_grp.classtype = 'struct' # MATLAB struct array convention + sub_grp.nsteps = len(steps) + for k, step in enumerate(steps, start=1): + sg = sub_grp.createGroup(f'step_{k}') + for sfield, sval in step.__dict__.items(): + if sfield in ('errlog', 'outlog'): + continue + _write_value(sval, sfield, sg, root_nc, f'{path}.step_{k}.{sfield}', verbose) + if verbose: + print(f' [struct] {path}.{vname} ({len(steps)} steps)') -def serialize_array_of_objects(list_name, address_of_child, group, NCData, verbose): - if verbose: - print(f"Serializing array of objects.") - - # Get the dimensions of the cell array - if len(np.shape(address_of_child)) > 1: - rows, cols = np.shape(address_of_child) - else: rows, cols = 1, np.shape(address_of_child)[0] - # Make subgroup to represent the array - name_of_subgroup = f"{str(rows)}x{str(cols)}_cell_array_of_objects" - subgroup = group.createGroup(name_of_subgroup) +# ===================================================================== +# Leaf writers +# ===================================================================== - # Save the name of the cell array - serialize_string('name_of_cell_array', list_name, subgroup, NCData, verbose) +def _write_numpy(arr: np.ndarray, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + if arr.size == 0: + _write_scalar_nan(vname, grp, root_nc) + return - # Save the dimensions of the cell array - rowsID = subgroup.createVariable('rows', int, ('int',)) - colsID = subgroup.createVariable('cols', int, ('int',)) - rowsID[:] = rows - colsID[:] = cols + # Booleans → int16 + if arr.dtype == bool: + arr = arr.astype(np.int16) + is_bool = True + else: + is_bool = False + # Squeeze trailing size-1 dimensions; keep 0-d as scalar + arr = np.squeeze(arr) + if arr.ndim == 0: + arr = arr.reshape(1) - # If this is a multidimensional cell array, iterate over rows here and cols in serialize_objects - if rows > 1: - for row in range(rows): - # Make a subgroup for each row - name_of_subgroup = f"Row_{row+1}_of_{rows}" - subgroup = group.createGroup(name_of_subgroup) - serialize_objects(address_of_child, subgroup, NCData, cols, verbose) + # Determine NetCDF dtype + if np.issubdtype(arr.dtype, np.integer): + nc_dtype = 'i8' else: - serialize_objects(address_of_child, subgroup, NCData, cols, verbose) - + nc_dtype = 'f8' + arr = arr.astype(np.float64) + + dims = _make_dims(arr.shape, grp, root_nc) + + var = grp.createVariable(vname, nc_dtype, dims, zlib=(arr.size > 1)) + if is_bool: + var.type_is = 'bool' + var[:] = arr + if verbose: - print(f"Successfully serialized array of objects: {list_name}") + print(f' [numeric] {path} shape={arr.shape} dtype={arr.dtype}') -def serialize_objects(address_of_child, group, NCData, cols, verbose): - for col in range(cols): - # Make subgroup to contain each col of array - name_of_subgroup = f'Col_{col+1}_of_{cols}' - subgroup = group.createGroup(name_of_subgroup) +def _write_scalar(val, vname: str, nc_dtype: str, grp, root_nc) -> None: + var = grp.createVariable(vname, nc_dtype, ('scalar',)) + var[:] = val - # index the current item - variable = address_of_child[col] - # Get the kind of object we're working with: - # see if it's a solution instance - if isinstance(variable, solution) and len(variable.steps) != 0: - pass - # this needs more work... - - # see if it's a general class -- assume ISSM classes all have __dict__ - elif hasattr(variable, '__dict__'): - # Handle class instances - serialize_class_instance(variable, subgroup, NCData, verbose) - else: - print('ERROR: Cell arrays of mixed types are not yet supported in read_NCData!') - print('Deserialization will not be able to complete!') - # Handle regular data structures that are already supported - serialize_var(variable_name, variable, subgroup, NCData, verbose) - - -def serialize_class_instance(instance, group, NCData, verbose): - # get parent class name: - name = instance.__class__.__name__ - - # save the name of the class - serialize_string(variable_name='class_is_a', address_of_child=name, group=group, NCData=NCData, verbose = verbose) - - # make subgroup to contain attributes - name_of_subgroup = 'Properties_of_' + name - subgroup = group.createGroup(name_of_subgroup) - - # get attributes - keys = instance.__dict__.keys() - - for name in keys: - serialize_var(name, instance.__dict__[name], subgroup, NCData, verbose) - - - - -def serialize_var(variable_name, address_of_child, group, NCData, verbose = False): - # There are lots of different variable types that we need to handle from the md class - - # This first conditional statement will catch numpy arrays of any dimension and save them - if isinstance(address_of_child, np.ndarray): - serialize_numpy_array(variable_name, address_of_child, group, NCData, verbose=verbose) - - # check if it's an int - elif isinstance(address_of_child, int) or isinstance(address_of_child, np.integer): - variable = group.createVariable(variable_name, int, ('int',)) - variable[:] = address_of_child - - # or a float - elif isinstance(address_of_child, float) or isinstance(address_of_child, np.floating): - variable = group.createVariable(variable_name, float, ('float',)) - variable[:] = address_of_child - - # or a string - elif isinstance(address_of_child, str): - serialize_string(variable_name, address_of_child, group, NCData, verbose=verbose) - - #or a bool - elif isinstance(address_of_child, bool) or isinstance(address_of_child, np.bool_): - # NetCDF can't handle bool types like True/False so we convert all to int 1/0 and add an attribute named units with value 'bool' - variable = group.createVariable(variable_name, int, ('int',)) - variable[:] = int(address_of_child) - variable.units = "bool" - - # or an empty list - elif isinstance(address_of_child, list) and len(address_of_child)==0: - variable = group.createVariable(variable_name, int, ('int',)) - - # or a list of strings -- this needs work as it can only handle a list of 1 string - elif isinstance(address_of_child,list) and isinstance(address_of_child[0],str): - for string in address_of_child: - serialize_string(variable_name, string, group, list=True, NCData=NCData, verbose=verbose) - - # or a regular list - elif isinstance(address_of_child, list): - variable = group.createVariable(variable_name, type(address_of_child[0]), ('Unlim',)) - variable[:] = address_of_child - - # anything else... (will likely need to add more cases; ie helpers.OrderedStruct) - else: +def _write_scalar_nan(vname: str, grp, root_nc) -> None: + var = grp.createVariable(vname, 'f8', ('scalar',)) + var[:] = np.nan + var.type_is = 'empty' + + +def _write_string(s: str, vname: str, grp) -> None: + if len(s) == 0: + dim_name = 'char0' try: - variable = group.createVariable(variable_name, type(address_of_child), ('Unlim',)) - variable[:] = address_of_child - print(f'Unrecognized variable was saved {variable_name}') - except TypeError: pass # this would mean that we have an object, so we just let this continue to feed thru the recursive function above - except Exception as e: - print(f'There was error with {variable_name} in {group}') - print("The error message is:") - print(e) - print('Datatype given: ' + str(type(address_of_child))) - - if verbose: - print(f'Successfully transferred data from {variable_name} to the NCData') - - -def serialize_string(variable_name, address_of_child, group, list=False, NCData=None, verbose = False): - # NCData and strings dont get along.. we have to do it 'custom': - # if we hand it an address we need to do it this way: - if list: - """ - Convert a list of strings to a numpy.char_array with utf-8 encoded elements - and size rows x cols with each row the same # of cols and save to NCData - as char array. - """ + grp.createDimension(dim_name, 0) + except Exception: + pass + v = grp.createVariable(vname, 'S1', (dim_name,)) + v.type_is = 'string' + return + n = len(s) + dim_name = f'char{n}' + try: + grp.createDimension(dim_name, n) + except Exception: + pass + v = grp.createVariable(vname, 'S1', (dim_name,)) + from netCDF4 import stringtochar + v[:] = stringtochar(np.array([s], dtype=f'S{n}')) + v.type_is = 'string' + + +def _write_cell_strings(strings: list, vname: str, grp) -> None: + if not strings: + dim_name = 'char0' try: - strings = address_of_child - # get dims of array to save - rows = len(strings) - cols = len(max(strings, key = len)) - - # Define dimensions for the strings - rows_name = 'rows' + str(rows) - cols_name = 'cols' + str(cols) - try: - group.createDimension(rows_name, rows) - except: pass - - try: - group.createDimension(cols_name, cols) - except: pass - - # Create a variable to store the strings - string_var = group.createVariable(str(variable_name), 'S1', (rows_name, cols_name)) - - # break the list into a list of lists of words with the same length as the longest word: - # make words same sizes by adding spaces - modded_strings = [word + ' ' * (len(max(strings, key=len)) - len(word)) for word in strings] - # encoded words into list of encoded lists - new_list = [[s.encode('utf-8') for s in word] for word in modded_strings] - - # make numpy char array with dims rows x cols - arr = np.chararray((rows, cols)) - - # fill array with list of encoded lists - for i in range(len(new_list)): - arr[i] = new_list[i] - - # save array to NCData file - string_var[:] = arr - - if verbose: - print(f'Saved {len(modded_strings)} strings to {variable_name}') - - except Exception as e: - print(f'Error: {e}') - + grp.createDimension(dim_name, 0) + except Exception: + pass + v = grp.createVariable(vname, 'S1', (dim_name,)) + v.type_is = 'cell_of_strings' + return + + nrows = len(strings) + max_len = max(len(s) for s in strings) + if max_len == 0: + max_len = 1 + + rows_name = f'strrows_{vname}' + cols_name = f'strcols_{vname}' + try: + grp.createDimension(rows_name, nrows) + except Exception: + pass + try: + grp.createDimension(cols_name, max_len) + except Exception: + pass + + v = grp.createVariable(vname, 'S1', (rows_name, cols_name)) + v.type_is = 'cell_of_strings' + + arr = np.chararray((nrows, max_len)) + for i, s in enumerate(strings): + padded = s.ljust(max_len) + arr[i] = list(padded.encode('utf-8')) + v[:] = arr + + +def _write_bool(arr: np.ndarray, vname: str, grp, root_nc) -> None: + if arr.size == 1: + v = grp.createVariable(vname, 'i2', ('scalar',)) + v[:] = arr[0] else: - the_string_to_save = address_of_child - length_of_the_string = len(the_string_to_save) - numpy_datatype = 'S' + str(length_of_the_string) - str_out = netCDF4.stringtochar(np.array([the_string_to_save], dtype=numpy_datatype)) - - # we'll need to make a new dimension for the string if it doesn't already exist - name_of_dimension = 'char' + str(length_of_the_string) - try: - group.createDimension(name_of_dimension, length_of_the_string) - except: pass - # this is another band-aid to the results sub classes... + n = arr.size + dim_name = f'booldim{n}' try: - # now we can make a variable in this dimension: - string = group.createVariable(variable_name, 'S1', (name_of_dimension)) - #finally we can write the variable: - string[:] = str_out - #except RuntimeError: pass - except Exception as e: - print(f'There was an error saving a string from {variable_name}') - print(e) - - -def serialize_numpy_array(variable_name, address_of_child, group, NCData, verbose = False): - # to make a nested array in NCData, we have to get the dimensions of the array, - # create corresponding dimensions in the NCData file, then we can make a variable - # in the NCData with dimensions identical to those in the original array - - # start by getting the data type at the lowest level in the array: - typeis = address_of_child.dtype - - # catch boolean arrays here - if typeis == bool: - # sometimes an array has just 1 element in it, we account for those cases here: - if len(address_of_child) == 1: - variable = group.createVariable(variable_name, int, ('int',)) - variable[:] = int(address_of_child) - variable.units = "bool" - else: - # make the dimensions - dimensions = [] - for dimension in np.shape(address_of_child): - dimensions.append(str('dim' + str(dimension))) - # if the dimension already exists we can't have a duplicate - try: - group.createDimension(str('dim' + str(dimension)), dimension) - except: pass # this would mean that the dimension already exists - - # create the variable: - variable = group.createVariable(variable_name, int, tuple(dimensions)) - # write the variable: - variable[:] = address_of_child.astype(int) - variable.units = "bool" - - # handle all other datatypes here - else: - # sometimes an array has just 1 element in it, we account for those cases here: - if len(address_of_child) == 1: - if typeis is np.dtype('float64'): - variable = group.createVariable(variable_name, typeis, ('float',)) - variable[:] = address_of_child[0] - elif typeis is np.dtype('int64'): - variable = group.createVariable(variable_name, typeis, ('int',)) - variable[:] = address_of_child[0] - else: - print(f'Encountered single datatype from {variable_name} that was not float64 or int64, saving under unlimited dimension, may cause errors.') - variable = group.createVariable(variable_name, typeis, ('Unlim',)) - variable[:] = address_of_child[0] - - # This catches all arrays/lists: + grp.createDimension(dim_name, n) + except Exception: + pass + v = grp.createVariable(vname, 'i2', (dim_name,)) + v[:] = arr.flatten() + v.type_is = 'bool' + + +def _write_cell_of_objects(val: list, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + sub_grp = grp.createGroup(vname) + sub_grp.classtype = 'cell_of_objects' + sub_grp.nrows = 1 + sub_grp.ncols = len(val) + + for c, elem in enumerate(val, start=1): + eg = sub_grp.createGroup(f'item_1_{c}') + eg.classtype = type(elem).__name__ + if hasattr(elem, '__dict__'): + _write_obj(elem, eg, root_nc, f'{path}[{c}]', verbose) + + if verbose: + print(f' [cell] {path} ({len(val)} objects)') + + +def _write_dict(val: dict, vname: str, grp, root_nc, path: str, verbose: bool) -> None: + """Store a dict as a sub-group with one string variable per key.""" + sub_grp = grp.createGroup(vname) + sub_grp.classtype = 'dict' + for k, v in val.items(): + _write_value(v, str(k), sub_grp, root_nc, f'{path}.{k}', verbose) + + +# ===================================================================== +# Dimension helper +# ===================================================================== + +def _make_dims(shape: tuple, grp, root_nc) -> tuple: + """Return a tuple of dimension-name strings, creating dims as needed.""" + dim_names = [] + for n in shape: + if n == 1: + dim_names.append('scalar') else: - # make the dimensions - dimensions = [] - for dimension in np.shape(address_of_child): - dimensions.append(str('dim' + str(dimension))) - # if the dimension already exists we can't have a duplicate + dname = f'dim{n}' + # Try to create in local group; if it exists anywhere, reuse it + for target in (grp, root_nc): try: - group.createDimension(str('dim' + str(dimension)), dimension) - except: pass # this would mean that the dimension already exists - - # create the variable: - variable = group.createVariable(variable_name, typeis, tuple(dimensions)) - - # write the variable: - variable[:] = address_of_child - - \ No newline at end of file + target.createDimension(dname, n) + break + except Exception: + pass # already exists – that's fine + dim_names.append(dname) + return tuple(dim_names)