nrncore_io.cpp

Go to the documentation of this file.

#include "nrncore_io.h"
#include "nrncore_write/data/cell_group.h"
#include "nrncore_write/callbacks/nrncore_callbacks.h"
 
#include <cstdlib>
#include "nrnmpi.h"
#include "section.h"
#include "hocdec.h"
#include "ocfile.h"     // for idDirExist and makePath
#include "nrnran123.h"  // globalindex written to globals.dat
#include "cvodeobj.h"
#include "netcvode.h"  // for nrnbbcore_vecplay_write
#include "vrecitem.h"  // for nrnbbcore_vecplay_write
#include <fstream>
#include <sstream>
#include "nrnsection_mapping.h"
 
extern short* nrn_is_artificial_;
extern int* bbcore_dparam_size;
extern bbcore_write_t* nrn_bbcore_write_;
extern NetCvode* net_cvode_instance;
extern void (*nrnthread_v_transfer_)(NrnThread*);
 
int chkpnt;
const char* bbcore_write_version = "1.5";  // Generalize POINTER to allow pointing to any RANGE
                                           // variable
 
/// create directory with given path
void create_dir_path(const std::string& path) {
    // only one rank needs to create directory
    if (nrnmpi_myid == 0) {
        if (!isDirExist(path)) {
            if (!makePath(path)) {
                hoc_execerror(path.c_str(), "directory did not exist and makePath for it failed");
            }
        }
    }
    // rest of the ranks should wait before continue simulation
#ifdef NRNMPI
    nrnmpi_barrier();
#endif
}
 
std::string get_write_path() {
    std::string path(".");  // default path
    if (ifarg(1)) {
        path = hoc_gargstr(1);
    }
    return path;
}
 
std::string get_filename(const std::string& path, std::string file_name) {
    std::string fname(path + '/' + file_name);
    nrn_assert(fname.size() < 1024);
    return fname;
}
 
 
void write_memb_mech_types(const char* fname) {
    if (nrnmpi_myid > 0) {
        return;
    }  // only rank 0 writes this file
    std::ofstream fs(fname);
    if (!fs.good()) {
        hoc_execerror("nrncore_write write_mem_mech_types could not open for writing: %s\n", fname);
    }
    write_memb_mech_types_direct(fs);
}
 
 
// format is name value
// with last line of 0 0
// In case of an array, the line is name[num] with num lines following with
// one value per line.  Values are %.20g format.
void write_globals(const char* fname) {
    if (nrnmpi_myid > 0) {
        return;
    }  // only rank 0 writes this file
 
    FILE* f = fopen(fname, "w");
    if (!f) {
        hoc_execerror("nrncore_write write_globals could not open for writing: %s\n", fname);
    }
 
    fprintf(f, "%s\n", bbcore_write_version);
    const char* name;
    int size;            // 0 means scalar, is 0 will still allocated one element for val.
    double* val = NULL;  // Allocated by new in get_global_item, must be delete [] here.
    // Note that it is possible for get_global_dbl_item to return NULL but
    // name, size, and val must still be handled if val != NULL
    for (void* sp = NULL;;) {
        sp = get_global_dbl_item(sp, name, size, val);
        if (val) {
            if (size) {
                fprintf(f, "%s[%d]\n", name, size);
                for (int i = 0; i < size; ++i) {
                    fprintf(f, "%.20g\n", val[i]);
                }
            } else {
                fprintf(f, "%s %.20g\n", name, val[0]);
            }
            delete[] val;
            val = NULL;
        }
        if (!sp) {
            break;
        }
    }
    fprintf(f, "0 0\n");
    fprintf(f, "secondorder %d\n", secondorder);
    fprintf(f, "Random123_globalindex %d\n", nrnran123_get_globalindex());
    fprintf(f, "_nrnunit_use_legacy_ %d\n", _nrnunit_use_legacy_);
 
    fclose(f);
}
 
 
void write_nrnthread(const char* path, NrnThread& nt, CellGroup& cg) {
    char fname[1000];
    if (cg.n_output <= 0) {
        return;
    }
    assert(cg.group_id >= 0);
    nrn_assert(snprintf(fname, 1000, "%s/%d_1.dat", path, cg.group_id) < 1000);
    FILE* f = fopen(fname, "wb");
    if (!f) {
        hoc_execerror("nrncore_write write_nrnthread could not open for writing:", fname);
    }
    fprintf(f, "%s\n", bbcore_write_version);
 
    // nrnthread_dat1(int tid, int& n_presyn, int& n_netcon, int*& output_gid, int*& netcon_srcgid);
    fprintf(f, "%d npresyn\n", cg.n_presyn);
    fprintf(f, "%d nnetcon\n", cg.n_netcon);
    writeint(cg.output_gid, cg.n_presyn);
    writeint(cg.netcon_srcgid, cg.n_netcon);
 
    if (cg.output_gid) {
        delete[] cg.output_gid;
        cg.output_gid = NULL;
    }
    if (cg.netcon_srcgid) {
        delete[] cg.netcon_srcgid;
        cg.netcon_srcgid = NULL;
    }
    fclose(f);
 
    nrn_assert(snprintf(fname, 1000, "%s/%d_2.dat", path, cg.group_id) < 1000);
    f = fopen(fname, "w");
    if (!f) {
        hoc_execerror("nrncore_write write_nrnthread could not open for writing:", fname);
    }
 
    fprintf(f, "%s\n", bbcore_write_version);
 
    // sizes and total data count
    int ngid, n_real_gid, nnode, ndiam, nmech, *tml_index, *ml_nodecount, nidata, nvdata, nweight;
    nrnthread_dat2_1(nt.id,
                     ngid,
                     n_real_gid,
                     nnode,
                     ndiam,
                     nmech,
                     tml_index,
                     ml_nodecount,
                     nidata,
                     nvdata,
                     nweight);
 
    fprintf(f, "%d ngid\n", ngid);
    fprintf(f, "%d n_real_gid\n", n_real_gid);
    fprintf(f, "%d nnode\n", nnode);
    fprintf(f, "%d ndiam\n", ndiam);
    fprintf(f, "%d nmech\n", nmech);
 
    for (int i = 0; i < nmech; ++i) {
        fprintf(f, "%d\n", tml_index[i]);
        fprintf(f, "%d\n", ml_nodecount[i]);
    }
    delete[] tml_index;
    delete[] ml_nodecount;
 
    fprintf(f, "%d nidata\n", 0);
    fprintf(f, "%d nvdata\n", nvdata);
    fprintf(f, "%d nweight\n", nweight);
 
    // data
    int* v_parent_index = NULL;
    double *a = NULL, *b = NULL, *area = NULL, *v = NULL, *diamvec = NULL;
    nrnthread_dat2_2(nt.id, v_parent_index, a, b, area, v, diamvec);
    assert(cg.n_real_output == nt.ncell);
    writeint(nt._v_parent_index, nt.end);
    writedbl(nt._actual_a, nt.end);
    writedbl(nt._actual_b, nt.end);
    writedbl(nt._actual_area, nt.end);
    writedbl(nt._actual_v, nt.end);
    if (cg.ndiam) {
        writedbl(diamvec, nt.end);
        delete[] diamvec;
    }
 
    // mechanism data
    int dsz_inst = 0;
    MlWithArt& mla = cg.mlwithart;
    for (size_t i = 0; i < mla.size(); ++i) {
        int type = mla[i].first;
        int *nodeindices = NULL, *pdata = NULL;
        double* data = NULL;
        std::vector<int> pointer2type;
        nrnthread_dat2_mech(nt.id, i, dsz_inst, nodeindices, data, pdata, pointer2type);
        Memb_list* ml = mla[i].second;
        int n = ml->nodecount;
        int sz = nrn_prop_param_size_[type];
        if (nodeindices) {
            writeint(nodeindices, n);
        }
        writedbl(data, n * sz);
        if (nrn_is_artificial_[type]) {
            delete[] data;
        }
        sz = bbcore_dparam_size[type];
        if (pdata) {
            ++dsz_inst;
            writeint(pdata, n * sz);
            delete[] pdata;
            sz = pointer2type.size();
            fprintf(f, "%d npointer\n", int(sz));
            if (sz > 0) {
                writeint(pointer2type.data(), sz);
            }
        }
    }
 
    int *output_vindex, *netcon_pnttype, *netcon_pntindex;
    double *output_threshold, *weights, *delays;
    nrnthread_dat2_3(nt.id,
                     nweight,
                     output_vindex,
                     output_threshold,
                     netcon_pnttype,
                     netcon_pntindex,
                     weights,
                     delays);
    writeint(output_vindex, cg.n_presyn);
    delete[] output_vindex;
    writedbl(output_threshold, cg.n_real_output);
    delete[] output_threshold;
 
    // connections
    int n = cg.n_netcon;
    // printf("n_netcon=%d nweight=%d\n", n, nweight);
    writeint(netcon_pnttype, n);
    delete[] netcon_pnttype;
    writeint(netcon_pntindex, n);
    delete[] netcon_pntindex;
    writedbl(weights, nweight);
    delete[] weights;
    writedbl(delays, n);
    delete[] delays;
 
    // special handling for BBCOREPOINTER
    // how many mechanisms require it
    nrnthread_dat2_corepointer(nt.id, n);
    fprintf(f, "%d bbcorepointer\n", n);
    // for each of those, what is the mech type and data size
    // and what is the data
    for (size_t i = 0; i < mla.size(); ++i) {
        int type = mla[i].first;
        if (nrn_bbcore_write_[type]) {
            int icnt, dcnt, *iArray;
            double* dArray;
            nrnthread_dat2_corepointer_mech(nt.id, type, icnt, dcnt, iArray, dArray);
            fprintf(f, "%d\n", type);
            fprintf(f, "%d\n%d\n", icnt, dcnt);
            if (icnt) {
                writeint(iArray, icnt);
                delete[] iArray;
            }
            if (dcnt) {
                writedbl(dArray, dcnt);
                delete[] dArray;
            }
        }
    }
 
    nrnbbcore_vecplay_write(f, nt);
 
    fclose(f);
}
 
 
void writeint_(int* p, size_t size, FILE* f) {
    fprintf(f, "chkpnt %d\n", chkpnt++);
    size_t n = fwrite(p, sizeof(int), size, f);
    assert(n == size);
}
 
void writedbl_(double* p, size_t size, FILE* f) {
    fprintf(f, "chkpnt %d\n", chkpnt++);
    size_t n = fwrite(p, sizeof(double), size, f);
    assert(n == size);
}
 
#define writeint(p, size) writeint_(p, size, f)
#define writedbl(p, size) writedbl_(p, size, f)
 
void write_contiguous_art_data(double** data, int nitem, int szitem, FILE* f) {
    fprintf(f, "chkpnt %d\n", chkpnt++);
    // the assumption is that an fwrite of nitem groups of szitem doubles can be
    // fread as a single group of nitem*szitem doubles.
    for (int i = 0; i < nitem; ++i) {
        size_t n = fwrite(data[i], sizeof(double), szitem, f);
        assert(n == szitem);
    }
}
 
double* contiguous_art_data(double** data, int nitem, int szitem) {
    double* d1 = new double[nitem * szitem];
    int k = 0;
    for (int i = 0; i < nitem; ++i) {
        for (int j = 0; j < szitem; ++j) {
            d1[k++] = data[i][j];
        }
    }
    return d1;
}
 
 
void nrnbbcore_vecplay_write(FILE* f, NrnThread& nt) {
    // Get the indices in NetCvode.fixed_play_ for this thread
    // error if not a VecPlayContinuous with no discon vector
    std::vector<int> indices;
    nrnthread_dat2_vecplay(nt.id, indices);
    fprintf(f, "%d VecPlay instances\n", int(indices.size()));
    for (auto i: indices) {
        int vptype, mtype, ix, sz;
        double *yvec, *tvec;
        // the 'if' is not necessary as item i is certainly in this thread
        int unused = 0;
        if (nrnthread_dat2_vecplay_inst(
                nt.id, i, vptype, mtype, ix, sz, yvec, tvec, unused, unused, unused)) {
            fprintf(f, "%d\n", vptype);
            fprintf(f, "%d\n", mtype);
            fprintf(f, "%d\n", ix);
            fprintf(f, "%d\n", sz);
            writedbl(yvec, sz);
            writedbl(tvec, sz);
        }
    }
}
 
 
static void fgets_no_newline(char* s, int size, FILE* f) {
    if (fgets(s, size, f) == NULL) {
        fclose(f);
        hoc_execerror("Error reading line in files.dat", strerror(errno));
    }
    int n = strlen(s);
    if (n && s[n - 1] == '\n') {
        s[n - 1] = '\0';
    }
}
 
/** Write all dataset ids to files.dat.
 *
 * Format of the files.dat file is:
 *
 *     version string
 *     -1 (if model uses gap junction)
 *     n (number of datasets) in format %10d
 *     id1
 *     id2
 *     ...
 *     idN
 */
void write_nrnthread_task(const char* path, CellGroup* cgs, bool append) {
    // ids of datasets that will be created
    std::vector<int> iSend;
 
    // ignore empty nrnthread (has -1 id)
    for (int iInt = 0; iInt < nrn_nthread; ++iInt) {
        if (cgs[iInt].group_id >= 0) {
            iSend.push_back(cgs[iInt].group_id);
        }
    }
 
    // receive and displacement buffers for mpi
    std::vector<int> iRecv, iDispl;
 
    if (nrnmpi_myid == 0) {
        iRecv.resize(nrnmpi_numprocs);
        iDispl.resize(nrnmpi_numprocs);
    }
 
    // number of datasets on the current rank
    int num_datasets = iSend.size();
 
#ifdef NRNMPI
    // gather number of datasets from each task
    if (nrnmpi_numprocs > 1) {
        nrnmpi_int_gather(&num_datasets, begin_ptr(iRecv), 1, 0);
    } else {
        iRecv[0] = num_datasets;
    }
#else
    iRecv[0] = num_datasets;
#endif
 
    // total number of datasets across all ranks
    int iSumThread = 0;
 
    // calculate mpi displacements
    if (nrnmpi_myid == 0) {
        for (int iInt = 0; iInt < nrnmpi_numprocs; ++iInt) {
            iDispl[iInt] = iSumThread;
            iSumThread += iRecv[iInt];
        }
    }
 
    // buffer for receiving all dataset ids
    std::vector<int> iRecvVec(iSumThread);
 
#ifdef NRNMPI
    // gather ids into the array with correspondent offsets
    if (nrnmpi_numprocs > 1) {
        nrnmpi_int_gatherv(begin_ptr(iSend),
                           num_datasets,
                           begin_ptr(iRecvVec),
                           begin_ptr(iRecv),
                           begin_ptr(iDispl),
                           0);
    } else {
        for (int iInt = 0; iInt < num_datasets; ++iInt) {
            iRecvVec[iInt] = iSend[iInt];
        }
    }
#else
    for (int iInt = 0; iInt < num_datasets; ++iInt) {
        iRecvVec[iInt] = iSend[iInt];
    }
#endif
 
    /// Writing the file with task, correspondent number of threads and list of correspondent first
    /// gids
    if (nrnmpi_myid == 0) {
        // If append is false, begin a new files.dat (overwrite old if exists).
        // If append is true, append groupids to existing files.dat.
        //   Note: The number of groupids (2nd or 3rd line) has to be
        //   overwritten wih the total number so far. To avoid copying
        //   old to new, we allocate 10 chars for that number.
 
        std::stringstream ss;
        ss << path << "/files.dat";
 
        std::string filename = ss.str();
 
        FILE* fp = NULL;
        if (append == false) {  // start a new file
            fp = fopen(filename.c_str(), "w");
            if (!fp) {
                hoc_execerror("nrncore_write: could not open for writing:", filename.c_str());
            }
        } else {  // modify groupid number and append to existing file
            fp = fopen(filename.c_str(), "r+");
            if (!fp) {
                hoc_execerror("nrncore_write append: could not open for modifying:",
                              filename.c_str());
            }
        }
 
        constexpr int max_line_len = 20;
        char line[max_line_len];  // All lines are actually no larger than %10d.
 
        if (append) {
            // verify same version
            fgets_no_newline(line, max_line_len, fp);
            // unfortunately line has the newline
            size_t n = strlen(bbcore_write_version);
            if ((strlen(line) != n) || strncmp(line, bbcore_write_version, n) != 0) {
                fclose(fp);
                hoc_execerror("nrncore_write append: existing files.dat has inconsisten version:",
                              line);
            }
        } else {
            fprintf(fp, "%s\n", bbcore_write_version);
        }
 
        // notify coreneuron that this model involves gap junctions
        if (nrnthread_v_transfer_) {
            if (append) {
                fgets_no_newline(line, max_line_len, fp);
                if (strcmp(line, "-1") != 0) {
                    fclose(fp);
                    hoc_execerror(
                        "nrncore_write append: existing files.dat does not have a gap junction "
                        "indicator\n",
                        NULL);
                }
            } else {
                fprintf(fp, "-1\n");
            }
        }
 
        // total number of datasets
        if (append) {
            // this is the one that needs the space to get a new value
            long pos = ftell(fp);
            fgets_no_newline(line, max_line_len, fp);
            int oldval = 0;
            if (sscanf(line, "%d", &oldval) != 1) {
                fclose(fp);
                hoc_execerror("nrncore_write append: error reading number of groupids", NULL);
            }
            if (oldval == -1) {
                fclose(fp);
                hoc_execerror(
                    "nrncore_write append: existing files.dat has gap junction indicator where we "
                    "expected a groupgid count.",
                    NULL);
            }
            iSumThread += oldval;
            fseek(fp, pos, SEEK_SET);
        }
        fprintf(fp, "%10d\n", iSumThread);
 
        if (append) {
            // Start writing the groupids starting at the end of the file.
            fseek(fp, 0, SEEK_END);
        }
 
        // write all dataset ids
        for (int i = 0; i < iRecvVec.size(); ++i) {
            fprintf(fp, "%d\n", iRecvVec[i]);
        }
 
        fclose(fp);
    }
}
 
/** @brief dump mapping information to gid_3.dat file */
void nrn_write_mapping_info(const char* path, int gid, NrnMappingInfo& minfo) {
    /** full path of mapping file */
    std::stringstream ss;
    ss << path << "/" << gid << "_3.dat";
 
    std::string fname(ss.str());
    FILE* f = fopen(fname.c_str(), "w");
 
    if (!f) {
        hoc_execerror("nrnbbcore_write could not open for writing:", fname.c_str());
    }
 
    fprintf(f, "%s\n", bbcore_write_version);
 
    /** number of gids in NrnThread */
    fprintf(f, "%zd\n", minfo.size());
 
    /** all cells mapping information in NrnThread */
    for (size_t i = 0; i < minfo.size(); i++) {
        CellMapping* c = minfo.mapping[i];
 
        /** gid, #section, #compartments,  #sectionlists */
        fprintf(f, "%d %d %d %zd\n", c->gid, c->num_sections(), c->num_segments(), c->size());
 
        for (size_t j = 0; j < c->size(); j++) {
            SecMapping* s = c->secmapping[j];
            /** section list name, number of sections, number of segments */
            fprintf(f, "%s %d %zd\n", s->name.c_str(), s->nsec, s->size());
 
            /** section - segment mapping */
            if (s->size()) {
                writeint(&(s->sections.front()), s->size());
                writeint(&(s->segments.front()), s->size());
            }
        }
    }
    fclose(f);
}