Diff of /src/jags.cc [c629bb] .. [fda875] Maximize Restore

  Switch to unified view

a/src/jags.cc b/src/jags.cc
...
...
52
    error("Invalid integer parameter");
52
    error("Invalid integer parameter");
53
    }
53
    }
54
    
54
    
55
    SEXP intarg;
55
    SEXP intarg;
56
    PROTECT(intarg = AS_INTEGER(arg));
56
    PROTECT(intarg = AS_INTEGER(arg));
57
    int i = INTEGER_POINTER(intarg)[0];
57
    int i = INTEGER(intarg)[0];
58
    UNPROTECT(1);
58
    UNPROTECT(1); //intarg
59
    return i;
59
    return i;
60
}
60
}
61
61
62
static char const *stringArg(SEXP arg, unsigned int i = 0)
62
static char const *stringArg(SEXP arg, unsigned int i = 0)
63
{
63
{
...
...
120
}
120
}
121
121
122
/* Write data from an R list into a JAGS data table */
122
/* Write data from an R list into a JAGS data table */
123
static void writeDataTable(SEXP data, map<string,SArray> &table)
123
static void writeDataTable(SEXP data, map<string,SArray> &table)
124
{
124
{
125
    SEXP names;
126
    PROTECT(names = getAttrib(data, R_NamesSymbol));
125
    SEXP names = getAttrib(data, R_NamesSymbol);
127
    if (!isNewList(data)) {
126
    if (!isNewList(data)) {
128
    error("data must be a list");
127
    error("data must be a list");
129
    }
128
    }
130
    if (length(names) != length(data)) {
129
    if (length(names) != length(data)) {
131
    error("data must be a named list");
130
    error("data must be a named list");
132
    }
131
    }
133
    int N = length(data);
134
132
135
    for (int i = 0; i < N; ++i) {
133
    for (int i = 0; i < length(data); ++i) {
136
  SEXP e, e2, dim;
134
  SEXP e;
137
    PROTECT(e = VECTOR_ELT(data, i));
135
    PROTECT(e = AS_NUMERIC(VECTOR_ELT(data, i)));
138
  PROTECT(dim = GET_DIM(e)); 
136
  if (length(e) > 0) {
139
  PROTECT(e2 = AS_NUMERIC(e));
137
140
    //Replace R missing values in e2 with JAGS missing values
138
        //Replace R missing values in e with JAGS missing values
141
  int elength = length(e2);
142
    for (int j = 0; j < elength; ++j) {
139
        for (int j = 0; j < length(e); ++j) {
143
        if (ISNA(NUMERIC_POINTER(e2)[j])) {
140
      if (ISNA(NUMERIC_POINTER(e)[j])) {
144
        NUMERIC_POINTER(e2)[j] = JAGS_NA;
141
            NUMERIC_POINTER(e)[j] = JAGS_NA;
142
      }
143
      }
145
        }
144
        
146
  }
147
148
    string ename = CHAR(STRING_ELT(names, i));
145
        string ename = CHAR(STRING_ELT(names, i));
149
146
150
  int ndim = length(dim);
147
      SEXP dim = getAttrib(VECTOR_ELT(data, i), R_DimSymbol); 
151
  if (ndim == 0) {
148
      if (dim == R_NilValue) {
152
      // Scalar or vector entry. Skip vectors of length zero
149
      // Scalar or vector entry.
153
      if (elength > 0) {
154
        SArray sarray(vector<unsigned int>(1, elength));
150
        SArray sarray(vector<unsigned int>(1, length(e)));
155
        setSArrayValue(sarray, e2);
151
        setSArrayValue(sarray, e);
156
        table.insert(pair<string,SArray>(ename, sarray));
152
        table.insert(pair<string,SArray>(ename, sarray));
157
        }
153
        }
158
  }
159
    else {
154
        else {
160
        // Array entry
155
      // Array entry
156
      int ndim = length(dim);
161
        vector<unsigned int> idim(ndim);
157
      vector<unsigned int> idim(ndim);
162
      SEXP dim2;
163
      PROTECT(dim2 = AS_INTEGER(dim));
164
        for (int j = 0; j < ndim; ++j) {
158
      for (int j = 0; j < ndim; ++j) {
165
        idim[j] = INTEGER_POINTER(dim2)[j];
159
            idim[j] = INTEGER(dim)[j];
166
      }
160
      }
167
      UNPROTECT(1);
168
        SArray sarray(idim);
161
      SArray sarray(idim);
169
        setSArrayValue(sarray, e2);
162
      setSArrayValue(sarray, e);
170
        table.insert(pair<string,SArray>(ename,sarray));
163
      table.insert(pair<string,SArray>(ename,sarray));
171
  }
172
  UNPROTECT(3);
173
    }
164
      }
165
  }
174
    UNPROTECT(1);
166
  UNPROTECT(1); //e
167
    }
175
}
168
}
176
169
177
static Range makeRange(SEXP lower, SEXP upper)
170
static Range makeRange(SEXP lower, SEXP upper)
178
{
171
{
179
    if (lower == R_NilValue || upper == R_NilValue) {
172
    if (lower == R_NilValue || upper == R_NilValue) {
...
...
188
    PROTECT(il = AS_INTEGER(lower));
181
    PROTECT(il = AS_INTEGER(lower));
189
    PROTECT(iu = AS_INTEGER(upper));
182
    PROTECT(iu = AS_INTEGER(upper));
190
    vector<int> lvec(n), uvec(n);
183
    vector<int> lvec(n), uvec(n);
191
    copy(INTEGER(il), INTEGER(il) + n, lvec.begin());
184
    copy(INTEGER(il), INTEGER(il) + n, lvec.begin());
192
    copy(INTEGER(iu), INTEGER(iu) + n, uvec.begin());
185
    copy(INTEGER(iu), INTEGER(iu) + n, uvec.begin());
193
    UNPROTECT(2);
186
    UNPROTECT(2); //il, iu
194
187
195
    Range r;
188
    Range r;
196
    try {
189
    try {
197
    r = Range(lvec, uvec);
190
    r = Range(lvec, uvec);
198
    }
191
    }
...
...
237
        vector<unsigned int> const &idim = p->second.dim(false);
230
        vector<unsigned int> const &idim = p->second.dim(false);
238
        unsigned int ndim = idim.size();
231
        unsigned int ndim = idim.size();
239
        SEXP dim;
232
        SEXP dim;
240
        PROTECT(dim = allocVector(INTSXP, ndim));
233
        PROTECT(dim = allocVector(INTSXP, ndim));
241
        for (unsigned int k = 0; k < ndim; ++k) {
234
        for (unsigned int k = 0; k < ndim; ++k) {
242
        INTEGER_POINTER(dim)[k] = idim[k];
235
        INTEGER(dim)[k] = idim[k];
243
        }
236
        }
244
237
245
        //Set names of the dimensions 
238
        //Set names of the dimensions 
246
        vector<string> const &names = p->second.dimNames();
239
        vector<string> const &names = p->second.dimNames();
247
        if (!names.empty()) {
240
        if (!names.empty()) {
...
...
546
    for (unsigned int n = 0; n < nchain; ++n) {
539
    for (unsigned int n = 0; n < nchain; ++n) {
547
        string srng;
540
        string srng;
548
        map<string,SArray> param_table;
541
        map<string,SArray> param_table;
549
        console->dumpState(param_table, srng, DUMP_PARAMETERS, n+1);
542
        console->dumpState(param_table, srng, DUMP_PARAMETERS, n+1);
550
        //Read the parameter values into an R list
543
        //Read the parameter values into an R list
551
        SEXP params, names;
544
        SEXP params;
552
        PROTECT(params = readDataTable(param_table));
545
        PROTECT(params = readDataTable(param_table));
553
        int nparam = length(params);
546
        int nparam = length(params);
554
        PROTECT(names = getAttrib(params, R_NamesSymbol));
547
        SEXP names = getAttrib(params, R_NamesSymbol);
555
        //Now we have to make a copy of the list with an extra element
548
        //Now we have to make a copy of the list with an extra element
556
        SEXP staten, namesn;
549
        SEXP staten, namesn;
557
        PROTECT(staten = allocVector(VECSXP, nparam + 1));
550
        PROTECT(staten = allocVector(VECSXP, nparam + 1));
558
        PROTECT(namesn = allocVector(STRSXP, nparam + 1));
551
        PROTECT(namesn = allocVector(STRSXP, nparam + 1));
559
        for (int j = 0; j < nparam; ++j) {
552
        for (int j = 0; j < nparam; ++j) {
560
        SET_ELEMENT(staten, j, VECTOR_ELT(params, j));
553
        SET_ELEMENT(staten, j, VECTOR_ELT(params, j));
561
        SET_STRING_ELT(namesn, j, STRING_ELT(names, j));
554
        SET_STRING_ELT(namesn, j, STRING_ELT(names, j));
562
        }
555
        }
563
        //Assign .RNG.name as the last element
556
        //Assign .RNG.name as the last element
564
        SEXP rngname;
557
        SEXP rngname;
565
      PROTECT(rngname = allocVector(STRSXP,1));
558
      PROTECT(rngname = mkString(srng.c_str()));
566
      SET_STRING_ELT(rngname, 0, mkChar(srng.c_str()));
567
        SET_ELEMENT(staten, nparam, rngname);
559
        SET_ELEMENT(staten, nparam, rngname);
568
        SET_STRING_ELT(namesn, nparam, mkChar(".RNG.name"));
560
        SET_STRING_ELT(namesn, nparam, mkChar(".RNG.name"));
569
        setAttrib(staten, R_NamesSymbol, namesn);
561
        setAttrib(staten, R_NamesSymbol, namesn);
570
        //And we're done with this chain
562
        //And we're done with this chain
571
        SET_ELEMENT(ans, n, staten);
563
        SET_ELEMENT(ans, n, staten);
572
        UNPROTECT(5); //rngname, namesn, statesn, names, params
564
        UNPROTECT(4); //rngname, namesn, staten, params
573
    }
565
    }
574
    UNPROTECT(1); //ans
566
    UNPROTECT(1); //ans
575
    return ans;
567
    return ans;
576
    }
568
    }
577
569
...
...
583
    SEXP varnames;
575
    SEXP varnames;
584
    PROTECT(varnames = allocVector(STRSXP,namevec.size()));
576
    PROTECT(varnames = allocVector(STRSXP,namevec.size()));
585
    for (unsigned int i = 0; i < namevec.size(); ++i) {
577
    for (unsigned int i = 0; i < namevec.size(); ++i) {
586
        SET_STRING_ELT(varnames, i, mkChar(namevec[i].c_str()));
578
        SET_STRING_ELT(varnames, i, mkChar(namevec[i].c_str()));
587
    }
579
    }
588
  UNPROTECT(1);
580
  UNPROTECT(1); //varnames
589
    return varnames;
581
    return varnames;
590
    }
582
    }
591
583
592
    SEXP get_samplers(SEXP ptr)
584
    SEXP get_samplers(SEXP ptr)
593
    {
585
    {
...
...
611
        SET_ELEMENT(node_list, i, e);
603
        SET_ELEMENT(node_list, i, e);
612
        SET_STRING_ELT(sampler_names, i, mkChar(samplers[i][0].c_str()));
604
        SET_STRING_ELT(sampler_names, i, mkChar(samplers[i][0].c_str()));
613
        UNPROTECT(1); //e
605
        UNPROTECT(1); //e
614
    }
606
    }
615
    setAttrib(node_list, R_NamesSymbol, sampler_names); 
607
    setAttrib(node_list, R_NamesSymbol, sampler_names); 
616
  UNPROTECT(2); //names, ans
608
  UNPROTECT(2); //node_list, sampler_names
617
    return node_list;
609
    return node_list;
618
    }
610
    }
619
611
620
    SEXP get_factories(SEXP type)
612
    SEXP get_factories(SEXP type)
621
    {
613
    {
622
    FactoryType ft = asFactoryType(type);
614
    FactoryType ft = asFactoryType(type);
623
    vector<pair<string, bool> > factories = Console::listFactories(ft);
615
    vector<pair<string, bool> > factories = Console::listFactories(ft);
624
      
625
    unsigned int n = factories.size();
616
    unsigned int n = factories.size();
617
618
  SEXP fac_list;
619
  PROTECT(fac_list = allocVector(VECSXP, 2));
620
626
    SEXP names, status;
621
    SEXP names, status;
627
    PROTECT(names = allocVector(STRSXP, n));
622
    PROTECT(names = allocVector(STRSXP, n));
628
    PROTECT(status = allocVector(LGLSXP, n));
623
    PROTECT(status = allocVector(LGLSXP, n));
629
    for (unsigned int i = 0; i < n; ++i) {
624
    for (unsigned int i = 0; i < n; ++i) {
630
        SET_STRING_ELT(names, i, mkChar(factories[i].first.c_str()));
625
        SET_STRING_ELT(names, i, mkChar(factories[i].first.c_str()));
631
        LOGICAL_POINTER(status)[i] = factories[i].second;
626
        LOGICAL_POINTER(status)[i] = factories[i].second;
632
    }
627
    }
633
628
634
  SEXP fac_list;
635
  PROTECT(fac_list = allocVector(VECSXP, 2));
636
    SET_ELEMENT(fac_list, 0, names);
629
    SET_ELEMENT(fac_list, 0, names);
637
    SET_ELEMENT(fac_list, 1, status);
630
    SET_ELEMENT(fac_list, 1, status);
638
    UNPROTECT(2); //names, status
631
    UNPROTECT(2); //names, status
639
632
640
    SEXP fac_names;
633
    SEXP fac_names;
...
...
655
    return R_NilValue;
648
    return R_NilValue;
656
    }
649
    }
657
650
658
    SEXP get_iter(SEXP ptr)
651
    SEXP get_iter(SEXP ptr)
659
    {
652
    {
660
  Console *console = ptrArg(ptr);
653
  return ScalarInteger(ptrArg(ptr)->iter());
661
  unsigned int iter = console->iter();
662
663
  SEXP ans;
664
  PROTECT(ans = allocVector(INTSXP, 1));
665
  INTEGER(ans)[0] = iter;
666
  UNPROTECT(1);
667
  return ans;
668
    }
654
    }
669
    
655
    
670
    SEXP get_nchain(SEXP ptr)
656
    SEXP get_nchain(SEXP ptr)
671
    {
657
    {
672
  Console *console = ptrArg(ptr);
658
  return ScalarInteger(ptrArg(ptr)->nchain());
673
  unsigned int nchain = console->nchain();
674
    
675
  SEXP ans;
676
  PROTECT(ans = allocVector(INTSXP,1));
677
  INTEGER(ans)[0] = nchain;
678
  UNPROTECT(1);
679
  return ans;
680
    }
659
    }
681
660
682
    SEXP load_module(SEXP name)
661
    SEXP load_module(SEXP name)
683
    {
662
    {
684
    return ScalarLogical(Console::loadModule(stringArg(name)));
663
    return ScalarLogical(Console::loadModule(stringArg(name)));