diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java index cffbaad71f9c59b41e3cae935605ccf4bf441156..1d88093b56227120073c643e15d546084b029479 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java @@ -12,9 +12,9 @@ */ package com.oracle.truffle.r.library.stats; -import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.nullValue; -import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.missingValue; import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.instanceOf; +import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.missingValue; +import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.nullValue; import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN; import com.oracle.truffle.api.dsl.Cached; @@ -25,8 +25,6 @@ import com.oracle.truffle.r.nodes.attributes.SetAttributeNode; import com.oracle.truffle.r.nodes.attributes.SpecialAttributesFunctions.GetDimAttributeNode; import com.oracle.truffle.r.nodes.attributes.SpecialAttributesFunctions.SetClassAttributeNode; import com.oracle.truffle.r.nodes.builtin.RExternalBuiltinNode; -import com.oracle.truffle.r.runtime.data.nodes.ReadAccessor; -import com.oracle.truffle.r.runtime.data.nodes.VectorReadAccess; import com.oracle.truffle.r.runtime.RError; import com.oracle.truffle.r.runtime.RRuntime; import com.oracle.truffle.r.runtime.data.RDataFactory; @@ -34,9 +32,12 @@ import com.oracle.truffle.r.runtime.data.RDoubleVector; import com.oracle.truffle.r.runtime.data.RList; import com.oracle.truffle.r.runtime.data.RStringVector; import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess.RandomIterator; import com.oracle.truffle.r.runtime.ops.na.NACheck; public abstract class Cdist extends RExternalBuiltinNode.Arg4 { + private static final NACheck naCheck = NACheck.create(); @Child private GetFixedAttributeNode getNamesAttrNode = GetFixedAttributeNode.createNames(); @@ -49,9 +50,10 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { casts.arg(3).asDoubleVector().findFirst(); } - @Specialization(guards = "method == cachedMethod") - protected RDoubleVector cdist(RAbstractDoubleVector x, @SuppressWarnings("unused") int method, RList list, double p, @SuppressWarnings("unused") @Cached("method") int cachedMethod, - @Cached("create()") VectorReadAccess.Double xAccess, + @Specialization(guards = {"method == cachedMethod", "xAccess.supports(x)"}) + protected RDoubleVector cdist(RAbstractDoubleVector x, @SuppressWarnings("unused") int method, RList list, double p, + @Cached("method") @SuppressWarnings("unused") int cachedMethod, + @Cached("x.access()") VectorAccess xAccess, @Cached("getMethod(method)") Method methodObj, @Cached("create()") SetAttributeNode setAttrNode, @Cached("create()") SetClassAttributeNode setClassAttrNode, @@ -60,8 +62,10 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { int nc = getDimNode.ncols(x); int n = nr * (nr - 1) / 2; /* avoid int overflow for N ~ 50,000 */ double[] ans = new double[n]; - RDoubleVector xm = x.materialize(); - rdistance(new ReadAccessor.Double(x, xAccess), nr, nc, ans, false, methodObj, p); + + try (RandomIterator xIter = xAccess.randomAccess(x)) { + rdistance(xAccess, xIter, nr, nc, ans, false, methodObj, p); + } RDoubleVector result = RDataFactory.createDoubleVector(ans, naCheck.neverSeenNA()); DynamicObject resultAttrs = result.initAttributes(); @@ -81,6 +85,14 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { return result; } + @Specialization(replaces = "cdist") + protected RDoubleVector cdistGeneric(RAbstractDoubleVector x, int method, RList list, double p, + @Cached("create()") SetAttributeNode setAttrNode, + @Cached("create()") SetClassAttributeNode setClassAttrNode, + @Cached("create()") GetDimAttributeNode getDimNode) { + return cdist(x, method, list, p, method, x.slowPathAccess(), getMethod(method), setAttrNode, setClassAttrNode, getDimNode); + } + private static boolean bothNonNAN(double a, double b) { return !RRuntime.isNAorNaN(a) && !RRuntime.isNAorNaN(b); } @@ -96,7 +108,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { return Method.values()[method - 1]; } - private void rdistance(ReadAccessor.Double xAccess, int nr, int nc, double[] d, boolean diag, Method method, double p) { + private void rdistance(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, double[] d, boolean diag, Method method, double p) { int ij; /* can exceed 2^31 - 1, but Java can't handle that */ // if (method == Method.MINKOWSKI) { @@ -109,7 +121,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { naCheck.enable(true); for (int j = 0; j <= nr; j++) { for (int i = j + dc; i < nr; i++) { - double r = method.dist(xAccess, nr, nc, i, j, p); + double r = method.dist(xAccess, xIter, nr, nc, i, j, p); naCheck.check(r); d[ij++] = r; } @@ -119,7 +131,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { public enum Method { EUCLIDEAN { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; double dev; @@ -130,8 +142,8 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { count = 0; dist = 0; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - dev = (xAccess.getDataAt(i1) - xAccess.getDataAt(i2)); + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + dev = (xAccess.getDouble(xIter, i1) - xAccess.getDouble(xIter, i2)); if (!RRuntime.isNAorNaN(dev)) { dist += dev * dev; count++; @@ -152,7 +164,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { }, MAXIMUM { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; double dev; @@ -163,8 +175,8 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { count = 0; dist = -Double.MAX_VALUE; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - dev = Math.abs(xAccess.getDataAt(i1) - xAccess.getDataAt(i2)); + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + dev = Math.abs(xAccess.getDouble(xIter, i1) - xAccess.getDouble(xIter, i2)); if (!RRuntime.isNAorNaN(dev)) { if (dev > dist) { dist = dev; @@ -184,7 +196,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { }, MANHATTAN { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; double dev; @@ -195,8 +207,8 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { count = 0; dist = 0; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - dev = Math.abs(xAccess.getDataAt(i1) - xAccess.getDataAt(i2)); + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + dev = Math.abs(xAccess.getDouble(xIter, i1) - xAccess.getDouble(xIter, i2)); if (!RRuntime.isNAorNaN(dev)) { dist += dev; count++; @@ -217,7 +229,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { }, CANBERRA { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; double dev; @@ -230,9 +242,9 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { count = 0; dist = 0; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - sum = Math.abs(xAccess.getDataAt(i1) + xAccess.getDataAt(i2)); - diff = Math.abs(xAccess.getDataAt(i1) - xAccess.getDataAt(i2)); + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + sum = Math.abs(xAccess.getDouble(xIter, i1) + xAccess.getDouble(xIter, i2)); + diff = Math.abs(xAccess.getDouble(xIter, i1) - xAccess.getDouble(xIter, i2)); if (sum > DBL_MIN || diff > DBL_MIN) { dev = diff / sum; if (!RRuntime.isNAorNaN(dev) || @@ -258,7 +270,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { }, BINARY { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; int total; @@ -271,13 +283,13 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { dist = 0; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - if (!bothFinite(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + if (!bothFinite(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { RError.warning(RError.SHOW_CALLER2, RError.Message.GENERIC, "treating non-finite values as NA"); } else { - if (xAccess.getDataAt(i1) != 0. || xAccess.getDataAt(i2) != 0.) { + if (xAccess.getDouble(xIter, i1) != 0. || xAccess.getDouble(xIter, i2) != 0.) { count++; - if (!(xAccess.getDataAt(i1) != 0. && xAccess.getDataAt(i2) != 0.)) { + if (!(xAccess.getDouble(xIter, i1) != 0. && xAccess.getDouble(xIter, i2) != 0.)) { dist++; } } @@ -300,7 +312,7 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { }, MINKOWSKI { @Override - public double dist(ReadAccessor.Double xAccess, int nr, int nc, final int i1in, final int i2in, double p) { + public double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, final int i1in, final int i2in, double p) { int i1 = i1in; int i2 = i2in; double dev; @@ -311,8 +323,8 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { count = 0; dist = 0; for (j = 0; j < nc; j++) { - if (bothNonNAN(xAccess.getDataAt(i1), xAccess.getDataAt(i2))) { - dev = (xAccess.getDataAt(i1) - xAccess.getDataAt(i2)); + if (bothNonNAN(xAccess.getDouble(xIter, i1), xAccess.getDouble(xIter, i2))) { + dev = (xAccess.getDouble(xIter, i1) - xAccess.getDouble(xIter, i2)); if (!RRuntime.isNAorNaN(dev)) { dist += Math.pow(Math.abs(dev), p); count++; @@ -331,6 +343,6 @@ public abstract class Cdist extends RExternalBuiltinNode.Arg4 { } }; - public abstract double dist(ReadAccessor.Double xAccess, int nr, int nc, int i1, int i2, double p); + public abstract double dist(VectorAccess xAccess, RandomIterator xIter, int nr, int nc, int i1, int i2, double p); } } diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cutree.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cutree.java index 0a07a8b78c67ee7965a500d2cb61c1561e9feffc..8c0052139f253fd4233b127aeefe7f0cedaeb1df 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cutree.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cutree.java @@ -20,7 +20,8 @@ import com.oracle.truffle.r.nodes.builtin.RExternalBuiltinNode; import com.oracle.truffle.r.runtime.data.RDataFactory; import com.oracle.truffle.r.runtime.data.RIntVector; import com.oracle.truffle.r.runtime.data.model.RAbstractIntVector; -import com.oracle.truffle.r.runtime.data.nodes.VectorReadAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess.RandomIterator; // translated from library/stats/src/hclust_utils.c @@ -32,10 +33,10 @@ public abstract class Cutree extends RExternalBuiltinNode.Arg2 { casts.arg(1).mustNotBeMissing().mapIf(nullValue(), emptyIntegerVector()).asIntegerVector(); } - @Specialization + @Specialization(guards = {"mergeAccess.supports(merge)", "whichAccess.supports(which)"}) protected RIntVector cutree(RAbstractIntVector merge, RAbstractIntVector which, - @Cached("create()") VectorReadAccess.Int mergeAccess, - @Cached("create()") VectorReadAccess.Int whichAccess, + @Cached("merge.access()") VectorAccess mergeAccess, + @Cached("which.access()") VectorAccess whichAccess, @Cached("create()") GetDimAttributeNode getDimNode) { int whichLen = which.getLength(); @@ -59,92 +60,96 @@ public abstract class Cutree extends RExternalBuiltinNode.Arg2 { int[] z = new int[n]; int[] iAns = new int[n * whichLen]; - Object mergeStore = mergeAccess.getDataStore(merge); - Object whichStore = whichAccess.getDataStore(which); + try (RandomIterator mergeIter = mergeAccess.randomAccess(merge); RandomIterator whichIter = whichAccess.randomAccess(which)) { - // for (k = 1; k <= n; k++) { - for (k = 0; k < n; k++) { - sing[k] = true; /* is k-th obs. still alone in cluster ? */ - mNr[k] = 0; /* containing last merge-step number of k-th obs. */ - } + // for (k = 1; k <= n; k++) { + for (k = 0; k < n; k++) { + sing[k] = true; /* is k-th obs. still alone in cluster ? */ + mNr[k] = 0; /* containing last merge-step number of k-th obs. */ + } - for (k = 1; k <= n - 1; k++) { - /* k-th merge, from n-k+1 to n-k atoms: (m1,m2) = merge[ k , ] */ - m1 = mergeAccess.getDataAt(merge, mergeStore, k - 1); - m2 = mergeAccess.getDataAt(merge, mergeStore, n - 1 + k - 1); - - if (m1 < 0 && m2 < 0) { /* merging atoms [-m1] and [-m2] */ - mNr[adj(-m1)] = mNr[adj(-m2)] = k; - sing[adj(-m1)] = sing[adj(-m2)] = false; - } else if (m1 < 0 || m2 < 0) { /* the other >= 0 */ - if (m1 < 0) { - j = -m1; - m1 = m2; - } else { - j = -m2; - } - /* merging atom j & cluster m1 */ - for (l = 1; l <= n; l++) { - if (mNr[adj(l)] == m1) { - mNr[adj(l)] = k; + for (k = 1; k <= n - 1; k++) { + /* k-th merge, from n-k+1 to n-k atoms: (m1,m2) = merge[ k , ] */ + m1 = mergeAccess.getInt(mergeIter, k - 1); + m2 = mergeAccess.getInt(mergeIter, n - 1 + k - 1); + + if (m1 < 0 && m2 < 0) { /* merging atoms [-m1] and [-m2] */ + mNr[adj(-m1)] = mNr[adj(-m2)] = k; + sing[adj(-m1)] = sing[adj(-m2)] = false; + } else if (m1 < 0 || m2 < 0) { /* the other >= 0 */ + if (m1 < 0) { + j = -m1; + m1 = m2; + } else { + j = -m2; } - } - mNr[adj(j)] = k; - sing[adj(j)] = false; - } else { /* both m1, m2 >= 0 */ - for (l = 1; l <= n; l++) { - if (mNr[adj(l)] == m1 || mNr[adj(l)] == m2) { - mNr[adj(l)] = k; + /* merging atom j & cluster m1 */ + for (l = 1; l <= n; l++) { + if (mNr[adj(l)] == m1) { + mNr[adj(l)] = k; + } + } + mNr[adj(j)] = k; + sing[adj(j)] = false; + } else { /* both m1, m2 >= 0 */ + for (l = 1; l <= n; l++) { + if (mNr[adj(l)] == m1 || mNr[adj(l)] == m2) { + mNr[adj(l)] = k; + } } } - } - /* - * does this k-th merge belong to a desired group size which[j] ? if yes, find j (maybe - * multiple ones): - */ - foundJ = false; - for (j = 0; j < whichLen; j++) { - if (whichAccess.getDataAt(which, whichStore, j) == n - k) { - if (!foundJ) { /* first match (and usually only one) */ - foundJ = true; - // for (l = 1; l <= n; l++) - for (l = 0; l < n; l++) { - z[l] = 0; - } - nclust = 0; - mm = j * n; /* may want to copy this column of ans[] */ - for (l = 1, m1 = mm; l <= n; l++, m1++) { - if (sing[adj(l)]) { - iAns[m1] = ++nclust; - } else { - if (z[adj(mNr[adj(l)])] == 0) { - z[adj(mNr[adj(l)])] = ++nclust; + /* + * does this k-th merge belong to a desired group size which[j] ? if yes, find j + * (maybe multiple ones): + */ + foundJ = false; + for (j = 0; j < whichLen; j++) { + if (whichAccess.getInt(whichIter, j) == n - k) { + if (!foundJ) { /* first match (and usually only one) */ + foundJ = true; + // for (l = 1; l <= n; l++) + for (l = 0; l < n; l++) { + z[l] = 0; + } + nclust = 0; + mm = j * n; /* may want to copy this column of ans[] */ + for (l = 1, m1 = mm; l <= n; l++, m1++) { + if (sing[adj(l)]) { + iAns[m1] = ++nclust; + } else { + if (z[adj(mNr[adj(l)])] == 0) { + z[adj(mNr[adj(l)])] = ++nclust; + } + iAns[m1] = z[adj(mNr[adj(l)])]; } - iAns[m1] = z[adj(mNr[adj(l)])]; + } + } else { /* found_j: another which[j] == n-k : copy column */ + for (l = 1, m1 = j * n, m2 = mm; l <= n; l++, m1++, m2++) { + iAns[m1] = iAns[m2]; } } - } else { /* found_j: another which[j] == n-k : copy column */ - for (l = 1, m1 = j * n, m2 = mm; l <= n; l++, m1++, m2++) { - iAns[m1] = iAns[m2]; - } + } /* if ( match ) */ + } /* for(j .. which[j] ) */ + } /* for(k ..) {merge} */ + + /* Dealing with trivial case which[] = n : */ + for (j = 0; j < whichLen; j++) { + if (whichAccess.getInt(whichIter, j) == n) { + for (l = 1, m1 = j * n; l <= n; l++, m1++) { + iAns[m1] = l; } - } /* if ( match ) */ - } /* for(j .. which[j] ) */ - } /* for(k ..) {merge} */ - - /* Dealing with trivial case which[] = n : */ - for (j = 0; j < whichLen; j++) { - if (whichAccess.getDataAt(which, whichStore, j) == n) { - for (l = 1, m1 = j * n; l <= n; l++, m1++) { - iAns[m1] = l; } } } - RIntVector result = RDataFactory.createIntVector(iAns, RDataFactory.COMPLETE_VECTOR, new int[]{n, whichLen}); - return result; + return RDataFactory.createIntVector(iAns, RDataFactory.COMPLETE_VECTOR, new int[]{n, whichLen}); + } + @Specialization(replaces = "cutree") + protected RIntVector cutreeGeneric(RAbstractIntVector merge, RAbstractIntVector which, + @Cached("create()") GetDimAttributeNode getDimNode) { + return cutree(merge, which, merge.slowPathAccess(), which.slowPathAccess(), getDimNode); } private static int adj(int i) { diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DoubleCentre.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DoubleCentre.java index 62883376e19d6cccdc6e3ab014c7765086f15e4e..a817389e28340b1a2cca0112feefc0260960de6f 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DoubleCentre.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DoubleCentre.java @@ -10,56 +10,63 @@ */ package com.oracle.truffle.r.library.stats; -import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.missingValue; -import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.nullValue; - import com.oracle.truffle.api.dsl.Cached; import com.oracle.truffle.api.dsl.Specialization; import com.oracle.truffle.r.nodes.attributes.SpecialAttributesFunctions.GetDimAttributeNode; import com.oracle.truffle.r.nodes.builtin.RExternalBuiltinNode; import com.oracle.truffle.r.runtime.RError; -import com.oracle.truffle.r.runtime.data.RDoubleVector; import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; -import com.oracle.truffle.r.runtime.data.nodes.SetDataAt; -import com.oracle.truffle.r.runtime.data.nodes.VectorReadAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess.RandomIterator; +import com.oracle.truffle.r.runtime.data.nodes.VectorReuse; public abstract class DoubleCentre extends RExternalBuiltinNode.Arg1 { static { Casts casts = new Casts(DoubleCentre.class); - casts.arg(0).mustBe(missingValue().not()).mustBe(nullValue().not(), RError.Message.MACRO_CAN_BE_APPLIED_TO, "REAL()", "numeric", "NULL").asDoubleVector(); + casts.arg(0).mustNotBeNull(RError.Message.MACRO_CAN_BE_APPLIED_TO, "REAL()", "numeric", "NULL").mustNotBeMissing().asDoubleVector(); } - @Specialization - protected RDoubleVector doubleCentre(RAbstractDoubleVector aVecAbs, - @Cached("create()") VectorReadAccess.Double aAccess, - @Cached("create()") SetDataAt.Double aSetter, + @Specialization(guards = {"aAccess.supports(a)", "reuse.supports(a)"}) + protected RAbstractDoubleVector doubleCentre(RAbstractDoubleVector a, + @Cached("a.access()") VectorAccess aAccess, + @Cached("createNonShared(a)") VectorReuse reuse, @Cached("create()") GetDimAttributeNode getDimNode) { - RDoubleVector aVec = aVecAbs.materialize(); - int n = getDimNode.nrows(aVec); - Object aStore = aAccess.getDataStore(aVec); - for (int i = 0; i < n; i++) { - double sum = 0; - for (int j = 0; j < n; j++) { - sum += aAccess.getDataAt(aVec, aStore, i + j * n); - } - sum /= n; - for (int j = 0; j < n; j++) { - double val = aAccess.getDataAt(aVec, aStore, i + j * n); - aSetter.setDataAt(aVec, aStore, i + j * n, val - sum); - } - } - for (int j = 0; j < n; j++) { - double sum = 0; - for (int i = 0; i < n; i++) { - sum += aAccess.getDataAt(aVec, aStore, i + j * n); - } - sum /= n; - for (int i = 0; i < n; i++) { - double val = aAccess.getDataAt(aVec, aStore, i + j * n); - aSetter.setDataAt(aVec, aStore, i + j * n, val - sum); + int n = getDimNode.nrows(a); + + try (RandomIterator aIter = aAccess.randomAccess(a)) { + RAbstractDoubleVector result = reuse.getResult(a); + VectorAccess resultAccess = reuse.access(result); + try (RandomIterator resultIter = resultAccess.randomAccess(result)) { + for (int i = 0; i < n; i++) { + double sum = 0; + for (int j = 0; j < n; j++) { + sum += aAccess.getDouble(aIter, i + j * n); + } + sum /= n; + for (int j = 0; j < n; j++) { + resultAccess.setDouble(resultIter, i + j * n, aAccess.getDouble(aIter, i + j * n) - sum); + } + } + for (int j = 0; j < n; j++) { + double sum = 0; + for (int i = 0; i < n; i++) { + sum += resultAccess.getDouble(aIter, i + j * n); + } + sum /= n; + for (int i = 0; i < n; i++) { + resultAccess.setDouble(resultIter, i + j * n, resultAccess.getDouble(aIter, i + j * n) - sum); + } + } } + return result; } - return aVec; + } + + @Specialization(replaces = "doubleCentre") + protected RAbstractDoubleVector doubleCentreGeneric(RAbstractDoubleVector a, + @Cached("createNonSharedGeneric()") VectorReuse reuse, + @Cached("create()") GetDimAttributeNode getDimNode) { + return doubleCentre(a, a.slowPathAccess(), reuse, getDimNode); } }