https://github.com/cran/sbioPN
Tip revision: cb9bccd9cb9d3b21276ae307fb246f69810e9c43 authored by Roberto Bertolusso on 15 March 2014, 00:00:00 UTC
version 1.1.0
version 1.1.0
Tip revision: cb9bccd
sGillespieDirectCR.c
// -*- mode: C; c-indent-level: 2; c-basic-offset: 2; tab-width: 8 -*-
//
// Copyright (C) 2009-2014 Roberto Bertolusso and Marek Kimmel
//
// This file is part of sbioPN.
//
// sbioPN is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// sbioPN is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with sbioPN. If not, see <http://www.gnu.org/licenses/>.
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include <R_ext/Rdynload.h>
#include <R_ext/Visibility.h>
#include "helper.h"
#define RB_MEMORY
#define RB_TIME
#define RB_SUBTIME
#ifdef RB_TIME
#include <time.h>
#endif
struct transition {
// DO NOT CHANGE ORDER OF ELEMENTS!! (Optimized for total size of structure)
char cPreNZ_len;
char cSNZ_len;
char cType_of_h;
char cInGroup;
int iPositionInGroup;
// NOTE: 4 chars + 1 int use 8 bytes so no space is wasted in 64 bit machines
int *piPreNZ_VoxelPlace;
char *pcPreNZ_value;
int *piSNZ_VoxelPlace;
char *pcSNZ_value;
union {
double dTrCnt;
DL_FUNC *pDL_FUNC;
SEXP SEXP;
} h;
int *piHazardsToMod;
// Next int wastes 4 bytes per element.
int iHazardsToMod_count;
int iLocalVoxelPlace_offset;
};
typedef struct transition transition_el;
// This structure has 7 unused bytes per element.
struct tree_el {
char iGroup;
double dPartialAcumHazard;
struct tree_el *parent, *left, *right;
};
typedef struct tree_el node;
SEXP getListElement(SEXP list, const char *str);
/*
{
SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
int i, len = length(list);
for (i = 0; i < len; i++)
if(strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
elmt = VECTOR_ELT(list, i);
break;
}
return elmt;
}
*/
typedef enum {
HZ_DOUBLE,
HZ_CFUNCTION,
HZ_RFUNCTION
} HZ_type;
SEXP sGillespieDirectCR(SEXP model, SEXP T, SEXP delta, SEXP runs,
SEXP place, SEXP transition, SEXP rho)
{
SEXP sexpTmp;
int *piTmp;
char cTmp;
int iVoxelPlaceIdx;
int iVoxelTransitionIdx;
#ifdef RB_TIME
clock_t c0, c1;
c0 = clock();
clock_t clkStep = 0;
int iElSteps = 0, iTotalStepsOver50, iTotalStepsOver10;
#endif
#ifdef RB_MEMORY
double dUsedMemAcum = 0, dThisMem;
#endif
SEXP VoxelFamily = getListElement(model, "VoxelFamily");
if (VoxelFamily == R_NilValue)
error("VoxelFamily not provided!\n");
// Process groups of Voxels
int iVoxelFamily_len = length(VoxelFamily);
int iVoxelFamily;
int iVoxelTransitions = 0;
int iVoxelPlaces = 0;
// First: Find out how many real Places and Transitions there are.
for (iVoxelFamily = 0; iVoxelFamily < iVoxelFamily_len; iVoxelFamily++) {
SEXP this_VoxelFamily = VECTOR_ELT(VoxelFamily, iVoxelFamily);
if ((sexpTmp = getListElement(this_VoxelFamily, "pre")) == R_NilValue)
error("pre not provided!\n");
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iOrigTransitions = piTmp[0], iOrigPlaces = piTmp[1];
if ((sexpTmp = getListElement(this_VoxelFamily, "Voxels")) == R_NilValue)
error("Voxels not provided!\n");
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iVoxels = INTEGER(sexpTmp)[0];
UNPROTECT_PTR(sexpTmp);
iVoxelPlaces += iOrigPlaces * iVoxels;
if ((sexpTmp = getListElement(this_VoxelFamily, "TransitionInVoxelSubset")) == R_NilValue)
error("TransitionInVoxelSubset not provided!\n");
SEXP TransitionInVoxelSubset;
PROTECT(TransitionInVoxelSubset = coerceVector(sexpTmp, INTSXP));
int *piTransitionInVoxelSubset = INTEGER(TransitionInVoxelSubset);
SEXP VoxelSubset;
if ((VoxelSubset = getListElement(this_VoxelFamily, "VoxelSubset")) == R_NilValue)
error("VoxelSubset not provided!\n");
int iOrigTransition;
for (iOrigTransition = 0; iOrigTransition < iOrigTransitions; iOrigTransition++) {
iVoxelTransitions += length(VECTOR_ELT(VoxelSubset, piTransitionInVoxelSubset[iOrigTransition] - 1));
}
UNPROTECT_PTR(TransitionInVoxelSubset);
SEXP JumpingPlace;
if ((JumpingPlace = getListElement(this_VoxelFamily, "JumpingPlace")) != R_NilValue) {
SEXP JumpingPattern;
if ((JumpingPattern = getListElement(this_VoxelFamily, "JumpingPattern")) == R_NilValue)
error("JumpingPattern not provided!\n");
int iJumpingPlace;
for (iJumpingPlace = 0; iJumpingPlace < length(JumpingPlace); iJumpingPlace++) {
sexpTmp = VECTOR_ELT(VECTOR_ELT(JumpingPlace, iJumpingPlace), 2);
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iJumpingPattern = INTEGER(sexpTmp)[0] - 1;
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(JumpingPattern, iJumpingPattern);
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iJumps = piTmp[0];
iVoxelTransitions += iJumps;
}
}
}
Rprintf("VoxelTransitions: %d\tVoxelPlaces: %d\n", iVoxelTransitions, iVoxelPlaces);
transition_el *pVoxelTransition = (transition_el *) R_alloc(iVoxelTransitions, sizeof(transition_el));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(transition_el) * iVoxelTransitions)/1e6);
Rprintf("pVoxelTransition: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitions, sizeof(transition_el));
#endif
int iVoxelTransition;
int iVoxelTransitionsPreNZ_len = 0;
int iVoxelTransitionsSNZ_len = 0;
int iVoxelPlacesPreNZ_len = 0;
int *piVoxelPlacePreNZ_len = (int *) R_alloc(iVoxelPlaces, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelPlaces)/1e6);
Rprintf("piVoxelPlacePreNZ_len: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlaces, sizeof(int));
#endif
int **ppiVoxelPlacePreNZ_VoxelTransition = (int **) R_alloc(iVoxelPlaces, sizeof(int *));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int *) * iVoxelPlaces)/1e6);
Rprintf("ppiVoxelPlacePreNZ_VoxelTransition: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlaces, sizeof(int *));
#endif
int iVoxelPlacesSNZ_len = 0;
int *piVoxelPlaceSNZ_len = (int *) R_alloc(iVoxelPlaces, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelPlaces)/1e6);
Rprintf("piVoxelPlaceSNZ_len: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlaces, sizeof(int));
#endif
int **ppiVoxelPlaceSNZ_VoxelTransition = (int **) R_alloc(iVoxelPlaces, sizeof(int *));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int *) * iVoxelPlaces)/1e6);
Rprintf("ppiVoxelPlaceSNZ_VoxelTransition: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlaces, sizeof(int *));
#endif
int iVoxelPlace;
for(iVoxelPlace = 0; iVoxelPlace < iVoxelPlaces; iVoxelPlace++) {
piVoxelPlacePreNZ_len[iVoxelPlace] = piVoxelPlaceSNZ_len[iVoxelPlace] = 0;
}
int *piPreNZ_VoxelPlace = NULL, *piSNZ_VoxelPlace = NULL;
char *pcPreNZ_value = NULL, *pcSNZ_value = NULL;
int iVoxelCurrentTransition = 0;
int iVoxelPlace_offset = 0;
// Second:
for (iVoxelFamily = 0; iVoxelFamily < iVoxelFamily_len; iVoxelFamily++) {
//Rprintf("Voxel group: %d\n", iVoxelFamily);
SEXP this_VoxelFamily = VECTOR_ELT(VoxelFamily, iVoxelFamily);
sexpTmp = getListElement(this_VoxelFamily, "pre");
SEXP pre;
PROTECT(pre = coerceVector(sexpTmp, INTSXP));
int *piPre = INTEGER(pre);
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iOrigTransitions = piTmp[0], iOrigPlaces = piTmp[1];
//Rprintf("pre: rows: %d\tcols: %d\n", iOrigTransitions, iOrigPlaces);
sexpTmp = getListElement(this_VoxelFamily, "Voxels");
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iVoxels = INTEGER(sexpTmp)[0];
UNPROTECT_PTR(sexpTmp);
sexpTmp = getListElement(this_VoxelFamily, "post");
SEXP post;
PROTECT(post = coerceVector(sexpTmp, INTSXP));
int *piPost = INTEGER(post);
piPreNZ_VoxelPlace = (int *) Realloc(piPreNZ_VoxelPlace, iOrigPlaces, int);
pcPreNZ_value = (char *) Realloc(pcPreNZ_value, iOrigPlaces, char);
piSNZ_VoxelPlace = (int *) Realloc(piSNZ_VoxelPlace, iOrigPlaces, int);
pcSNZ_value = (char *) Realloc(pcSNZ_value, iOrigPlaces, char);
sexpTmp = getListElement(this_VoxelFamily, "TransitionInVoxelSubset");
SEXP TransitionInVoxelSubset;
PROTECT(TransitionInVoxelSubset = coerceVector(sexpTmp, INTSXP));
int *piTransitionInVoxelSubset = INTEGER(TransitionInVoxelSubset);
SEXP VoxelSubset;
VoxelSubset = getListElement(this_VoxelFamily, "VoxelSubset");
SEXP h;
h = getListElement(this_VoxelFamily, "h");
int iOrigTransition, iOrigPlace;
for (iOrigTransition = 0; iOrigTransition < iOrigTransitions; iOrigTransition++) {
char cPreNZ_len = 0;
char cSNZ_len = 0;
for (iOrigPlace = 0; iOrigPlace < iOrigPlaces; iOrigPlace++) {
if ((cTmp = piPre[iOrigTransition + iOrigTransitions * iOrigPlace])) {
piPreNZ_VoxelPlace[(int) cPreNZ_len] = iOrigPlace;
pcPreNZ_value[(int) cPreNZ_len++] = cTmp;
}
if ((cTmp = piPost[iOrigTransition + iOrigTransitions * iOrigPlace] -
piPre[iOrigTransition + iOrigTransitions * iOrigPlace])) {
piSNZ_VoxelPlace[(int) cSNZ_len] = iOrigPlace;
pcSNZ_value[(int) cSNZ_len++] = cTmp;
}
}
// For the current transition, find the subset of Voxels in which it applies.
sexpTmp = VECTOR_ELT(VoxelSubset, piTransitionInVoxelSubset[iOrigTransition] - 1);
int iVoxelSubset_len = length(sexpTmp);
SEXP VoxelSubset;
PROTECT(VoxelSubset = coerceVector(sexpTmp, INTSXP));
int *piVoxelSubset = INTEGER(VoxelSubset); // Vector with subset of Voxels
int iVoxelSubset;
// Create a corresponding transition on each of the Voxels of this subset.
for (iVoxelSubset = 0; iVoxelSubset < iVoxelSubset_len; iVoxelSubset++) {
int iVoxel = piVoxelSubset[iVoxelSubset] - 1;
int iLocalVoxelPlace_offset = iVoxelPlace_offset + iOrigPlaces * iVoxel;
//Rprintf("Voxel: %d\tTransition: %d\tVoxelTransition: %d\n", iVoxel, iOrigTransition, iVoxelCurrentTransition);
iVoxelTransitionsPreNZ_len += cPreNZ_len;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < cPreNZ_len; iVoxelPlaceIdx++) {
iVoxelPlace = iLocalVoxelPlace_offset + piPreNZ_VoxelPlace[iVoxelPlaceIdx];
piVoxelPlacePreNZ_len[iVoxelPlace]++;
iVoxelPlacesPreNZ_len++;
}
iVoxelTransitionsSNZ_len += cSNZ_len;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < cSNZ_len; iVoxelPlaceIdx++) {
iVoxelPlace = iLocalVoxelPlace_offset + piSNZ_VoxelPlace[iVoxelPlaceIdx];
piVoxelPlaceSNZ_len[iVoxelPlace]++;
iVoxelPlacesSNZ_len++;
}
iVoxelCurrentTransition++;
}
UNPROTECT_PTR(VoxelSubset);
}
UNPROTECT_PTR(TransitionInVoxelSubset);
UNPROTECT_PTR(post);
UNPROTECT_PTR(pre);
// Manage jumps.
SEXP JumpingPlace;
if ((JumpingPlace = getListElement(this_VoxelFamily, "JumpingPlace")) != R_NilValue) {
SEXP JumpingPattern;
JumpingPattern = getListElement(this_VoxelFamily, "JumpingPattern");
int iJumpingPlace;
for (iJumpingPlace = 0; iJumpingPlace < length(JumpingPlace); iJumpingPlace++) {
SEXP this_JumpingPlace = VECTOR_ELT(JumpingPlace, iJumpingPlace);
sexpTmp = VECTOR_ELT(this_JumpingPlace, 0);
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iOrigPlace = INTEGER(sexpTmp)[0] - 1;
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(this_JumpingPlace, 2);
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iJumpingPattern = INTEGER(sexpTmp)[0] - 1;
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(JumpingPattern, iJumpingPattern);
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iJumps = piTmp[0];
SEXP jump;
PROTECT(jump = coerceVector(sexpTmp, INTSXP));
int *piJump = INTEGER(jump);
int iJump;
for (iJump = 0; iJump < iJumps; iJump++) {
int iFromVoxelPlace = iVoxelPlace_offset + iOrigPlaces * (piJump[iJump] - 1) + iOrigPlace;
int iToVoxelPlace = iVoxelPlace_offset + iOrigPlaces * (piJump[iJump + iJumps] - 1) + iOrigPlace;
char cToVoxelPlacePost_value = piJump[iJump + iJumps*2];
iVoxelTransitionsPreNZ_len++;
piVoxelPlacePreNZ_len[iFromVoxelPlace]++;
iVoxelPlacesPreNZ_len++;
iVoxelTransitionsSNZ_len++;
piVoxelPlaceSNZ_len[iFromVoxelPlace]++;
iVoxelPlacesSNZ_len++;
if (cToVoxelPlacePost_value) {
iVoxelTransitionsSNZ_len++;
piVoxelPlaceSNZ_len[iToVoxelPlace]++;
iVoxelPlacesSNZ_len++;
}
iVoxelCurrentTransition++;
//Rprintf("%d\t%d\n", iFromVoxelPlace, iToVoxelPlace);
}
UNPROTECT_PTR(jump);
}
}
iVoxelPlace_offset += iOrigPlaces * iVoxels;
}
int *piThisPreNZ = (int *) R_alloc(iVoxelPlacesPreNZ_len, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelPlacesPreNZ_len)/1e6);
Rprintf("ppiVoxelPlacePreNZ_VoxelTransition elements: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlacesPreNZ_len, sizeof(int));
#endif
int *piThisSNZ = (int *) R_alloc(iVoxelPlacesSNZ_len, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelPlacesSNZ_len)/1e6);
Rprintf("ppiVoxelPlaceSNZ_VoxelTransition elements: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelPlacesSNZ_len, sizeof(int));
#endif
for(iVoxelPlace = 0; iVoxelPlace < iVoxelPlaces; iVoxelPlace++) {
ppiVoxelPlacePreNZ_VoxelTransition[iVoxelPlace] = piThisPreNZ;
ppiVoxelPlaceSNZ_VoxelTransition[iVoxelPlace] = piThisSNZ;
piThisPreNZ += piVoxelPlacePreNZ_len[iVoxelPlace];
piThisSNZ += piVoxelPlaceSNZ_len[iVoxelPlace];
piVoxelPlacePreNZ_len[iVoxelPlace] = 0;
piVoxelPlaceSNZ_len[iVoxelPlace] = 0;
}
piThisPreNZ = (int *) R_alloc(iVoxelTransitionsPreNZ_len, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelTransitionsPreNZ_len)/1e6);
Rprintf("piPreNZ_VoxelPlace elements: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitionsPreNZ_len, sizeof(int));
#endif
char *pcThisPreNZ_value = (char *) R_alloc(iVoxelTransitionsPreNZ_len, sizeof(char));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(char) * iVoxelTransitionsPreNZ_len)/1e6);
Rprintf("pcPreNZ_value total: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitionsPreNZ_len, sizeof(char));
#endif
piThisSNZ = (int *) R_alloc(iVoxelTransitionsSNZ_len, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelTransitionsSNZ_len)/1e6);
Rprintf("piSNZ_VoxelPlace elements: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitionsSNZ_len, sizeof(int));
#endif
char *pcThisSNZ_value = (char *) R_alloc(iVoxelTransitionsSNZ_len, sizeof(char));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(char) * iVoxelTransitionsSNZ_len)/1e6);
Rprintf("pcSNZ_value total: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitionsSNZ_len, sizeof(char));
#endif
iVoxelCurrentTransition = 0;
iVoxelPlace_offset = 0;
// Third:
for (iVoxelFamily = 0; iVoxelFamily < iVoxelFamily_len; iVoxelFamily++) {
//Rprintf("Voxel group: %d\n", iVoxelFamily);
SEXP this_VoxelFamily = VECTOR_ELT(VoxelFamily, iVoxelFamily);
sexpTmp = getListElement(this_VoxelFamily, "pre");
SEXP pre;
PROTECT(pre = coerceVector(sexpTmp, INTSXP));
int *piPre = INTEGER(pre);
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iOrigTransitions = piTmp[0], iOrigPlaces = piTmp[1];
//Rprintf("pre: rows: %d\tcols: %d\n", iOrigTransitions, iOrigPlaces);
sexpTmp = getListElement(this_VoxelFamily, "Voxels");
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iVoxels = INTEGER(sexpTmp)[0];
UNPROTECT_PTR(sexpTmp);
sexpTmp = getListElement(this_VoxelFamily, "post");
SEXP post;
PROTECT(post = coerceVector(sexpTmp, INTSXP));
int *piPost = INTEGER(post);
piPreNZ_VoxelPlace = (int *) Realloc(piPreNZ_VoxelPlace, iOrigPlaces, int);
pcPreNZ_value = (char *) Realloc(pcPreNZ_value, iOrigPlaces, char);
piSNZ_VoxelPlace = (int *) Realloc(piSNZ_VoxelPlace, iOrigPlaces, int);
pcSNZ_value = (char *) Realloc(pcSNZ_value, iOrigPlaces, char);
sexpTmp = getListElement(this_VoxelFamily, "TransitionInVoxelSubset");
SEXP TransitionInVoxelSubset;
PROTECT(TransitionInVoxelSubset = coerceVector(sexpTmp, INTSXP));
int *piTransitionInVoxelSubset = INTEGER(TransitionInVoxelSubset);
SEXP VoxelSubset;
VoxelSubset = getListElement(this_VoxelFamily, "VoxelSubset");
SEXP h;
h = getListElement(this_VoxelFamily, "h");
int iOrigTransition, iOrigPlace;
for (iOrigTransition = 0; iOrigTransition < iOrigTransitions; iOrigTransition++) {
char cPreNZ_len = 0;
char cSNZ_len = 0;
for (iOrigPlace = 0; iOrigPlace < iOrigPlaces; iOrigPlace++) {
if ((cTmp = piPre[iOrigTransition + iOrigTransitions * iOrigPlace])) {
piPreNZ_VoxelPlace[(int) cPreNZ_len] = iOrigPlace;
pcPreNZ_value[(int) cPreNZ_len++] = cTmp;
}
if ((cTmp = piPost[iOrigTransition + iOrigTransitions * iOrigPlace] -
piPre[iOrigTransition + iOrigTransitions * iOrigPlace])) {
piSNZ_VoxelPlace[(int) cSNZ_len] = iOrigPlace;
pcSNZ_value[(int) cSNZ_len++] = cTmp;
}
}
// For the current transition, find the subset of Voxels in which it applies.
sexpTmp = VECTOR_ELT(VoxelSubset, piTransitionInVoxelSubset[iOrigTransition] - 1);
int iVoxelSubset_len = length(sexpTmp);
SEXP VoxelSubset;
PROTECT(VoxelSubset = coerceVector(sexpTmp, INTSXP));
int *piVoxelSubset = INTEGER(VoxelSubset); // Vector with subset of Voxels
int iVoxelSubset;
// Create a corresponding transition on each of the Voxels of this subset.
for (iVoxelSubset = 0; iVoxelSubset < iVoxelSubset_len; iVoxelSubset++) {
int iVoxel = piVoxelSubset[iVoxelSubset] - 1;
int iLocalVoxelPlace_offset = iVoxelPlace_offset + iOrigPlaces * iVoxel;
//Rprintf("Voxel: %d\tTransition: %d\tVoxelTransition: %d\n", iVoxel, iOrigTransition, iVoxelCurrentTransition);
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelCurrentTransition];
if (inherits(sexpTmp = VECTOR_ELT(h, iOrigTransition), "NativeSymbol")) {
pThisVoxelTransition->h.pDL_FUNC = (void *) R_ExternalPtrAddr(sexpTmp);
pThisVoxelTransition->cType_of_h = HZ_CFUNCTION;
} else if (isNumeric(sexpTmp)){
pThisVoxelTransition->h.dTrCnt = REAL(sexpTmp)[0];
pThisVoxelTransition->cType_of_h = HZ_DOUBLE;
} else if (isFunction(sexpTmp)) {
//SET_VECTOR_ELT(sexpFunction, iOrigTransition, lang1(sexpTmp));
//piHzType[iOrigTransition] = HZ_RFUNCTION;
;
} else {
error("Unrecognized transition function type\n");
}
pThisVoxelTransition->iLocalVoxelPlace_offset = iLocalVoxelPlace_offset;
pThisVoxelTransition->cPreNZ_len = cPreNZ_len;
pThisVoxelTransition->piPreNZ_VoxelPlace = piThisPreNZ;
piThisPreNZ += cPreNZ_len;
pThisVoxelTransition->pcPreNZ_value = pcThisPreNZ_value;
pcThisPreNZ_value += cPreNZ_len;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < cPreNZ_len; iVoxelPlaceIdx++) {
iVoxelPlace = iLocalVoxelPlace_offset + piPreNZ_VoxelPlace[iVoxelPlaceIdx];
pThisVoxelTransition->piPreNZ_VoxelPlace[iVoxelPlaceIdx] = iVoxelPlace;
pThisVoxelTransition->pcPreNZ_value[iVoxelPlaceIdx] = pcPreNZ_value[iVoxelPlaceIdx];
ppiVoxelPlacePreNZ_VoxelTransition[iVoxelPlace][piVoxelPlacePreNZ_len[iVoxelPlace]++] = iVoxelCurrentTransition;
}
pThisVoxelTransition->cSNZ_len = cSNZ_len;
pThisVoxelTransition->piSNZ_VoxelPlace = piThisSNZ;
piThisSNZ += cSNZ_len;
pThisVoxelTransition->pcSNZ_value = pcThisSNZ_value;
pcThisSNZ_value += cSNZ_len;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < cSNZ_len; iVoxelPlaceIdx++) {
iVoxelPlace = iLocalVoxelPlace_offset + piSNZ_VoxelPlace[iVoxelPlaceIdx];
pThisVoxelTransition->piSNZ_VoxelPlace[iVoxelPlaceIdx] = iVoxelPlace;
pThisVoxelTransition->pcSNZ_value[iVoxelPlaceIdx] = pcSNZ_value[iVoxelPlaceIdx];
ppiVoxelPlaceSNZ_VoxelTransition[iVoxelPlace][piVoxelPlaceSNZ_len[iVoxelPlace]++] = iVoxelCurrentTransition;
}
iVoxelCurrentTransition++;
}
UNPROTECT_PTR(VoxelSubset);
}
UNPROTECT_PTR(TransitionInVoxelSubset);
UNPROTECT_PTR(post);
UNPROTECT_PTR(pre);
// Manage jumps.
SEXP JumpingPlace;
if ((JumpingPlace = getListElement(this_VoxelFamily, "JumpingPlace")) != R_NilValue) {
SEXP JumpingPattern;
JumpingPattern = getListElement(this_VoxelFamily, "JumpingPattern");
int iJumpingPlace;
for (iJumpingPlace = 0; iJumpingPlace < length(JumpingPlace); iJumpingPlace++) {
SEXP this_JumpingPlace = VECTOR_ELT(JumpingPlace, iJumpingPlace);
sexpTmp = VECTOR_ELT(this_JumpingPlace, 0);
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iOrigPlace = INTEGER(sexpTmp)[0] - 1;
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(this_JumpingPlace, 1);
PROTECT(sexpTmp = coerceVector(sexpTmp, REALSXP));
double dH = REAL(sexpTmp)[0];
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(this_JumpingPlace, 2);
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iJumpingPattern = INTEGER(sexpTmp)[0] - 1;
UNPROTECT_PTR(sexpTmp);
sexpTmp = VECTOR_ELT(JumpingPattern, iJumpingPattern);
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iJumps = piTmp[0];
SEXP jump;
PROTECT(jump = coerceVector(sexpTmp, INTSXP));
int *piJump = INTEGER(jump);
int iJump;
for (iJump = 0; iJump < iJumps; iJump++) {
int iFromVoxelPlace = iVoxelPlace_offset + iOrigPlaces * (piJump[iJump] - 1) + iOrigPlace;
int iToVoxelPlace = iVoxelPlace_offset + iOrigPlaces * (piJump[iJump + iJumps] - 1) + iOrigPlace;
char cToVoxelPlacePost_value = piJump[iJump + iJumps*2];
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelCurrentTransition];
pThisVoxelTransition->cPreNZ_len = 1;
pThisVoxelTransition->piPreNZ_VoxelPlace = piThisPreNZ;
piThisPreNZ++;
pThisVoxelTransition->pcPreNZ_value = pcThisPreNZ_value;
pcThisPreNZ_value++;
pThisVoxelTransition->h.dTrCnt = dH;
pThisVoxelTransition->cType_of_h = HZ_DOUBLE;
pThisVoxelTransition->iLocalVoxelPlace_offset = iVoxelPlace_offset + iOrigPlaces * (piJump[iJump] - 1);
pThisVoxelTransition->piPreNZ_VoxelPlace[0] = iFromVoxelPlace;
pThisVoxelTransition->pcPreNZ_value[0] = 1;
ppiVoxelPlacePreNZ_VoxelTransition[iFromVoxelPlace][piVoxelPlacePreNZ_len[iFromVoxelPlace]++] = iVoxelCurrentTransition;
pThisVoxelTransition->cSNZ_len = 1;
pThisVoxelTransition->piSNZ_VoxelPlace = piThisSNZ;
piThisSNZ++;
pThisVoxelTransition->pcSNZ_value = pcThisSNZ_value;
pcThisSNZ_value++;
pThisVoxelTransition->piSNZ_VoxelPlace[0] = iFromVoxelPlace;
pThisVoxelTransition->pcSNZ_value[0] = -1;
ppiVoxelPlaceSNZ_VoxelTransition[iFromVoxelPlace][piVoxelPlaceSNZ_len[iFromVoxelPlace]++] = iVoxelCurrentTransition;
if (cToVoxelPlacePost_value) {
pThisVoxelTransition->cSNZ_len++;
piThisSNZ++;
pcThisSNZ_value++;
pThisVoxelTransition->piSNZ_VoxelPlace[1] = iToVoxelPlace;
pThisVoxelTransition->pcSNZ_value[1] = cToVoxelPlacePost_value;
ppiVoxelPlaceSNZ_VoxelTransition[iToVoxelPlace][piVoxelPlaceSNZ_len[iToVoxelPlace]++] = iVoxelCurrentTransition;
}
//Rprintf("%d\t%d\n", iFromVoxelPlace, iToVoxelPlace);
iVoxelCurrentTransition++;
}
UNPROTECT_PTR(jump);
}
}
iVoxelPlace_offset += iOrigPlaces * iVoxels;
}
Free(piPreNZ_VoxelPlace);
Free(pcPreNZ_value);
Free(piSNZ_VoxelPlace);
Free(pcSNZ_value);
// Rprintf("VoxelTransitions: %d\tVoxelCurrentTransition: %d\n", iVoxelTransitions, iVoxelCurrentTransition);
// For the initial calculation of all hazards...
int *piHazardsToMod_start = (int *) R_alloc(iVoxelTransitions, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iVoxelTransitions)/1e6);
Rprintf("piHazardsToMod_start: %.2f MB (%d x %d bytes)\n", dThisMem, iVoxelTransitions, sizeof(int));
#endif
int iHazardsToMod_count = 0;
// Identify the transitions that need the hazards recalculated
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelTransition];
int iHazardToCompTot = 0;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < pThisVoxelTransition->cSNZ_len; iVoxelPlaceIdx++) {
int iVoxelPlace = pThisVoxelTransition->piSNZ_VoxelPlace[iVoxelPlaceIdx];
int iVoxelTransitionIdx2, iVoxelTransition2;
for(iVoxelTransitionIdx2 = 0; iVoxelTransitionIdx2 < piVoxelPlacePreNZ_len[iVoxelPlace]; iVoxelTransitionIdx2++) {
iVoxelTransition2 = ppiVoxelPlacePreNZ_VoxelTransition[iVoxelPlace][iVoxelTransitionIdx2];
int iAddThis = TRUE;
int iVoxelTransitionIdx3;
for (iVoxelTransitionIdx3 = 0; iVoxelTransitionIdx3 < iHazardToCompTot; iVoxelTransitionIdx3++) {
if(piHazardsToMod_start[iVoxelTransitionIdx3] == iVoxelTransition2) {
iAddThis = FALSE;
break;
}
}
if (iAddThis)
piHazardsToMod_start[iHazardToCompTot++] = iVoxelTransition2;
}
}
iHazardsToMod_count += iHazardToCompTot;
}
int *piThisHazardsToMod = (int *) R_alloc(iHazardsToMod_count, sizeof(int));
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iHazardsToMod_count)/1e6);
Rprintf("piHazardsToMod elements: %.2f MB (%d x %d bytes)\n", dThisMem, iHazardsToMod_count, sizeof(int));
#endif
// Identify the transitions that need the hazards recalculated
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelTransition];
int iHazardToCompTot = 0;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < pThisVoxelTransition->cSNZ_len; iVoxelPlaceIdx++) {
int iVoxelPlace = pThisVoxelTransition->piSNZ_VoxelPlace[iVoxelPlaceIdx];
int iVoxelTransitionIdx2, iVoxelTransition2;
for(iVoxelTransitionIdx2 = 0; iVoxelTransitionIdx2 < piVoxelPlacePreNZ_len[iVoxelPlace]; iVoxelTransitionIdx2++) {
iVoxelTransition2 = ppiVoxelPlacePreNZ_VoxelTransition[iVoxelPlace][iVoxelTransitionIdx2];
int iAddThis = TRUE;
int iVoxelTransitionIdx3;
for (iVoxelTransitionIdx3 = 0; iVoxelTransitionIdx3 < iHazardToCompTot; iVoxelTransitionIdx3++) {
if(piHazardsToMod_start[iVoxelTransitionIdx3] == iVoxelTransition2) {
iAddThis = FALSE;
break;
}
}
if (iAddThis)
piHazardsToMod_start[iHazardToCompTot++] = iVoxelTransition2;
}
}
pThisVoxelTransition->iHazardsToMod_count = iHazardToCompTot;
pThisVoxelTransition->piHazardsToMod = piThisHazardsToMod;
piThisHazardsToMod += iHazardToCompTot;
for (iVoxelTransitionIdx = 0; iVoxelTransitionIdx < iHazardToCompTot; iVoxelTransitionIdx++) {
pThisVoxelTransition->piHazardsToMod[iVoxelTransitionIdx] = piHazardsToMod_start[iVoxelTransitionIdx];
}
}
// For the initial calculation of all hazards...
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
piHazardsToMod_start[iVoxelTransition] = iVoxelTransition;
}
SEXP sexpCrntMarking;
PROTECT(sexpCrntMarking = allocVector(REALSXP, iVoxelPlaces));
double *pdCrntMarking = REAL(sexpCrntMarking);
double dDelta = *REAL(delta);
int iTotalSteps, iSectionSteps;
double dT = 0;
void *pCManage_time = 0;
SEXP sexpRManage_time = 0;
if (inherits(T, "NativeSymbol")) {
pCManage_time = (void *) R_ExternalPtrAddr(T);
dT = ((double(*)(double, double *)) pCManage_time)(-1, pdCrntMarking);
} else if (isNumeric(T)){
dT = *REAL(T);
} else if (isFunction(T)) {
PROTECT(sexpRManage_time = lang1(T));
defineVar(install("y"), sexpCrntMarking, rho);
PROTECT(sexpTmp = allocVector(REALSXP, 1));
*REAL(sexpTmp) = -1;
defineVar(install("StartTime"), sexpTmp, rho);
UNPROTECT_PTR(sexpTmp);
dT = *REAL(VECTOR_ELT(eval(sexpRManage_time, rho),0));
} else {
error("Unrecognized time function type\n");
}
iTotalSteps = iSectionSteps = (int)(dT / dDelta) + 1;
#ifdef RB_TIME
iTotalStepsOver50 = iTotalSteps/50;
iTotalStepsOver10 = iTotalSteps/10;
#endif
double *pdHazard = (double *) R_alloc(iVoxelTransitions, sizeof(double));
int iRun, iRuns = *INTEGER(runs);
SEXP sexpRun;
PROTECT(sexpRun = allocVector(VECSXP, iRuns));
int iTotalUsedRandomNumbers = 0;
// DiscTime Vector
SEXP sexpD_time;
PROTECT(sexpD_time = allocVector(REALSXP, iTotalSteps));
double *pdDiscTime = REAL(sexpD_time);
int k;
double dTmp = 0;
for (k = 0; k < iTotalSteps; k++) {
pdDiscTime[k] = dTmp;
dTmp += dDelta;
}
SEXP sexpMarkingRowNames;
PROTECT(sexpMarkingRowNames = allocVector(INTSXP, iTotalSteps));
piTmp = INTEGER(sexpMarkingRowNames);
for (k = 0; k < iTotalSteps; k++)
piTmp[k] = k+1;
double **ppdMarking = (double **) R_alloc(iVoxelPlaces, sizeof(double *));
int iLevels = 7;
int iGroups = pow(2, iLevels - 1);
// Group holding the transitions that lie between boundaries
int **ppiGroup = (int **) R_alloc(iGroups, sizeof(int *));
// Number of transition each group has
int *piGroupElm = (int *) R_alloc(iGroups, sizeof(int));
// Total propensity hazard for each group
int *piTotGroupTransitions = (int *) R_alloc(iGroups, sizeof(int));
int iGroup, iGroupsWithElements = 0;
for (iGroup = 0; iGroup < iGroups; iGroup++) {
//ppiGroup[iGroup] = NULL;
ppiGroup[iGroup] = (int *) R_alloc(iVoxelTransitions, sizeof(int));
}
node **ppnodeLevel = (node **) R_alloc(iLevels, sizeof(node *));
int iLevel, iNode;
int iNodesPerLevel = 1;
for (iLevel = 0; iLevel < iLevels; iLevel++) {
ppnodeLevel[iLevel] = (node *) R_alloc(iNodesPerLevel, sizeof(node));
iNodesPerLevel *= 2;
}
node *pnodeRoot = &ppnodeLevel[0][0];
pnodeRoot->parent = 0;
node *pnodeGroup = ppnodeLevel[iLevels-1];
iNodesPerLevel = 1;
for (iLevel = 0; iLevel < iLevels; iLevel++) {
for (iNode = 0; iNode < iNodesPerLevel; iNode++) {
if (iLevel < iLevels-1) {
ppnodeLevel[iLevel][iNode].iGroup = -1;
ppnodeLevel[iLevel][iNode].left = &ppnodeLevel[iLevel+1][iNode*2];
ppnodeLevel[iLevel][iNode].right = &ppnodeLevel[iLevel+1][iNode*2+1];
ppnodeLevel[iLevel+1][iNode*2].parent = ppnodeLevel[iLevel+1][iNode*2+1].parent =
&ppnodeLevel[iLevel][iNode];
} else {
ppnodeLevel[iLevel][iNode].iGroup = iNode;
ppnodeLevel[iLevel][iNode].left = ppnodeLevel[iLevel][iNode].right = 0;
}
}
iNodesPerLevel *= 2;
}
double dNewHazard = 0;
// Find minimum propensity
double dMinHazard = DBL_MAX;
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelTransition];
switch(pThisVoxelTransition->cType_of_h) {
case HZ_DOUBLE:
dNewHazard = pThisVoxelTransition->h.dTrCnt;
int iVoxelPlaceIdx;
for(iVoxelPlaceIdx = 0;
iVoxelPlaceIdx < pThisVoxelTransition->cPreNZ_len;
iVoxelPlaceIdx++) {
for (k = 0; k < pThisVoxelTransition->pcPreNZ_value[iVoxelPlaceIdx]; k++)
dNewHazard *= (pThisVoxelTransition->pcPreNZ_value[iVoxelPlaceIdx] - k) / (double)(k+1);
}
if (dNewHazard > 0 && dNewHazard < dMinHazard)
dMinHazard = dNewHazard;
break;
case HZ_CFUNCTION:
break;
case HZ_RFUNCTION:
break;
}
}
#ifdef RB_TIME
c1 = clock();
Rprintf ("Elapsed CPU time: %f\n", (float) (c1 - c0)/CLOCKS_PER_SEC);
c0 = c1;
#endif
GetRNGstate();
for (iRun = 0; iRun < iRuns; iRun++) {
int iUsedRandomNumbers = 0;
Rprintf("%d ", iRun+1);
// Totals for kind of transition vector
SEXP sexpTotXTransition;
PROTECT(sexpTotXTransition = allocVector(INTSXP, iVoxelTransitions));
int *piTotTransitions = INTEGER(sexpTotXTransition);
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
piTotTransitions[iVoxelTransition] = 0;
}
SEXP sexpMarking;
PROTECT(sexpMarking = allocVector(VECSXP, iVoxelPlaces));
// NOTE: Need to address next line (expansion of place names!)
//setAttrib(sexpMarking, R_NamesSymbol, place);
//setAttrib(sexpMarking, R_RowNamesSymbol, sexpMarkingRowNames);
//setAttrib(sexpMarking, R_ClassSymbol, ScalarString(mkChar("data.frame")));
// Setup initial state
int iVoxelCurrentPlace = 0;
for (iVoxelFamily = 0; iVoxelFamily < iVoxelFamily_len; iVoxelFamily++) {
SEXP this_VoxelFamily = VECTOR_ELT(VoxelFamily, iVoxelFamily);
sexpTmp = getListElement(this_VoxelFamily, "pre");
piTmp = INTEGER(getAttrib(sexpTmp, R_DimSymbol));
int iOrigPlaces = piTmp[1];
sexpTmp = getListElement(this_VoxelFamily, "Voxels");
PROTECT(sexpTmp = coerceVector(sexpTmp, INTSXP));
int iVoxels = INTEGER(sexpTmp)[0];
UNPROTECT_PTR(sexpTmp);
if ((sexpTmp = getListElement(this_VoxelFamily, "M")) == R_NilValue)
error("M not provided!\n");
SEXP M;
PROTECT(M = coerceVector(sexpTmp, REALSXP));
double *pdM = REAL(M);
Rprintf("\n");
int iVoxel, iOrigPlace;
for (iVoxel = 0; iVoxel < iVoxels; iVoxel++) {
for (iOrigPlace = 0; iOrigPlace < iOrigPlaces; iOrigPlace++) {
SET_VECTOR_ELT(sexpMarking, iVoxelCurrentPlace, sexpTmp = allocVector(REALSXP, iTotalSteps));
ppdMarking[iVoxelCurrentPlace] = REAL(sexpTmp);
pdCrntMarking[iVoxelCurrentPlace] = pdM[iVoxel + iVoxels * iOrigPlace];
iVoxelCurrentPlace++;
}
}
UNPROTECT_PTR(M);
}
for(iVoxelTransition = 0; iVoxelTransition < iVoxelTransitions; iVoxelTransition++) {
pdHazard[iVoxelTransition] = 0;
pVoxelTransition[iVoxelTransition].cInGroup = -1;
}
for (iGroup = 0; iGroup < iGroups; iGroup++) {
piGroupElm[iGroup] = 0;
piTotGroupTransitions[iGroup] = 0;
}
iNodesPerLevel = 1;
for (iLevel = 0; iLevel < iLevels; iLevel++) {
for (iNode = 0; iNode < iNodesPerLevel; iNode++) {
ppnodeLevel[iLevel][iNode].dPartialAcumHazard = 0;
}
iNodesPerLevel *= 2;
}
node *pnode;
double dTime = 0, dTarget = 0;
int iTotTransitions = 0;
int iStep = 0;
int iInterruptCnt = 10000000;
double dNewHazard = 0;
#ifdef RB_TIME
clkStep = clock();
iElSteps = 0;
#endif
do {
if (pCManage_time || sexpRManage_time) {
double dEnd = 0;
if (pCManage_time) {
dEnd = ((double(*)(double, double *)) pCManage_time)(dTarget, pdCrntMarking);
} else {
defineVar(install("y"), sexpCrntMarking, rho);
PROTECT(sexpTmp = allocVector(REALSXP, 1));
*REAL(sexpTmp) = dTarget;
defineVar(install("StartTime"), sexpTmp, rho);
UNPROTECT_PTR(sexpTmp);
sexpTmp = eval(sexpRManage_time, rho);
dEnd = *REAL(VECTOR_ELT(sexpTmp,0));
for(iVoxelPlace = 0; iVoxelPlace < iVoxelPlaces; iVoxelPlace++) {
pdCrntMarking[iVoxelPlace] = REAL(VECTOR_ELT(sexpTmp,1))[iVoxelPlace];
}
}
iSectionSteps = (int)(dEnd / dDelta) + 1;
}
for(iVoxelPlace = 0; iVoxelPlace < iVoxelPlaces; iVoxelPlace++) {
ppdMarking[iVoxelPlace][iStep] = pdCrntMarking[iVoxelPlace];
}
dTime = dTarget;
dTarget += dDelta;
// For the calculation of all hazards...
int *piHazardsToMod = piHazardsToMod_start;
int iHazardsToMod_count = iVoxelTransitions;
do {
// Get hazards only for the transitions associated with
// places whose quantities changed in the last step.
for(iVoxelTransitionIdx = 0; iVoxelTransitionIdx < iHazardsToMod_count; iVoxelTransitionIdx++) {
iVoxelTransition = piHazardsToMod[iVoxelTransitionIdx];
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelTransition];
switch(pThisVoxelTransition->cType_of_h) {
case HZ_DOUBLE:
dNewHazard = pThisVoxelTransition->h.dTrCnt;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < pThisVoxelTransition->cPreNZ_len; iVoxelPlaceIdx++) {
for (k = 0; k < pThisVoxelTransition->pcPreNZ_value[iVoxelPlaceIdx]; k++)
dNewHazard *= (pdCrntMarking[pThisVoxelTransition->piPreNZ_VoxelPlace[iVoxelPlaceIdx]] - k) / (double)(k+1);
}
break;
case HZ_CFUNCTION:
dNewHazard = ((double(*)(double, double *)) pThisVoxelTransition->h.pDL_FUNC)(dTime, pdCrntMarking + pThisVoxelTransition->iLocalVoxelPlace_offset);
break;
case HZ_RFUNCTION:
defineVar(install("y"), sexpCrntMarking, rho);
//dNewHazard = REAL(eval(VECTOR_ELT(sexpFunction, iVoxelTransition), rho))[0];
break;
}
double dDeltaHazard;
frexp(dNewHazard/dMinHazard, &iGroup);
if (iGroup-- > 0) {
// Transition belongs to a group
/*
// Check if the group has been initialized
if (!ppiGroup[iGroup]) {
ppiGroup[iGroup] = (int *) R_alloc(iVoxelTransitions, sizeof(int));
iGroupsWithElements++;
}
*/
if (iGroup == pThisVoxelTransition->cInGroup) {
// Transitions will stay in same group as it was
dDeltaHazard = dNewHazard - pdHazard[iVoxelTransition];
pnode = &pnodeGroup[iGroup];
do {
pnode->dPartialAcumHazard += dDeltaHazard;
} while ((pnode = pnode->parent));
} else if (pThisVoxelTransition->cInGroup != -1) {
// Transition was in another group and needs to be moved to the new one
int iOldGroup = pThisVoxelTransition->cInGroup;
int iOldPositionInGroup = pThisVoxelTransition->iPositionInGroup;
dDeltaHazard = -pdHazard[iVoxelTransition];
pnode = &pnodeGroup[iOldGroup];
do {
pnode->dPartialAcumHazard += dDeltaHazard;
} while ((pnode = pnode->parent));
piGroupElm[iOldGroup]--; // Old group will have one less element
// Now, piGroupElm[iOldGroup] is the index to last transition in group
if (iOldPositionInGroup != piGroupElm[iOldGroup]) {
// Transition is not the last in group,
// put the last transition in place of the one to be removed
ppiGroup[iOldGroup][iOldPositionInGroup] =
ppiGroup[iOldGroup][piGroupElm[iOldGroup]];
// Update position of previous last transition in group
pVoxelTransition[ppiGroup[iOldGroup][iOldPositionInGroup]].iPositionInGroup =
iOldPositionInGroup;
}
dDeltaHazard = dNewHazard;
pnode = &pnodeGroup[iGroup];
do {
pnode->dPartialAcumHazard += dDeltaHazard;
} while ((pnode = pnode->parent));
pThisVoxelTransition->cInGroup = iGroup;
pThisVoxelTransition->iPositionInGroup = piGroupElm[iGroup];
ppiGroup[iGroup][piGroupElm[iGroup]++] = iVoxelTransition;
} else if (pThisVoxelTransition->cInGroup == -1) { // Transition was in no group
dDeltaHazard = dNewHazard;
pnode = &pnodeGroup[iGroup];
do {
pnode->dPartialAcumHazard += dDeltaHazard;
} while ((pnode = pnode->parent));
pThisVoxelTransition->cInGroup = iGroup;
pThisVoxelTransition->iPositionInGroup = piGroupElm[iGroup];
ppiGroup[iGroup][piGroupElm[iGroup]++] = iVoxelTransition;
} else {
error("ERROR: Option not considered 1\n");
}
} else if (pThisVoxelTransition->cInGroup != -1) {
// Transition will not belong to any group and needs to be removed from old
int iOldGroup = pThisVoxelTransition->cInGroup;
int iOldPositionInGroup = pThisVoxelTransition->iPositionInGroup;
dDeltaHazard = -pdHazard[iVoxelTransition];
pnode = &pnodeGroup[iOldGroup];
do {
pnode->dPartialAcumHazard += dDeltaHazard;
} while ((pnode = pnode->parent));
piGroupElm[iOldGroup]--; // Old group will have one less element
// Now, piGroupElm[iOldGroup] is the index to last transition in group
if (iOldPositionInGroup != piGroupElm[iOldGroup]) {
// Transition is not the last in group,
// put the last transition in place of the one to be removed
ppiGroup[iOldGroup][iOldPositionInGroup] =
ppiGroup[iOldGroup][piGroupElm[iOldGroup]];
// Update position of previous last transition in group
pVoxelTransition[ppiGroup[iOldGroup][iOldPositionInGroup]].iPositionInGroup =
iOldPositionInGroup;
}
pThisVoxelTransition->cInGroup = -1;
}
pdHazard[iVoxelTransition] = dNewHazard;
}
// Get Time to transition
if (pnodeRoot->dPartialAcumHazard > 0) {
dTime += exp_rand() / pnodeRoot->dPartialAcumHazard;
//dTime -= log(unif_rand()) / pnodeRoot->dPartialAcumHazard;
iUsedRandomNumbers++;
} else
dTime = dT + 1;
while (dTime >= dTarget) {
++iStep;
// Update the state for the fixed incremented time.
for(iVoxelPlace = 0; iVoxelPlace < iVoxelPlaces; iVoxelPlace++)
ppdMarking[iVoxelPlace][iStep] = pdCrntMarking[iVoxelPlace];
if (iStep == iSectionSteps - 1)
goto EXIT_LOOP;
#ifdef RB_TIME
iElSteps++;
if ((iStep % iTotalStepsOver50) == 0) {
double clkTmp2 = clock();
double dSecondsToFinish = (double) (clkTmp2 - clkStep)/CLOCKS_PER_SEC/iElSteps*(iTotalSteps-iStep);
Rprintf ("\t");
PrintfTime(dSecondsToFinish);
clkStep = clkTmp2;
iElSteps = 0;
}
if ((iStep % iTotalStepsOver10) == 0) {
Rprintf(".\n");
}
#endif
dTarget += dDelta;
// Force check if user interrupted
iInterruptCnt = 1;
}
if (! --iInterruptCnt) {
// Allow user interruption
R_CheckUserInterrupt();
iInterruptCnt = 10000000;
}
do {
// Find group containing firing transition
double dRnd = unif_rand() * pnodeRoot->dPartialAcumHazard;
iUsedRandomNumbers++;
pnode = pnodeRoot;
do {
if (dRnd < pnode->left->dPartialAcumHazard) {
pnode = pnode->left;
} else {
dRnd -= pnode->left->dPartialAcumHazard;
pnode = pnode->right;
}
} while (pnode->left);
// Next check is because
// once in a while it is generated a number that goes past
// the last group or selects a group with zero elements
// due to accumulated truncation errors.
// Discard this random number and try again.
} while (piGroupElm[iGroup = pnode->iGroup] == 0);
int iGroupElm = piGroupElm[iGroup];
double dMaxInGroup = dMinHazard * pow(2, iGroup + 1);
// Find transition in group
while (1) {
if (! --iInterruptCnt) {
// Allow user interruption
R_CheckUserInterrupt();
iInterruptCnt = 10000000;
}
int iVoxelTransitionIdx = (int) (unif_rand() * iGroupElm);
iUsedRandomNumbers++;
iVoxelTransition = ppiGroup[iGroup][iVoxelTransitionIdx];
iUsedRandomNumbers++;
if (pdHazard[iVoxelTransition] > unif_rand() * dMaxInGroup) {
piTotTransitions[iVoxelTransition]++;
transition_el *pThisVoxelTransition = &pVoxelTransition[iVoxelTransition];
piHazardsToMod = pThisVoxelTransition->piHazardsToMod;
iHazardsToMod_count = pThisVoxelTransition->iHazardsToMod_count;
for(iVoxelPlaceIdx = 0; iVoxelPlaceIdx < pThisVoxelTransition->cSNZ_len; iVoxelPlaceIdx++) {
// Update the state
pdCrntMarking[pThisVoxelTransition->piSNZ_VoxelPlace[iVoxelPlaceIdx]] += pThisVoxelTransition->pcSNZ_value[iVoxelPlaceIdx];
}
break;
}
}
++iTotTransitions;
} while (TRUE);
EXIT_LOOP:;
Rprintf("|\n");
} while (iSectionSteps < iTotalSteps);
iTotalUsedRandomNumbers += iUsedRandomNumbers;
Rprintf("\t%d\t%d\t%d", iTotTransitions, iUsedRandomNumbers, iTotalUsedRandomNumbers);
#ifdef RB_SUBTIME
c1 = clock();
Rprintf ("\t To go: ");
PrintfTime((double) (c1 - c0)/CLOCKS_PER_SEC/(iRun+1)*(iRuns-iRun-1));
#endif
Rprintf ("\n");
SEXP sexpTotTransitions;
PROTECT(sexpTotTransitions = allocVector(INTSXP, 1));
INTEGER(sexpTotTransitions)[0] = iTotTransitions;
SEXP sexpThisRun;
PROTECT(sexpThisRun = allocVector(VECSXP, 3));
SET_VECTOR_ELT(sexpThisRun, 0, sexpMarking);
UNPROTECT_PTR(sexpMarking);
SET_VECTOR_ELT(sexpThisRun, 1, sexpTotXTransition);
UNPROTECT_PTR(sexpTotXTransition);
SET_VECTOR_ELT(sexpThisRun, 2, sexpTotTransitions);
UNPROTECT_PTR(sexpTotTransitions);
SEXP sexpNames;
PROTECT(sexpNames = allocVector(VECSXP, 3));
SET_VECTOR_ELT(sexpNames, 0, mkChar("M"));
SET_VECTOR_ELT(sexpNames, 1, mkChar("transitions"));
SET_VECTOR_ELT(sexpNames, 2, mkChar("tot.transitions"));
setAttrib(sexpThisRun, R_NamesSymbol, sexpNames);
UNPROTECT_PTR(sexpNames);
SET_VECTOR_ELT(sexpRun, iRun, sexpThisRun);
UNPROTECT_PTR(sexpThisRun);
}
PutRNGstate();
#ifdef RB_MEMORY
dUsedMemAcum += (dThisMem = ((double) sizeof(int) * iGroupsWithElements * iVoxelTransitions)/1e6);
Rprintf("ppiGroup elements (%d/%d): %.2f MB (%d x %d bytes)\n", iGroupsWithElements, iGroups, dThisMem, iGroupsWithElements * iVoxelTransitions, sizeof(int));
#endif
#ifdef RB_MEMORY
Rprintf("Total Memory used: %.2f MB\n", dUsedMemAcum);
#endif
SEXP sexpAns;
PROTECT(sexpAns = allocVector(VECSXP, 4));
SET_VECTOR_ELT(sexpAns, 0, place);
SET_VECTOR_ELT(sexpAns, 1, transition);
SET_VECTOR_ELT(sexpAns, 2, sexpD_time);
UNPROTECT_PTR(sexpD_time);
SET_VECTOR_ELT(sexpAns, 3, sexpRun);
UNPROTECT_PTR(sexpRun);
SEXP sexpNames;
PROTECT(sexpNames = allocVector(VECSXP, 4));
SET_VECTOR_ELT(sexpNames, 0, mkChar("place"));
SET_VECTOR_ELT(sexpNames, 1, mkChar("transition"));
SET_VECTOR_ELT(sexpNames, 2, mkChar("dt"));
SET_VECTOR_ELT(sexpNames, 3, mkChar("run"));
setAttrib(sexpAns, R_NamesSymbol, sexpNames);
UNPROTECT_PTR(sexpNames);
#ifdef RB_TIME
c1 = clock();
double dCpuTime = (double) (c1 - c0)/CLOCKS_PER_SEC;
Rprintf ("Elapsed CPU time: ");
PrintfTime(dCpuTime);
Rprintf ("\t(%fs)\n", dCpuTime);
#endif
if (sexpRManage_time)
UNPROTECT_PTR(sexpRManage_time);
// UNPROTECT_PTR(sexpFunction);
UNPROTECT_PTR(sexpMarkingRowNames);
UNPROTECT_PTR(sexpCrntMarking);
UNPROTECT_PTR(sexpAns);
return(sexpAns);
}
