both algorithms combined for awesome performance
481487
function simplify(points, tolerance, highestQuality) {

            if (points.length <= 1)
                return points;

            var sqTolerance = tolerance !== undefined ? tolerance * tolerance
                    : 1;

            points = highestQuality ? points : simplifyRadialDist(points,
                    sqTolerance);
            points = simplifyDouglasPeucker(points, sqTolerance);

            return points;
        }
362837
public void performStep() {
        // check for zeros
        for( int i = helper.x2-1; i >= helper.x1; i-- ) {
            if( helper.isZero(i) ) {
                helper.splits[helper.numSplits++] = i;
                helper.x1 = i+1;
                return;
            }
        }

        double lambda;

        if( followingScript ) {
            if( helper.steps > 10 ) {
                followingScript = false;
                return;
            } else {
                // Using the true eigenvalues will in general lead to the fastest convergence
                // typically takes 1 or 2 steps
                lambda = eigenvalues[helper.x2];
            }
        } else {
            // the current eigenvalue isn't working so try something else
            lambda = helper.computeShift();
        }

        // similar transforms
        helper.performImplicitSingleStep(lambda,false);
    }
953917
public void computeOptimalStrategy() {
		long heavyMin, revHeavyMin, leftMin, revLeftMin, rightMin, revRightMin;
		long min = -1;
		int strategy = -1;
		int parent1 = -1;
		int parent2 = -1;
		boolean[] nodeTypeLeft1 = it1.nodeType[LEFT];
		boolean[] nodeTypeLeft2 = it2.nodeType[LEFT];
		boolean[] nodeTypeRigt1 = it1.nodeType[RIGHT];
		boolean[] nodeTypeRight2 = it2.nodeType[RIGHT];
		boolean[] nodeTypeHeavy1 = it1.nodeType[HEAVY];
		boolean[] nodeTypeHeavy2 = it2.nodeType[HEAVY];
		int[] post2size1 = it1.info[POST2_SIZE];
		int[] post2size2 = it2.info[POST2_SIZE];
		int[] post2descSum1 = it1.info[POST2_DESC_SUM];
		int[] post2descSum2 = it2.info[POST2_DESC_SUM];
		int[] post2krSum1 = it1.info[POST2_KR_SUM];
		int[] post2krSum2 = it2.info[POST2_KR_SUM];
		int[] post2revkrSum1 = it1.info[POST2_REV_KR_SUM];
		int[] post2revkrSum2 = it2.info[POST2_REV_KR_SUM];
		int[] post2parent1 = it1.info[POST2_PARENT];
		int[] post2parent2 = it2.info[POST2_PARENT];

		str = new int[size1][size2];

		// v represents nodes of left input tree in postorder
		// w represents nodes of right input tree in postorder
		for (int v = 0; v < size1; v++) {
			Arrays.fill(costW[0], 0);
			Arrays.fill(costW[1], 0);
			Arrays.fill(costW[2], 0);
			for (int w = 0; w < size2; w++) {
				if (post2size2[w] == 1) {
					// putTree zeros into arrays
					costW[LEFT][w] = 0;
					costW[RIGHT][w] = 0;
					costW[HEAVY][w] = 0;
				}
				if (post2size1[v] == 1) {
					// putTree zeros into arrays
					costV[LEFT][v][w] = 0;
					costV[RIGHT][v][w] = 0;
					costV[HEAVY][v][w] = 0;
				}

				// TODO: some things below may be putTree to outer loop

				// count the minimum + get the strategy
				heavyMin = (long) post2size1[v] * (long) post2descSum2[w]
						+ costV[HEAVY][v][w];
				revHeavyMin = (long) post2size2[w] * (long) post2descSum1[v]
						+ costW[HEAVY][w];

				leftMin = (long) post2size1[v] * (long) post2krSum2[w]
						+ costV[LEFT][v][w];
				revLeftMin = (long) post2size2[w] * (long) post2krSum1[v]
						+ costW[LEFT][w];

				rightMin = (long) post2size1[v] * (long) post2revkrSum2[w]
						+ costV[RIGHT][v][w];
				revRightMin = (long) post2size2[w] * (long) post2revkrSum1[v]
						+ costW[RIGHT][w];

				long[] mins = { leftMin, rightMin, heavyMin, Long.MAX_VALUE,
						revLeftMin, revRightMin, revHeavyMin };

				min = leftMin;
				strategy = 0;
				for (int i = 1; i <= 6; i++) {
					if (mins[i] < min) {
						min = mins[i];
						strategy = i;
					}
				}

				// store the strategy for the minimal cost
				str[v][w] = strategy;

				// fill the cost arrays
				parent1 = post2parent1[v];
				if (parent1 != -1) {
					costV[HEAVY][parent1][w] += nodeTypeHeavy1[v] ? costV[HEAVY][v][w]
							: min;
					costV[RIGHT][parent1][w] += nodeTypeRigt1[v] ? costV[RIGHT][v][w]
							: min;
					costV[LEFT][parent1][w] += nodeTypeLeft1[v] ? costV[LEFT][v][w]
							: min;
				}
				parent2 = post2parent2[w];
				if (parent2 != -1) {
					costW[HEAVY][parent2] += nodeTypeHeavy2[w] ? costW[HEAVY][w]
							: min;
					costW[LEFT][parent2] += nodeTypeLeft2[w] ? costW[LEFT][w]
							: min;
					costW[RIGHT][parent2] += nodeTypeRight2[w] ? costW[RIGHT][w]
							: min;
				}
			}
		}
	}
893788
public void shuttlesort(int from[], int to[], int low, int high) {
        if (high - low < 2) {
            return;
        }
        int middle = (low + high) / 2;
        shuttlesort(to, from, low, middle);
        shuttlesort(to, from, middle, high);

        int p = low;
        int q = middle;

        /* This is an optional short-cut; at each recursive call,
        check to see if the elements in this subset are already
        ordered.  If so, no further comparisons are needed; the
        sub-array can just be copied.  The array must be copied rather
        than assigned otherwise sister calls in the recursion might
        get out of sync.  When the number of elements is three they
        are partitioned so that the first set, [low, mid), has one
        element and and the second, [mid, high), has two. We skip the
        optimisation when the number of elements is three or less as
        the first compare in the normal merge will produce the same
        sequence of steps. This optimisation seems to be worthwhile
        for partially ordered lists but some analysis is needed to
        find out how the performance drops to Nlog(N) as the initial
        order diminishes - it may drop very quickly.  */

        if (high - low >= 4 && compare(from[middle - 1], from[middle]) <= 0) {
            System.arraycopy(from, low, to, low, high - low);
            return;
        }

        // A normal merge.

        for (int i = low; i < high; i++) {
            if (q >= high || (p < middle && compare(from[p], from[q]) <= 0)) {
                to[i] = from[p++];
            } else {
                to[i] = from[q++];
            }
        }
    }
160264
private float iterationSorSafe(GrayF32 image1, int x, int y, int pixelIndex) {
		float w = SOR_RELAXATION;

		float uf;
		float vf;
		float ui = initFlowX.data[pixelIndex];
		float vi = initFlowY.data[pixelIndex];

		float u = flowX.data[pixelIndex];
		float v = flowY.data[pixelIndex];

		float I1 = image1.data[pixelIndex];
		float I2 = warpImage2.data[pixelIndex];

		float I2x = warpDeriv2X.data[pixelIndex];
		float I2y = warpDeriv2Y.data[pixelIndex];

		float AU = A_safe(x,y,flowX);
		float AV = A_safe(x,y,flowY);

		flowX.data[pixelIndex] = uf = (1-w)*u + w*((I1-I2+I2x*ui - I2y*(v-vi))*I2x + alpha2*AU)/(I2x*I2x + alpha2);
		flowY.data[pixelIndex] = vf = (1-w)*v + w*((I1-I2+I2y*vi - I2x*(uf-ui))*I2y + alpha2*AV)/(I2y*I2y + alpha2);

		return (uf - u)*(uf - u) + (vf - v)*(vf - v);
	}
937710
def optimize_auto(self,max_iters=10000,verbose=True):
        
        self.Z.fix(warning=False)
        self.kern.fix(warning=False)
        self.kern_row.fix(warning=False)
        self.Zr.fix(warning=False)
        self.Xr.fix(warning=False)
        self.optimize(max_iters=int(0.1*max_iters),messages=verbose)
        self.unfix()
        self.optimize(max_iters=max_iters,messages=verbose)
459961
static void computePPXP( DMatrixRMaj P , DMatrixRMaj P_plus, DMatrixRMaj X , int offsetX,  DMatrixRMaj output ) {
		final int N = P.numRows;
		output.reshape(N,2);

		// diag(X) * P <-- N x 2
		for (int i = 0, index=0; i < N; i++, index += 2) {
			double x = -X.data[index+offsetX];
			output.data[index] = x*P.data[index];
			output.data[index+1] = x*P.data[index+1];
		}

		// A = P_plus*( diag(x)*P ) <-- 2 x 2
		double a00=0,a01=0,a10=0,a11=0;
		for (int i = 0, index=0; i < N; i++, index += 2) {
			a00 += P_plus.data[i]*output.data[index];
			a01 += P_plus.data[i]*output.data[index+1];
			a10 += P_plus.data[i+N]*output.data[index];
			a11 += P_plus.data[i+N]*output.data[index+1];
		}

		// P*A <-- N x 2
		for (int i = 0, index=0; i < N; i++, index += 2) {
			output.data[index]   = P.data[index]*a00 + P.data[index+1]*a10;
			output.data[index+1] = P.data[index]*a01+ P.data[index+1]*a11;
		}
	}
989483
public void process(ImagePyramid<GrayF32> image1 , ImagePyramid<GrayF32> image2 )
	{
		// Process the pyramid from low resolution to high resolution
		boolean first = true;
		for( int i = image1.getNumLayers()-1; i >= 0; i-- ) {
			GrayF32 layer1 = image1.getLayer(i);
			GrayF32 layer2 = image2.getLayer(i);

			resizeForLayer(layer1.width,layer2.height);

			// compute image derivatives
			gradient.process(layer1,deriv1X,deriv1Y);
			gradient.process(layer2,deriv2X,deriv2Y);
			hessian.process(deriv2X,deriv2Y,deriv2XX,deriv2YY,deriv2XY);

			if( !first ) {
				// interpolate initial flow from previous layer
				interpolateFlowScale(layer1.width, layer1.height);
			} else {
				// for the very first layer there is no information on flow so set everything to 0
				first = false;

				flowU.reshape(layer1.width, layer1.height);
				flowV.reshape(layer1.width, layer1.height);

				ImageMiscOps.fill(flowU,0);
				ImageMiscOps.fill(flowV,0);
			}

			// compute flow for this layer
			processLayer(layer1,layer2,deriv1X,deriv1Y,deriv2X,deriv2Y,deriv2XX,deriv2YY,deriv2XY);
		}
	}
1373057
private void score0(double[] data, double[] preds, int treeIdx) {
    Key[] keys = _output._treeKeys[treeIdx];
    for( int c=0; c<keys.length; c++ ) {
      if (keys[c] != null) {
        double pred = DKV.get(keys[c]).<CompressedTree>get().score(data,_output._domains);
        assert (!Double.isInfinite(pred));
        preds[keys.length == 1 ? 0 : c + 1] += pred;
      }
    }
  }
449152
def multifuncrs(X, funs, eps=1E-6,
                nswp=10,
                kickrank=5,
                y0=None,
                rmax=999999,  # TODO:infinity \
                kicktype='amr-two',        \
                pcatype='svd',             \
                trunctype='fro',           \
                d2=1,                      \
                do_qr=False,               \
                verb=1):
    

    dtype = np.float64
    if len([x for x in X if x.is_complex]) > 0:
        dtype = np.complex128

    y = y0
    wasrand = False

    nx = len(X)
    d = X[0].d
    n = X[0].n
    rx = np.transpose(np.array([ttx.r for ttx in X]))
    #crx = [tt.tensor.to_list(ttx) for x in X]
    #crx = zip(*crx)
    crx = np.transpose(np.array([tt.tensor.to_list(ttx)
                                 for ttx in X], dtype=np.object))
    crx = np.empty((nx, d), dtype=np.object)
    i = 0
    for ttx in X:
        v = tt.tensor.to_list(ttx)
        j = 0
        for w in v:
            crx[i, j] = w
            j = j + 1
        i = i + 1
    crx = crx.T
    if y is None:
        ry = d2 * np.ones((d + 1,), dtype=np.int32)
        ry[0] = 1
        y = tt.rand(n, d, ry)
        wasrand = True

    ry = y.r
    cry = tt.tensor.to_list(y)

    Ry = np.zeros((d + 1, ), dtype=np.object)
    Ry[0] = np.array([[1.0]], dtype=dtype)
    Ry[d] = np.array([[1.0]], dtype=dtype)
    Rx = np.zeros((d + 1, nx), dtype=np.object)
    Rx[0, :] = np.ones(nx, dtype=dtype)
    Rx[d, :] = np.ones(nx, dtype=dtype)

    block_order = [+d, -d]

    # orth
    for i in range(0, d - 1):
        cr = cry[i]
        cr = reshape(cr, (ry[i] * n[i], ry[i + 1]))
        cr, rv = np.linalg.qr(cr)
        cr2 = cry[i + 1]
        cr2 = reshape(cr2, (ry[i + 1], n[i + 1] * ry[i + 2]))
        cr2 = np.dot(rv, cr2)  # matrix multiplication
        ry[i + 1] = cr.shape[1]
        cr = reshape(cr, (ry[i], n[i], ry[i + 1]))
        cry[i + 1] = reshape(cr2, (ry[i + 1], n[i + 1], ry[i + 2]))
        cry[i] = cr

        Ry[i + 1] = np.dot(Ry[i], reshape(cr, (ry[i], n[i] * ry[i + 1])))
        Ry[i + 1] = reshape(Ry[i + 1], (ry[i] * n[i], ry[i + 1]))
        curind = []
        if wasrand:
            # EVERY DAY I'M SHUFFLIN'
            curind = np.random.permutation(n[i] * ry[i])[:ry[i + 1]]
        else:
            curind = maxvol(Ry[i + 1])
        Ry[i + 1] = Ry[i + 1][curind, :]
        for j in range(0, nx):
            try:
                Rx[i + 1, j] = reshape(crx[i, j],
                                       (rx[i, j], n[i] * rx[i + 1, j]))
            except:
                pass
            Rx[i + 1, j] = np.dot(Rx[i, j], Rx[i + 1, j])
            Rx[i + 1, j] = reshape(Rx[i + 1, j], (ry[i] * n[i], rx[i + 1, j]))
            Rx[i + 1, j] = Rx[i + 1, j][curind, :]

    d2 = ry[d]
    ry[d] = 1
    cry[d - 1] = np.transpose(cry[d - 1], [2, 0, 1])  # permute

    last_sweep = False
    swp = 1

    dy = np.zeros((d, ))
    max_dy = 0

    cur_order = copy.copy(block_order)
    order_index = 1
    i = d - 1
    # can't use 'dir' identifier in python
    dirn = int(math.copysign(1, cur_order[order_index]))

    # DMRG sweeps
    while swp <= nswp or dirn > 0:

        oldy = reshape(cry[i], (d2 * ry[i] * n[i] * ry[i + 1],))

        if not last_sweep:
            # compute the X superblocks
            curbl = np.zeros((ry[i] * n[i] * ry[i + 1], nx), dtype=dtype)
            for j in range(0, nx):
                cr = reshape(crx[i, j], (rx[i, j], n[i] * rx[i + 1, j]))
                cr = np.dot(Rx[i, j], cr)
                cr = reshape(cr, (ry[i] * n[i], rx[i + 1, j]))
                cr = np.dot(cr, Rx[i + 1, j])
                curbl[:, j] = cr.flatten('F')
            # call the function
            newy = funs(curbl)
            # multiply with inverted Ry
            newy = reshape(newy, (ry[i], n[i] * ry[i + 1] * d2))
            newy = np.linalg.solve(Ry[i], newy)  # y = R \ y
            newy = reshape(newy, (ry[i] * n[i] * ry[i + 1], d2))
            newy = reshape(np.transpose(newy), (d2 * ry[i] * n[i], ry[i + 1]))
            newy = np.transpose(np.linalg.solve(
                np.transpose(Ry[i + 1]), np.transpose(newy)))  # y=y/R
            newy = reshape(newy, (d2 * ry[i] * n[i] * ry[i + 1],))
        else:
            newy = oldy

        dy[i] = np.linalg.norm(newy - oldy) / np.linalg.norm(newy)
        max_dy = max(max_dy, dy[i])

        # truncation
        if dirn > 0:  # left-to-right
            newy = reshape(newy, (d2, ry[i] * n[i] * ry[i + 1]))
            newy = reshape(np.transpose(newy), (ry[i] * n[i], ry[i + 1] * d2))
        else:
            newy = reshape(newy, (d2 * ry[i], n[i] * ry[i + 1]))

        r = 0  # defines a variable in global scope

        if kickrank >= 0:
            u, s, v = np.linalg.svd(newy, full_matrices=False)
            v = np.conj(np.transpose(v))
            if trunctype == "fro" or last_sweep:
                r = my_chop2(s, eps / math.sqrt(d) * np.linalg.norm(s))
            else:
                # truncate taking into account the (r+1) overhead in the cross
                # (T.S.: what?)
                cums = abs(s * np.arange(2, len(s) + 2)) ** 2
                cums = np.cumsum(cums[::-1])[::-1]
                cums = cums / cums[0]
                ff = [i for i in range(len(cums)) if cums[i] < eps ** 2 / d]
                if len(ff) == 0:
                    r = len(s)
                else:
                    r = np.amin(ff)
            r = min(r, rmax, len(s))
        else:
            if dirn > 0:
                u, v = np.linalg.qr(newy)
                v = np.conj(np.transpose(v))
                r = u.shape[1]
                s = np.ones((r, ))
            else:
                v, u = np.linalg.qr(np.transpose(newy))
                v = np.conj(v)
                u = np.transpose(u)
                r = u.shape[1]
                s = np.ones((r, ))

        if verb > 1:
            print('=multifuncrs=   block %d{%d}, dy: %3.3e, r: %d' % (i, dirn, dy[i], r))

        # kicks and interfaces
        if dirn > 0 and i < d - 1:
            u = u[:, :r]
            v = np.dot(v[:, :r], np.diag(s[:r]))

            # kick
            radd = 0
            rv = 1
            if not last_sweep and kickrank > 0:
                uk = None
                if kicktype == 'amr-two':
                    # AMR(two)-like kick.

                    # compute the X superblocks
                    ind2 = np.unique(np.random.randint(
                        0, ry[i + 2] * n[i + 1], ry[i + 1]))
                    #ind2 = np.unique(np.floor(np.random.rand(ry[i + 1]) * (ry[i + 2] * n[i + 1])))
                    rkick = len(ind2)
                    curbl = np.zeros((ry[i] * n[i] * rkick, nx), dtype=dtype)
                    for j in range(nx):
                        cr1 = reshape(
                            crx[i, j], (rx[i, j], n[i] * rx[i + 1, j]))
                        cr1 = np.dot(Rx[i, j], cr1)
                        cr1 = reshape(cr1, (ry[i] * n[i], rx[i + 1, j]))
                        cr2 = reshape(
                            crx[i + 1, j], (rx[i + 1, j] * n[i + 1], rx[i + 2, j]))
                        cr2 = np.dot(cr2, Rx[i + 2, j])
                        cr2 = reshape(
                            cr2, (rx[i + 1, j], n[i + 1] * ry[i + 2]))
                        cr2 = cr2[:, ind2]
                        curbl[:, j] = reshape(
                            np.dot(cr1, cr2), (ry[i] * n[i] * rkick,))
                    # call the function
                    uk = funs(curbl)
                    uk = reshape(uk, (ry[i], n[i] * rkick * d2))
                    uk = np.linalg.solve(Ry[i], uk)
                    uk = reshape(uk, (ry[i] * n[i], rkick * d2))
                    if pcatype == 'svd':
                        uk, sk, vk = np.linalg.svd(uk, full_matrices=False)
                        vk = np.conj(np.transpose(vk))
                        uk = uk[:, :min(kickrank, uk.shape[1])]
                    else:
                        # uk = uchol(np.transpose(uk), kickrank + 1) # TODO
                        uk = uk[:, :max(uk.shape[1] - kickrank + 1, 1):-1]
                else:
                    uk = np.random.rand(ry[i] * n[i], kickrank)
                u, rv = np.linalg.qr(np.concatenate((u, uk), axis=1))
                radd = uk.shape[1]
            v = np.concatenate(
                (v, np.zeros((ry[i + 1] * d2, radd), dtype=dtype)), axis=1)
            v = np.dot(rv, np.conj(np.transpose(v)))
            r = u.shape[1]

            cr2 = cry[i + 1]
            cr2 = reshape(cr2, (ry[i + 1], n[i + 1] * ry[i + 2]))
            v = reshape(v, (r * ry[i + 1], d2))
            v = reshape(np.transpose(v), (d2 * r, ry[i + 1]))
            v = np.dot(v, cr2)

            ry[i + 1] = r

            u = reshape(u, (ry[i], n[i], r))
            v = reshape(v, (d2, r, n[i + 1], ry[i + 2]))

            cry[i] = u
            cry[i + 1] = v

            Ry[i + 1] = np.dot(Ry[i], reshape(u, (ry[i], n[i] * ry[i + 1])))
            Ry[i + 1] = reshape(Ry[i + 1], (ry[i] * n[i], ry[i + 1]))
            curind = maxvol(Ry[i + 1])
            Ry[i + 1] = Ry[i + 1][curind, :]
            for j in range(nx):
                Rx[i + 1, j] = reshape(crx[i, j],
                                       (rx[i, j], n[i] * rx[i + 1, j]))
                Rx[i + 1, j] = np.dot(Rx[i, j], Rx[i + 1, j])
                Rx[i + 1, j] = reshape(Rx[i + 1, j],
                                       (ry[i] * n[i], rx[i + 1, j]))
                Rx[i + 1, j] = Rx[i + 1, j][curind, :]
        elif dirn < 0 and i > 0:
            u = np.dot(u[:, :r], np.diag(s[:r]))
            v = np.conj(v[:, :r])

            radd = 0
            rv = 1
            if not last_sweep and kickrank > 0:
                if kicktype == 'amr-two':
                    # compute the X superblocks
                    ind2 = np.unique(np.random.randint(
                        0, ry[i - 1] * n[i - 1], ry[i]))
                    rkick = len(ind2)
                    curbl = np.zeros(
                        (rkick * n[i] * ry[i + 1], nx), dtype=dtype)
                    for j in range(nx):
                        cr1 = reshape(
                            crx[i, j], (rx[i, j] * n[i], rx[i + 1, j]))
                        cr1 = np.dot(cr1, Rx[i + 1, j])
                        cr1 = reshape(cr1, (rx[i, j], n[i] * ry[i + 1]))
                        cr2 = reshape(
                            crx[i - 1, j], (rx[i - 1, j], n[i - 1] * rx[i, j]))
                        cr2 = np.dot(Rx[i - 1, j], cr2)
                        cr2 = reshape(cr2, (ry[i - 1] * n[i - 1], rx[i, j]))
                        cr2 = cr2[ind2, :]
                        curbl[:, j] = reshape(
                            np.dot(cr2, cr1), (rkick * n[i] * ry[i + 1],))
                    # calling the function
                    uk = funs(curbl)
                    uk = reshape(uk, (rkick * n[i] * ry[i + 1], d2))
                    uk = reshape(np.transpose(
                        uk), (d2 * rkick * n[i], ry[i + 1]))
                    uk = np.transpose(np.linalg.solve(
                        np.transpose(Ry[i + 1]), np.transpose(uk)))
                    uk = reshape(uk, (d2 * rkick, n[i] * ry[i + 1]))
                    if pcatype == 'svd':
                        vk, sk, uk = np.linalg.svd(uk, full_matrices=False)
                        uk = np.conj(np.transpose(uk))
                        # TODO: refactor
                        uk = uk[:, :min(kickrank, uk.shape[1])]
                    else:
                        # uk = uchol(uk, kickrank + 1) # TODO
                        uk = uk[:, :max(uk.shape[1] - kickrank + 1, 1):-1]
                else:
                    uk = np.random.rand(n[i] * ry[i + 1], kickrank)
                v, rv = np.linalg.qr(np.concatenate((v, uk), axis=1))
                radd = uk.shape[1]
            u = np.concatenate(
                (u, np.zeros((d2 * ry[i], radd), dtype=dtype)), axis=1)
            u = np.dot(u, np.transpose(rv))
            r = v.shape[1]
            cr2 = cry[i - 1]
            cr2 = reshape(cr2, (ry[i - 1] * n[i - 1], ry[i]))
            u = reshape(u, (d2, ry[i] * r))
            u = reshape(np.transpose(u), (ry[i], r * d2))
            u = np.dot(cr2, u)

            u = reshape(u, (ry[i - 1] * n[i - 1] * r, d2))
            u = reshape(np.transpose(u), (d2, ry[i - 1], n[i - 1], r))
            v = reshape(np.transpose(v), (r, n[i], ry[i + 1]))

            ry[i] = r
            cry[i - 1] = u
            cry[i] = v

            Ry[i] = np.dot(reshape(v, (ry[i] * n[i], ry[i + 1])), Ry[i + 1])
            Ry[i] = reshape(Ry[i], (ry[i], n[i] * ry[i + 1]))
            curind = maxvol(np.transpose(Ry[i]))
            Ry[i] = Ry[i][:, curind]
            for j in range(nx):
                Rx[i, j] = reshape(crx[i, j], (rx[i, j] * n[i], rx[i + 1, j]))
                Rx[i, j] = np.dot(Rx[i, j], Rx[i + 1, j])
                Rx[i, j] = reshape(Rx[i, j], (rx[i, j], n[i] * ry[i + 1]))
                Rx[i, j] = Rx[i, j][:, curind]
        elif dirn > 0 and i == d - 1:
            newy = np.dot(np.dot(u[:, :r], np.diag(s[:r])),
                          np.conj(np.transpose(v[:, :r])))
            newy = reshape(newy, (ry[i] * n[i] * ry[i + 1], d2))
            cry[i] = reshape(np.transpose(newy), (d2, ry[i], n[i], ry[i + 1]))
        elif dirn < 0 and i == 0:
            newy = np.dot(np.dot(u[:, :r], np.diag(s[:r])),
                          np.conj(np.transpose(v[:, :r])))
            newy = reshape(newy, (d2, ry[i], n[i], ry[i + 1]))
            cry[i] = newy

        i = i + dirn
        cur_order[order_index] = cur_order[order_index] - dirn
        if cur_order[order_index] == 0:
            order_index = order_index + 1
            if verb > 0:
                print('=multifuncrs= sweep %d{%d}, max_dy: %3.3e, erank: %g' % (swp, order_index, max_dy,
                                                                                math.sqrt(np.dot(ry[:d], n * ry[1:]) / np.sum(n))))

            if last_sweep:
                break
            if max_dy < eps and dirn < 0:
                last_sweep = True
                kickrank = 0

            if order_index >= len(cur_order):
                cur_order = copy.copy(block_order)
                order_index = 0
                if last_sweep:
                    cur_order = [d - 1]

                max_dy = 0
                swp = swp + 1

            dirn = int(math.copysign(1, cur_order[order_index]))
            i = i + dirn

    cry[d - 1] = np.transpose(cry[d - 1][:, :, :, 0], [1, 2, 0])
    y = tt.tensor.from_list(cry)
    return y
1158312
private static void compute(Kernel kernel, int[] in, int[] out, int width, int height, boolean alpha, int edge)
    {
        final float[] matrix = kernel.getMatrix();
        final int cols = kernel.getWidth();
        final int cols2 = cols / 2;

        for (int y = 0; y < height; y++)
        {
            int index = y;
            final int ioffset = y * width;
            for (int x = 0; x < width; x++)
            {
                compute(matrix, in, out, x, index, ioffset, cols2, width, alpha, edge);
                index += height;
            }
        }
    }
1366201
def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
        
        sp, sw = 0., 0.
        upb = lwinfo.max() + self.wgt[t].lw.max()
        if hasattr(self.model, 'upper_bound_trans'):
            upb += self.model.upper_bound_trans(t + 1)
        # Loop over n, to avoid having in memory a NxN matrix
        for n in range(self.N):
            omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                            + self.model.logpt(t + 1, self.X[t][n],
                                               info.hist.X[ti]))
            sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
            sw += np.sum(omegan)
        return sp / sw
931314
def bitwise_or(self, t):
        
        """
        This implementation combines the approaches used by 'WYSINWYX: what you see is not what you execute'
        paper and 'Signedness-Agnostic Program Analysis: Precise Integer Bounds for Low-Level Code'. The
        first paper provides an sound way to approximate the stride, whereas the second provides a way
        to calculate the or operation using wrapping intervals.
        Note that, even though according Warren's work 'Hacker's delight', one should follow different
        approaches to calculate the minimun/maximum values of an or operations according on the type
        of the operands (signed/unsigned). On the other other hand, by splitting the wrapping-intervals
        at the south pole, we can safely and soundly only use the Warren's functions for unsigned
        integers.
        """
        s = self
        result_interval = list()

        for u in s._ssplit():
            for v in t._ssplit():
                w = u.bits
                # u |w v
                if u.is_integer:
                    s_t = StridedInterval._ntz(v.stride)
                elif v.is_integer:
                    s_t = StridedInterval._ntz(u.stride)
                else:
                    s_t = min(StridedInterval._ntz(u.stride), StridedInterval._ntz(v.stride))

                if u.is_integer and u.lower_bound == 0:
                    new_stride = v.stride
                elif v.is_integer and v.lower_bound == 0:
                    new_stride = u.stride
                else:
                    new_stride = 2 ** s_t
                mask = (1 << s_t) - 1
                r = (u.lower_bound & mask) | (v.lower_bound & mask)
                m = (2 ** w) - 1
                low_bound = WarrenMethods.min_or(u.lower_bound & (~mask & m), u.upper_bound & (~mask & m), v.lower_bound & (~mask & m), v.upper_bound & (~mask & m), w)
                upper_bound = WarrenMethods.max_or(u.lower_bound & (~mask & m), u.upper_bound & (~mask & m), v.lower_bound & (~mask & m), v.upper_bound & (~mask & m), w)
                if low_bound == upper_bound:
                    new_stride = 0

                new_interval = StridedInterval(lower_bound=((low_bound & (~mask & m)) | r), upper_bound=((upper_bound & (~mask & m)) | r), bits=w, stride=new_stride)
                result_interval.append(new_interval)

        return StridedInterval.least_upper_bound(*result_interval).normalize()
1276795
public void process( T image1 , T image2 )
	{
		// declare image data structures
		if( pyr1 == null || pyr1.getInputWidth() != image1.width || pyr1.getInputHeight() != image1.height ) {
			pyr1 = UtilDenseOpticalFlow.standardPyramid(image1.width, image1.height, scale, sigma, 5, maxLayers, GrayF32.class);
			pyr2 = UtilDenseOpticalFlow.standardPyramid(image1.width, image1.height, scale, sigma, 5, maxLayers, GrayF32.class);

			pyr1.initialize(image1.width,image1.height);
			pyr2.initialize(image1.width,image1.height);
		}

		norm1.reshape(image1.width, image1.height);
		norm2.reshape(image1.width, image1.height);

		// normalize input image to make sure alpha is image independent
		imageNormalization(image1, image2, norm1, norm2);

		// create image pyramid
		pyr1.process(norm1);
		pyr2.process(norm2);

		// compute flow from pyramid
		process(pyr1, pyr2);
	}
1125860
private Pair<Gradient, INDArray> getGradientsAndDelta(INDArray preOut, LayerWorkspaceMgr workspaceMgr) {
        ILossFunction lossFunction = layerConf().getLossFn();
        INDArray labels2d = getLabels2d(workspaceMgr, ArrayType.BP_WORKING_MEM);
        INDArray delta = lossFunction.computeGradient(labels2d, preOut, layerConf().getActivationFn(), maskArray);
        org.deeplearning4j.nn.conf.ocnn.OCNNOutputLayer conf = ( org.deeplearning4j.nn.conf.ocnn.OCNNOutputLayer) conf().getLayer();


        if(conf.getLastEpochSinceRUpdated() == 0 && epochCount == 0) {
            INDArray currentR = doOutput(false,workspaceMgr);
            if(window == null) {
                window = Nd4j.createUninitializedDetached(conf.getWindowSize()).assign(0.0);
            }

            if(batchWindowSizeIndex < window.length() - currentR.length()) {
                window.put(new INDArrayIndex[]{NDArrayIndex.interval(batchWindowSizeIndex,batchWindowSizeIndex + currentR.length())},currentR);
            }
            else if(batchWindowSizeIndex < window.length()) {
                int windowIdx = (int) window.length() - batchWindowSizeIndex;
                window.put(new INDArrayIndex[]{NDArrayIndex.interval(window.length() - windowIdx,window.length())},currentR.get(NDArrayIndex.interval(0,windowIdx)));

            }

            batchWindowSizeIndex += currentR.length();
            conf.setLastEpochSinceRUpdated(epochCount);
        }
        else if(conf.getLastEpochSinceRUpdated()  != epochCount) {
            double percentile = window.percentileNumber(100.0 * conf.getNu()).doubleValue();
            getParam(R_KEY).putScalar(0,percentile);
            conf.setLastEpochSinceRUpdated(epochCount);
            batchWindowSizeIndex = 0;
        }
        else {
            //track a running average per minibatch per epoch
            //calculate the average r value quantl=ile
            //once the epoch changes

            INDArray currentR = doOutput(false,workspaceMgr);
            window.put(new INDArrayIndex[]{NDArrayIndex.interval(batchWindowSizeIndex,batchWindowSizeIndex + currentR.length())},currentR);
        }


        Gradient gradient = new DefaultGradient();
        INDArray vGradView = gradientViews.get(V_KEY);
        double oneDivNu = 1.0 / layerConf().getNu();
        INDArray xTimesV = input.mmul(getParam(V_KEY));
        INDArray derivW = layerConf().getActivationFn().getActivation(xTimesV.dup(),true).negi();
        derivW = derivW.muliColumnVector(delta).mean(0).muli(oneDivNu).addi(getParam(W_KEY));
        gradient.setGradientFor(W_KEY,gradientViews.get(W_KEY).assign(derivW));

        //dG -> sigmoid derivative

        INDArray firstVertDerivV =  layerConf().getActivationFn()
                .backprop(xTimesV.dup(),Nd4j.ones(input.dataType(), xTimesV.shape()))
                .getFirst().muliRowVector(getParam(W_KEY).neg());
        firstVertDerivV = firstVertDerivV.muliColumnVector(delta)
                        .reshape('f',input.size(0),1,layerConf().getHiddenSize());
        INDArray secondTermDerivV = input.reshape('f',
                input.size(0),getParam(V_KEY).size(0),1);

        long[]  shape = new long[firstVertDerivV.shape().length];
        for(int i = 0; i < firstVertDerivV.rank(); i++) {
            shape[i] = Math.max(firstVertDerivV.size(i),secondTermDerivV.size(i));
        }

        INDArray firstDerivVBroadcast = Nd4j.createUninitialized(input.dataType(), shape);

        INDArray mulResult = firstVertDerivV.broadcast(firstDerivVBroadcast);
        int[] bcDims = {0,1};
        Broadcast.mul(mulResult, secondTermDerivV, mulResult, bcDims);

        INDArray derivV = mulResult
                .mean(0).muli(oneDivNu).addi(getParam(V_KEY));
        gradient.setGradientFor(V_KEY,vGradView.assign(derivV));



        INDArray derivR = Nd4j.scalar(delta.meanNumber()).muli(oneDivNu).addi(-1);
        gradient.setGradientFor(R_KEY,gradientViews.get(R_KEY).assign(derivR));
        clearNoiseWeightParams();

        delta = backpropDropOutIfPresent(delta);
        return new Pair<>(gradient, delta);
    }
208523
private double optimize(int k0, int kmax, double sumx, double sumx2, double xp, double yp, double sumxy, double sumy) {
    int k = kmax - k0 + 1;
    return (sumxy - xp * sumy - yp * sumx + k * xp * yp) / (sumx2 - 2 * sumx * xp + k * xp * xp);
    // return (-xp * yp * + yp * sumx - sumxy + xp * sumy) / (-xp * xp *
    // kmax - sumx2 + 2 * xp * sumx);
  }
752390
private double update(double[][] data, double T, int k, double[][] centroids, double[][] posteriori, double[] priori) {
        int n = data.length;
        int d = data[0].length;

        double D = 0.0;
        double H = 0.0;

        int iter = 0;
        double currentDistortion = Double.MAX_VALUE;
        double newDistortion = Double.MAX_VALUE / 2;
        while (iter < 100 && currentDistortion > newDistortion) {
            currentDistortion = newDistortion;

            D = Double.NaN;
            H = 0.0;

            if (tasks != null) {
                try {
                    D = 0.0;
                    for (UpdateThread t : tasks) {
                        t.k = k;
                        t.T = T;
                    }
                    
                    for (UpdateThread t : MulticoreExecutor.run(tasks)) {
                        D += t.D;
                        H += t.H;
                    }
                } catch (Exception ex) {
                    logger.error("Failed to run Deterministic Annealing on multi-core", ex);

                    D = Double.NaN;
                }
            }

            if (Double.isNaN(D)) {
                D = 0.0;
                double[] dist = new double[k];

                for (int i = 0; i < n; i++) {
                    double p = 0.0;

                    for (int j = 0; j < k; j++) {
                        dist[j] = Math.squaredDistance(data[i], centroids[j]);
                        posteriori[i][j] = priori[j] * Math.exp(-dist[j] / T);
                        p += posteriori[i][j];
                    }

                    double r = 0.0;
                    for (int j = 0; j < k; j++) {
                        posteriori[i][j] /= p;
                        D += posteriori[i][j] * dist[j];
                        r += -posteriori[i][j] * Math.log(posteriori[i][j]);
                    }
                    H += r;
                }
            }
            
            for (int i = 0; i < k; i++) {
                priori[i] = 0;
                for (int j = 0; j < n; j++) {
                    priori[i] += posteriori[j][i];
                }
                priori[i] /= n;
            }

            boolean parallel = false;
            if (ctasks != null) {
                try {
                    for (CentroidThread t : ctasks) {
                        t.k = k;
                    }
                    
                    MulticoreExecutor.run(ctasks);
                    parallel = true;
                } catch (Exception ex) {
                    logger.error("Failed to run Deterministic Annealing on multi-core", ex);

                    parallel = false;
                }
            }

            if (!parallel) {
                for (int i = 0; i < k; i++) {
                    Arrays.fill(centroids[i], 0.0);
                    for (int j = 0; j < d; j++) {
                        for (int m = 0; m < n; m++) {
                            centroids[i][j] += data[m][j] * posteriori[m][i];
                        }
                        centroids[i][j] /= (n * priori[i]);
                    }
                }
            }

            newDistortion = D - T * H;
            iter++;
        }

        logger.info(String.format("Deterministic Annealing clustering entropy after %3d iterations at temperature %.4f and k = %d: %.5f (soft distortion = %.5f )%n", iter, T, k / 2, H, D));

        return currentDistortion;
    }