#include <config.h> #include <model/NodeArray.h> #include <graph/ConstantNode.h> #include <graph/StochasticNode.h> #include <graph/AggNode.h> #include <sarray/RangeIterator.h> #include <graph/NodeError.h> #include <sarray/SArray.h> #include <util/nainf.h> #include <string> #include <stdexcept> #include <limits> using std::vector; using std::map; using std::string; using std::runtime_error; using std::logic_error; using std::set; using std::numeric_limits; static ConstantNode const *asConstant(Node const *node) { return dynamic_cast<ConstantNode const *>(node); } NodeArray::NodeArray(string const &name, vector<unsigned int> const &dim, unsigned int nchain) : _name(name), _range(dim), _nchain(nchain) { unsigned int length = _range.length(); _node_pointers = new Node *[length]; _offsets = new unsigned int[length]; for (unsigned int i = 0; i < length; i++) { _node_pointers[i] = 0; _offsets[i] = numeric_limits<unsigned int>::max(); } } NodeArray::~NodeArray() { delete [] _node_pointers; delete [] _offsets; } 00046 bool NodeArray::isEmpty(Range const &target_range) const { if (!_range.contains(target_range)) throw logic_error("Range error in NodeArray::isEmpty"); for (RangeIterator i(target_range); !i.atEnd(); i.nextLeft()) { if (_node_pointers[_range.leftOffset(i)] != 0) return false; } return true; } 00058 void NodeArray::insert(Node *node, Range const &target_range) { if (!node) { throw logic_error(string("Attempt to insert NULL node at ") + name() + print(target_range)); } if (node->dim() != target_range.dim(true)) { throw runtime_error(string("Cannot insert node into ") + name() + print(target_range) + ". Dimension mismatch"); } if (!_range.contains(target_range)) { throw runtime_error(string("Cannot insert node into ") + name() + print(target_range) + ". Range out of bounds"); } if (!isEmpty(target_range)) { throw runtime_error(string("Node ") + name() + print(target_range) + " overlaps previously defined nodes"); } /* Set the _node_pointers array and the offset array */ RangeIterator j(target_range); for (unsigned int k = 0; !j.atEnd(); j.nextLeft(), ++k) { unsigned int offset = _range.leftOffset(j); _node_pointers[offset] = node; _offsets[offset] = k; } /* Add to the graph */ _member_graph.add(node); } 00090 Node *NodeArray::find(Range const &target_range) const { if (!_range.contains(target_range)) { return 0; } unsigned int offset = _range.leftOffset(target_range.lower()); Node *node = _node_pointers[offset]; if (!node) return 0; if (node->dim() != target_range.dim(true)) return 0; RangeIterator j(target_range); for (unsigned int k = 0; !j.atEnd(); j.nextLeft(), ++k) { offset = _range.leftOffset(j); if (_node_pointers[offset] != node || _offsets[offset] != k) return 0; } return node; } 00114 Node *NodeArray::getSubset(Range const &target_range, Graph &graph) { //Check validity of target_range if (!_range.contains(target_range)) { throw runtime_error(string("Cannot get subset ") + _name + print(target_range) + ". Range out of bounds"); } /* If range corresponds to a set node, then return this */ Node *node = find(target_range); if (node) return node; /* If range corresponds to a previously created subset, then return this */ map<Range, Node *>::iterator p = _generated_nodes.find(target_range); if (p != _generated_nodes.end()) { return p->second; } /* Otherwise create an aggregate node */ vector<Node const *> nodes; vector<unsigned int> offsets; for (RangeIterator i(target_range); !i.atEnd(); i.nextLeft()) { unsigned int offset = _range.leftOffset(i); if (_node_pointers[offset] == 0) { return 0; } nodes.push_back(_node_pointers[offset]); offsets.push_back(_offsets[offset]); } node = new AggNode(target_range.dim(true), nodes, offsets); _generated_nodes.insert(std::pair<Range,Node*>(target_range, node)); graph.add(node); _member_graph.add(node); return node; } 00150 void NodeArray::setValue(SArray const &value, unsigned int chain) { if (!(_range == value.range())) { throw runtime_error(string("Dimension mismatch when setting value of node array ") + name()); } vector<double> const &x = value.value(); unsigned int N = value.length(); //Gather all the nodes for which a data value is supplied set<Node*> setnodes; for (unsigned int i = 0; i < _range.length(); ++i) { if (x[i] != JAGS_NA) { Node *node = _node_pointers[i]; if (node == 0) { string msg = "Attempt to set value of undefined node"; throw runtime_error(msg + name() + print(value.range().leftIndex(i))); } if (asStochastic(node) || asConstant(node)) { setnodes.insert(node); } else { throw NodeError(node, "Attempt to set value of non-variable node"); } } } set<Node*>::const_iterator p; double *node_value = new double[N]; for (p = setnodes.begin(); p != setnodes.end(); ++p) { //Step through each node Node *node = *p; //Get vector of values for this node for (unsigned int i = 0; i < N; ++i) { if (_node_pointers[i] == node) { if (_offsets[i] > node->length()) { throw logic_error("Invalid offset in NodeArray::setValue"); } else { node_value[_offsets[i]] = x[i]; } } } // If there are any missing values, they must all be missing bool missing = node_value[0] == JAGS_NA; for (unsigned int j = 1; j < node->length(); ++j) { if ((node_value[j] == JAGS_NA) != missing) { delete [] node_value; throw NodeError(node,"Values supplied for node are partially missing"); } } if (!missing) { node->setValue(node_value, node->length(), chain); } } delete [] node_value; } 00211 void NodeArray::getValue(SArray &value, unsigned int chain, bool (*condition)(Node const *)) const { if (!(_range == value.range())) { string msg("Dimension mismatch when getting value of node array "); msg.append(name()); throw runtime_error(msg); } unsigned int array_length = _range.length(); vector<double> array_value(array_length); for (unsigned int j = 0; j < array_length; ++j) { Node const *node = _node_pointers[j]; if (node && condition(node)) { array_value[j] = node->value(chain)[_offsets[j]]; } else { array_value[j] = JAGS_NA; } } value.setValue(array_value); } //FIXME: A lot of code overlap with setValue here. #include <iostream> 00238 void NodeArray::setData(SArray const &value, Graph &graph) { if (!(_range == value.range())) { throw runtime_error(string("Dimension mismatch when setting value of node array ") + name()); } unsigned int N = value.length(); vector<double> const &x = value.value(); //Gather all the nodes for which a data value is supplied set<Node*> setnodes; for (unsigned int i = 0; i < _range.length(); ++i) { if (x[i] != JAGS_NA) { if (_node_pointers[i] == 0) { //Insert a new constant node ConstantNode *cnode = new ConstantNode(x[i], _nchain); graph.add(cnode); insert(cnode, _range.leftIndex(i)); } else { //Existing node for which we must set value setnodes.insert(_node_pointers[i]); } } } set<Node*>::const_iterator p; for (p = setnodes.begin(); p != setnodes.end(); ++p) { //Step through each node Node *node = *p; vector<double> node_value(node->length()); //Get vector of values for this node for (unsigned int i = 0; i < N; ++i) { if (_node_pointers[i] == node) { if (_offsets[i] > node->length()) { throw logic_error("Invalid offset in NodeArray::setValue"); } else { node_value[_offsets[i]] = x[i]; } } } // If there are any missing values, they must all be missing bool missing = (node_value[0] == JAGS_NA); for (unsigned int j = 1; j < node->length(); ++j) { if ((node_value[j] == JAGS_NA) != missing) { throw NodeError(node,"Values supplied for node are partially missing"); } } if (!missing) { node->setObserved(node_value); } } } 00294 string const &NodeArray::name() const { return _name; } 00299 Range const &NodeArray::range() const { return _range; } /* Graph const &NodeArray::graph() const { return _graph; } */ bool NodeArray::findActiveIndices(vector<unsigned int> &ind, unsigned int k, vector<int> const &lower, vector<unsigned int> const &dim) const { /* We pay a heavy computational price for the flexibility of allowing users to insert multivariate nodes in arbritary ways into the NodeArray. Suppose we have an array of dimension [3,4,2,5], the lower index is [1,2,1,2] and the dimension of the node is [3,2]. Then the node could be inserted in 5 different ways. [1:3, 2:3, 1, 2] Active indices (0,1) [1:3, 2, 1:2, 2] (0,2) [1:3, 2, 1, 2:3] (0,3) [1, 2:4, 1:2, 2] (1,2) [1, 2:4, 1, 2:3] (1,3) We can't have active indices (2,3) because the node won't fit */ if (k == 0) ind[k] = 0; else ind[k] = ind[k-1] + 1; unsigned int m = ind.size(); unsigned int M = _range.ndim(false); for (;ind[k] + m <= M + k; ind[k] = ind[k] + 1) { if (k == m - 1) { vector<int> upper(lower); for (unsigned int l = 0; l < m; ++l) { upper[ind[l]] = upper[ind[l]] + dim[l] - 1; } Range test_range(lower, upper); if (_range.contains(test_range)) { Node *node = _node_pointers[_range.leftOffset(lower)]; unsigned int j = 0; bool ok = true; for (RangeIterator i(test_range); !i.atEnd(); i.nextLeft(), ++j) { unsigned int offset = _range.leftOffset(i); if (_node_pointers[offset] != node || _offsets[offset] != j) { ok = false; break; } } if (ok) return true; } } else { if (findActiveIndices(ind, k+1, lower, dim)) return true; } } return false; } 00368 Range NodeArray::getRange(Node const *node) const { if (!_member_graph.contains(node)) { return Range(); } //Look in the generated nodes first for (map<Range, Node *>::const_iterator p = _generated_nodes.begin(); p != _generated_nodes.end(); ++p) { if (node == p->second) return p->first; } /* Find the lower limit of the range. This is easy */ unsigned int ndim = _range.ndim(false); vector<int> lower(ndim); unsigned int j = 0; for (; j < _range.length(); ++j) { if (_node_pointers[j] == node) { lower = _range.leftIndex(j); break; } } if (j == _range.length()) { return Range(); } unsigned int m = node->dim().size(); vector<unsigned int> ind(m, 1); if (findActiveIndices(ind, 0, lower, node->dim())) { vector<int> upper = lower; for (unsigned int l = 0; l < m; ++l) { upper[ind[l]] = upper[ind[l]] + node->dim()[l] - 1; } return Range(lower, upper); } else { throw logic_error("Unable to find node range"); } } 00410 unsigned int NodeArray::nchain() const { return _nchain; }