4D simplex noise, better simplex rank ordering method 2012-03-09
347096
function(x, y, z, w) {
      var perm = this.perm;
      var grad4 = this.grad4;

      var n0, n1, n2, n3, n4; // Noise contributions from the five corners
      // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
      var s = (x + y + z + w) * F4; // Factor for 4D skewing
      var i = Math.floor(x + s);
      var j = Math.floor(y + s);
      var k = Math.floor(z + s);
      var l = Math.floor(w + s);
      var t = (i + j + k + l) * G4; // Factor for 4D unskewing
      var X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
      var Y0 = j - t;
      var Z0 = k - t;
      var W0 = l - t;
      var x0 = x - X0; // The x,y,z,w distances from the cell origin
      var y0 = y - Y0;
      var z0 = z - Z0;
      var w0 = w - W0;
      // For the 4D case, the simplex is a 4D shape I won't even try to describe.
      // To find out which of the 24 possible simplices we're in, we need to
      // determine the magnitude ordering of x0, y0, z0 and w0.
      // Six pair-wise comparisons are performed between each possible pair
      // of the four coordinates, and the results are used to rank the numbers.
      var rankx = 0;
      var ranky = 0;
      var rankz = 0;
      var rankw = 0;
      if (x0 > y0) rankx++;
      else ranky++;
      if (x0 > z0) rankx++;
      else rankz++;
      if (x0 > w0) rankx++;
      else rankw++;
      if (y0 > z0) ranky++;
      else rankz++;
      if (y0 > w0) ranky++;
      else rankw++;
      if (z0 > w0) rankz++;
      else rankw++;
      var i1, j1, k1, l1; // The integer offsets for the second simplex corner
      var i2, j2, k2, l2; // The integer offsets for the third simplex corner
      var i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
      // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
      // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
      // impossible. Only the 24 indices which have non-zero entries make any sense.
      // We use a thresholding to set the coordinates in turn from the largest magnitude.
      // Rank 3 denotes the largest coordinate.
      i1 = rankx >= 3 ? 1 : 0;
      j1 = ranky >= 3 ? 1 : 0;
      k1 = rankz >= 3 ? 1 : 0;
      l1 = rankw >= 3 ? 1 : 0;
      // Rank 2 denotes the second largest coordinate.
      i2 = rankx >= 2 ? 1 : 0;
      j2 = ranky >= 2 ? 1 : 0;
      k2 = rankz >= 2 ? 1 : 0;
      l2 = rankw >= 2 ? 1 : 0;
      // Rank 1 denotes the second smallest coordinate.
      i3 = rankx >= 1 ? 1 : 0;
      j3 = ranky >= 1 ? 1 : 0;
      k3 = rankz >= 1 ? 1 : 0;
      l3 = rankw >= 1 ? 1 : 0;
      // The fifth corner has all coordinate offsets = 1, so no need to compute that.
      var x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords
      var y1 = y0 - j1 + G4;
      var z1 = z0 - k1 + G4;
      var w1 = w0 - l1 + G4;
      var x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in (x,y,z,w) coords
      var y2 = y0 - j2 + 2.0 * G4;
      var z2 = z0 - k2 + 2.0 * G4;
      var w2 = w0 - l2 + 2.0 * G4;
      var x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in (x,y,z,w) coords
      var y3 = y0 - j3 + 3.0 * G4;
      var z3 = z0 - k3 + 3.0 * G4;
      var w3 = w0 - l3 + 3.0 * G4;
      var x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in (x,y,z,w) coords
      var y4 = y0 - 1.0 + 4.0 * G4;
      var z4 = z0 - 1.0 + 4.0 * G4;
      var w4 = w0 - 1.0 + 4.0 * G4;
      // Work out the hashed gradient indices of the five simplex corners
      var ii = i & 255;
      var jj = j & 255;
      var kk = k & 255;
      var ll = l & 255;
      // Calculate the contribution from the five corners
      var t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
      if (t0 < 0) n0 = 0.0;
      else {
        var gi0 = (perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32) * 4;
        t0 *= t0;
        n0 = t0 * t0 * (grad4[gi0] * x0 + grad4[gi0 + 1] * y0 + grad4[gi0 + 2] * z0 + grad4[gi0 + 3] * w0);
      }
      var t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
      if (t1 < 0) n1 = 0.0;
      else {
        var gi1 = (perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32) * 4;
        t1 *= t1;
        n1 = t1 * t1 * (grad4[gi1] * x1 + grad4[gi1 + 1] * y1 + grad4[gi1 + 2] * z1 + grad4[gi1 + 3] * w1);
      }
      var t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
      if (t2 < 0) n2 = 0.0;
      else {
        var gi2 = (perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32) * 4;
        t2 *= t2;
        n2 = t2 * t2 * (grad4[gi2] * x2 + grad4[gi2 + 1] * y2 + grad4[gi2 + 2] * z2 + grad4[gi2 + 3] * w2);
      }
      var t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
      if (t3 < 0) n3 = 0.0;
      else {
        var gi3 = (perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32) * 4;
        t3 *= t3;
        n3 = t3 * t3 * (grad4[gi3] * x3 + grad4[gi3 + 1] * y3 + grad4[gi3 + 2] * z3 + grad4[gi3 + 3] * w3);
      }
      var t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
      if (t4 < 0) n4 = 0.0;
      else {
        var gi4 = (perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32) * 4;
        t4 *= t4;
        n4 = t4 * t4 * (grad4[gi4] * x4 + grad4[gi4 + 1] * y4 + grad4[gi4 + 2] * z4 + grad4[gi4 + 3] * w4);
      }
      // Sum up and scale the result to cover the range [-1,1]
      return 27.0 * (n0 + n1 + n2 + n3 + n4);
    }
225582
function decision2DStumpTrain(data, labels, ix, options) {
    
    options = options || {};
    var numtries = options.numTries || 10;
    
    // choose a dimension at random and pick a best split
    var N = ix.length;
    
    var ri1= 0;
    var ri2= 1;
    if(data[0].length > 2) {
      // more than 2D data. Pick 2 random dimensions
      ri1= randi(0, data[0].length);
      ri2= randi(0, data[0].length);
      while(ri2 == ri1) ri2= randi(0, data[0].length); // must be distinct!
    }
    
    // evaluate class entropy of incoming data
    var H= entropy(labels, ix);
    var bestGain=0; 
    var bestw1, bestw2, bestthr;
    var dots= new Array(ix.length);
    for(var i=0;i<numtries;i++) {
        
        // pick random line parameters
        var alpha= randf(0, 2*Math.PI);
        var w1= Math.cos(alpha);
        var w2= Math.sin(alpha);
        
        // project data on this line and get the dot products
        for(var j=0;j<ix.length;j++) {
          dots[j]= w1*data[ix[j]][ri1] + w2*data[ix[j]][ri2];
        }
        
        // we are in a tricky situation because data dot product distribution
        // can be skewed. So we don't want to select just randomly between
        // min and max. But we also don't want to sort as that is too expensive
        // let's pick two random points and make the threshold be somewhere between them.
        // for skewed datasets, the selected points will with relatively high likelihood
        // be in the high-desnity regions, so the thresholds will make sense
        var ix1= ix[randi(0, N)];
        var ix2= ix[randi(0, N)];
        while(ix2==ix1) ix2= ix[randi(0, N)]; // enforce distinctness of ix2
        var a= Math.random();
        var dotthr= dots[ix1]*a + dots[ix2]*(1-a);
        
        // measure information gain we'd get from split with thr
        var l1=1, r1=1, lm1=1, rm1=1; //counts for Left and label 1, right and label 1, left and minus 1, right and minus 1
        for(var j=0;j<ix.length;j++) {
            if(dots[j] < dotthr) {
              if(labels[ix[j]]==1) l1++;
              else lm1++;
            } else {
              if(labels[ix[j]]==1) r1++;
              else rm1++;
            }
        }
        var t= l1+lm1; 
        l1=l1/t;
        lm1=lm1/t;
        t= r1+rm1;
        r1=r1/t;
        rm1= rm1/t;
        
        var LH= -l1*Math.log(l1) -lm1*Math.log(lm1); // left and right entropy
        var RH= -r1*Math.log(r1) -rm1*Math.log(rm1);
        
        var informationGain= H - LH - RH;
        //console.log("Considering split %f, entropy %f -> %f, %f. Gain %f", thr, H, LH, RH, informationGain);
        if(informationGain > bestGain || i === 0) {
            bestGain= informationGain;
            bestw1= w1;
            bestw2= w2;
            bestthr= dotthr;
        }
    }
    
    model= {};
    model.w1= bestw1;
    model.w2= bestw2;
    model.dotthr= bestthr;
    return model;
}
966214
def noise3(self, x, y, z, repeat, base=0.0):
		
		i = int(fmod(floor(x), repeat))
		j = int(fmod(floor(y), repeat))
		k = int(fmod(floor(z), repeat))
		ii = (i + 1) % repeat
		jj = (j + 1) % repeat
		kk = (k + 1) % repeat
		if base:
			i += base; j += base; k += base
			ii += base; jj += base; kk += base

		x -= floor(x); y -= floor(y); z -= floor(z)
		fx = x**3 * (x * (x * 6 - 15) + 10)
		fy = y**3 * (y * (y * 6 - 15) + 10)
		fz = z**3 * (z * (z * 6 - 15) + 10)

		perm = self.permutation
		A = perm[i]
		AA = perm[A + j]
		AB = perm[A + jj]
		B = perm[ii]
		BA = perm[B + j]
		BB = perm[B + jj]
		
		return lerp(fz, lerp(fy, lerp(fx, grad3(perm[AA + k], x, y, z),
										  grad3(perm[BA + k], x - 1, y, z)),
								 lerp(fx, grad3(perm[AB + k], x, y - 1, z),
										  grad3(perm[BB + k], x - 1, y - 1, z))),
						lerp(fy, lerp(fx, grad3(perm[AA + kk], x, y, z - 1),
										  grad3(perm[BA + kk], x - 1, y, z - 1)),
								 lerp(fx, grad3(perm[AB + kk], x, y - 1, z - 1),
										  grad3(perm[BB + kk], x - 1, y - 1, z - 1))))
712950
private List<Integer> getPermutation(Matrix4d transformation) {
		double rmsdThresholdSq = Math
				.pow(this.parameters.getRmsdThreshold(), 2);

		List<Point3d> centers = subunits.getOriginalCenters();
		List<Integer> seqClusterId = subunits.getClusterIds();

		List<Integer> permutations = new ArrayList<Integer>(centers.size());
		double[] dSqs = new double[centers.size()];
		boolean[] used = new boolean[centers.size()];
		Arrays.fill(used, false);

		for (int i = 0; i < centers.size(); i++) {
			Point3d tCenter = new Point3d(centers.get(i));
			transformation.transform(tCenter);
			int permutation = -1;
			double minDistSq = Double.MAX_VALUE;
			for (int j = 0; j < centers.size(); j++) {
				if (seqClusterId.get(i) == seqClusterId.get(j)) {
					if (!used[j]) {
						double dSq = tCenter.distanceSquared(centers.get(j));
						if (dSq < minDistSq && dSq <= rmsdThresholdSq) {
							minDistSq = dSq;
							permutation = j;
							dSqs[j] = dSq;
						}
					}
				}
			}
			// can't map to itself
			if (permutations.size() == permutation) {
				permutation = -1;
			}

			if (permutation != -1) {
				used[permutation] = true;
			}

			permutations.add(permutation);
		}

		return permutations;
	}
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
1045303
function getClosestPoints( simplex ){

        // see http://www.codezealot.org/archives/153
        // for algorithm details

        // we know that the position of the last point
        // is very close to the previous. (by nature of the distance test)
        // this won't give great results for the closest
        // points algorithm, so let's use the previous two
        var len = simplex.length
            ,last = simplex[ len - 2 ]
            ,prev = simplex[ len - 3 ]
            ,scratch = Physics.scratchpad()
            ,A = scratch.vector().clone( last.pt )
            // L = B - A
            ,L = scratch.vector().clone( prev.pt ).vsub( A )
            ,lambdaB
            ,lambdaA
            ;

        if ( L.equals(Physics.vector.zero) ){

            // oh.. it's a zero vector. So A and B are both the closest.
            // just use one of them
            return scratch.done({

                a: last.a,
                b: last.b
            });
        }

        lambdaB = - L.dot( A ) / L.normSq();
        lambdaA = 1 - lambdaB;

        if ( lambdaA <= 0 ){
            // woops.. that means the closest simplex point
            // isn't on the line it's point B itself
            return scratch.done({
                a: prev.a,
                b: prev.b
            });
        } else if ( lambdaB <= 0 ){
            // vice versa
            return scratch.done({
                a: last.a,
                b: last.b
            });
        }

        // guess we'd better do the math now...
        return scratch.done({
            // a closest = lambdaA * Aa + lambdaB * Ba
            a: A.clone( last.a ).mult( lambdaA ).vadd( L.clone( prev.a ).mult( lambdaB ) ).values(),
            // b closest = lambdaA * Ab + lambdaB * Bb
            b: A.clone( last.b ).mult( lambdaA ).vadd( L.clone( prev.b ).mult( lambdaB ) ).values()
        });
    }
1149227
function nelderMead(f, x0, parameters) {
        parameters = parameters || {};

        var maxIterations = parameters.maxIterations || x0.length * 200,
            nonZeroDelta = parameters.nonZeroDelta || 1.05,
            zeroDelta = parameters.zeroDelta || 0.001,
            minErrorDelta = parameters.minErrorDelta || 1e-6,
            minTolerance = parameters.minErrorDelta || 1e-5,
            rho = parameters.rho !== undefined ? parameters.rho : 1,
            chi = parameters.chi !== undefined ? parameters.chi : 2,
            psi = parameters.psi !== undefined ? parameters.psi : -0.5,
            sigma = parameters.sigma !== undefined ? parameters.sigma : 0.5,
            maxDiff;

        // initialize simplex.
        var N = x0.length,
            simplex = new Array(N + 1);
        simplex[0] = x0;
        simplex[0].fx = f(x0);
        simplex[0].id = 0;
        for (var i = 0; i < N; ++i) {
            var point = x0.slice();
            point[i] = point[i] ? point[i] * nonZeroDelta : zeroDelta;
            simplex[i + 1] = point;
            simplex[i + 1].fx = f(point);
            simplex[i + 1].id = i + 1;
        }

        function updateSimplex(value) {
            for (var i = 0; i < value.length; i++) {
                simplex[N][i] = value[i];
            }
            simplex[N].fx = value.fx;
        }

        var sortOrder = function sortOrder(a, b) {
            return a.fx - b.fx;
        };

        var centroid = x0.slice(),
            reflected = x0.slice(),
            contracted = x0.slice(),
            expanded = x0.slice();

        for (var iteration = 0; iteration < maxIterations; ++iteration) {
            simplex.sort(sortOrder);

            if (parameters.history) {
                // copy the simplex (since later iterations will mutate) and
                // sort it to have a consistent order between iterations
                var sortedSimplex = simplex.map(function (x) {
                    var state = x.slice();
                    state.fx = x.fx;
                    state.id = x.id;
                    return state;
                });
                sortedSimplex.sort(function (a, b) {
                    return a.id - b.id;
                });

                parameters.history.push({ x: simplex[0].slice(),
                    fx: simplex[0].fx,
                    simplex: sortedSimplex });
            }

            maxDiff = 0;
            for (i = 0; i < N; ++i) {
                maxDiff = Math.max(maxDiff, Math.abs(simplex[0][i] - simplex[1][i]));
            }

            if (Math.abs(simplex[0].fx - simplex[N].fx) < minErrorDelta && maxDiff < minTolerance) {
                break;
            }

            // compute the centroid of all but the worst point in the simplex
            for (i = 0; i < N; ++i) {
                centroid[i] = 0;
                for (var j = 0; j < N; ++j) {
                    centroid[i] += simplex[j][i];
                }
                centroid[i] /= N;
            }

            // reflect the worst point past the centroid  and compute loss at reflected
            // point
            var worst = simplex[N];
            weightedSum(reflected, 1 + rho, centroid, -rho, worst);
            reflected.fx = f(reflected);

            // if the reflected point is the best seen, then possibly expand
            if (reflected.fx < simplex[0].fx) {
                weightedSum(expanded, 1 + chi, centroid, -chi, worst);
                expanded.fx = f(expanded);
                if (expanded.fx < reflected.fx) {
                    updateSimplex(expanded);
                } else {
                    updateSimplex(reflected);
                }
            }

            // if the reflected point is worse than the second worst, we need to
            // contract
            else if (reflected.fx >= simplex[N - 1].fx) {
                    var shouldReduce = false;

                    if (reflected.fx > worst.fx) {
                        // do an inside contraction
                        weightedSum(contracted, 1 + psi, centroid, -psi, worst);
                        contracted.fx = f(contracted);
                        if (contracted.fx < worst.fx) {
                            updateSimplex(contracted);
                        } else {
                            shouldReduce = true;
                        }
                    } else {
                        // do an outside contraction
                        weightedSum(contracted, 1 - psi * rho, centroid, psi * rho, worst);
                        contracted.fx = f(contracted);
                        if (contracted.fx < reflected.fx) {
                            updateSimplex(contracted);
                        } else {
                            shouldReduce = true;
                        }
                    }

                    if (shouldReduce) {
                        // if we don't contract here, we're done
                        if (sigma >= 1) break;

                        // do a reduction
                        for (i = 1; i < simplex.length; ++i) {
                            weightedSum(simplex[i], 1 - sigma, simplex[0], sigma, simplex[i]);
                            simplex[i].fx = f(simplex[i]);
                        }
                    }
                } else {
                    updateSimplex(reflected);
                }
        }

        simplex.sort(sortOrder);
        return { fx: simplex[0].fx,
            x: simplex[0] };
    }
231654
def bilinear(px, py, band_array, gt):
    
    #import malib
    #band_array = malib.checkma(band_array)
    ndv = band_array.fill_value
    ny, nx = band_array.shape
    # Half raster cell widths
    hx = gt[1]/2.0
    hy = gt[5]/2.0
    # Calculate raster lower bound indices from point
    fx =(px -(gt[0] + hx))/gt[1]
    fy =(py -(gt[3] + hy))/gt[5]
    ix1 = int(np.floor(fx))
    iy1 = int(np.floor(fy))
    # Special case where point is on upper bounds
    if fx == float(nx - 1):
        ix1 -= 1
    if fy == float(ny - 1):
        iy1 -= 1
    # Upper bound indices on raster
    ix2 = ix1 + 1
    iy2 = iy1 + 1
    # Test array bounds to ensure point is within raster midpoints
    if(ix1 < 0) or(iy1 < 0) or(ix2 > nx - 1) or(iy2 > ny - 1):
        return ndv
    # Calculate differences from point to bounding raster midpoints
    dx1 = px -(gt[0] + ix1*gt[1] + hx)
    dy1 = py -(gt[3] + iy1*gt[5] + hy)
    dx2 =(gt[0] + ix2*gt[1] + hx) - px
    dy2 =(gt[3] + iy2*gt[5] + hy) - py
    # Use the differences to weigh the four raster values
    div = gt[1]*gt[5]
    return(band_array[iy1,ix1]*dx2*dy2/div +
            band_array[iy1,ix2]*dx1*dy2/div +
            band_array[iy2,ix1]*dx2*dy1/div +
            band_array[iy2,ix2]*dx1*dy1/div)
1218958
private void diffClique(float[] weights, int p) {
		for (int t = 0; t < orders.length; t++) {
			if (data[p][t] == -1)
				continue;

			if (orders[t] == 0) {
				if (golds[p] != preds[p]) {
					int ts = data[p][t] + golds[p];
					int ps = data[p][t] + preds[p];
					adjust(weights, ts, ps);
				}
			}
			if (p > 0 && orders[t] == 1) {
				int ts = data[p][t] + (golds[p - 1] * ysize + golds[p]);
				int ps = data[p][t] + (preds[p - 1] * ysize + preds[p]);
				adjust(weights, ts, ps);
			}
		}
	}
483798
function sqSegSegDist(x0, y0, x1, y1, x2, y2, x3, y3) {
    var ux = x1 - x0;
    var uy = y1 - y0;
    var vx = x3 - x2;
    var vy = y3 - y2;
    var wx = x0 - x2;
    var wy = y0 - y2;
    var a = ux * ux + uy * uy;
    var b = ux * vx + uy * vy;
    var c = vx * vx + vy * vy;
    var d = ux * wx + uy * wy;
    var e = vx * wx + vy * wy;
    var D = a * c - b * b;

    var sc, sN, tc, tN;
    var sD = D;
    var tD = D;

    if (D === 0) {
        sN = 0;
        sD = 1;
        tN = e;
        tD = c;
    } else {
        sN = b * e - c * d;
        tN = a * e - b * d;
        if (sN < 0) {
            sN = 0;
            tN = e;
            tD = c;
        } else if (sN > sD) {
            sN = sD;
            tN = e + b;
            tD = c;
        }
    }

    if (tN < 0.0) {
        tN = 0.0;
        if (-d < 0.0) sN = 0.0;
        else if (-d > a) sN = sD;
        else {
            sN = -d;
            sD = a;
        }
    } else if (tN > tD) {
        tN = tD;
        if ((-d + b) < 0.0) sN = 0;
        else if (-d + b > a) sN = sD;
        else {
            sN = -d + b;
            sD = a;
        }
    }

    sc = sN === 0 ? 0 : sN / sD;
    tc = tN === 0 ? 0 : tN / tD;

    var cx = (1 - sc) * x0 + sc * x1;
    var cy = (1 - sc) * y0 + sc * y1;
    var cx2 = (1 - tc) * x2 + tc * x3;
    var cy2 = (1 - tc) * y2 + tc * y3;
    var dx = cx2 - cx;
    var dy = cy2 - cy;

    return dx * dx + dy * dy;
}
231217
private static void gaussian(float a[][],
	                             int index[]) {
		int n = index.length;
		float c[] = new float[n];
		// Initialize the index
		for (int i = 0; i < n; ++i) {
			index[i] = i;
		}
		// Find the rescaling factors, one from each row
		for (int i = 0; i < n; ++i){
			float c1 = 0;
			for (int j = 0; j < n; ++j){
				float c0 = Math.abs(a[i][j]);
				if (c0 > c1) {
					c1 = c0;
				}
			}
			c[i] = c1;
		}
		// Search the pivoting element from each column
		int k = 0;
		for (int j = 0; j < n - 1; ++j){
			float pi1 = 0;
			for (int i = j; i < n; ++i){
				float pi0 = Math.abs(a[index[i]][j]);
				pi0 /= c[index[i]];
				if (pi0 > pi1) {

					pi1 = pi0;
					k = i;
				}
			}
			// Interchange rows according to the pivoting order
			int itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;
			for (int i = j + 1; i < n; ++i){
				float pj = a[index[i]][j] / a[index[j]][j];
				// Record pivoting ratios below the diagonal
				a[index[i]][j] = pj;

				// Modify other elements accordingly
				for (int l = j + 1; l < n; ++l){
					a[index[i]][l] -= pj * a[index[j]][l];
				}
			}
		}
	}
674938
private void diffDownClique(float[] weights, int[] orders, int p) {
		for (int t = 0; t < orders.length; t++) {
			if (data[1][p][t] == -1)
				continue;

			int base = data[1][p][t];
			if (orders[t] == 0) {
				if (golds[1][p] != preds[1][p]) {
					int ts = base + golds[1][p];
					int ps = base + preds[1][p];
					adjust(weights, ts, ps);
				}
			}
			if (orders[t] == -1) {
				if (golds[0][p] != preds[0][p] || golds[1][p] != preds[1][p]) {
					int ts = base + golds[0][p] * ysize[1] + golds[1][p];
					int ps = base + preds[0][p] * ysize[1] + preds[1][p];
					adjust(weights, ts, ps);
				}
			}
			if (p > 0) {
				if (orders[t] == 1) {
					if (golds[1][p - 1] != preds[1][p - 1]
							|| golds[1][p] != preds[1][p]) {
						int ts = base + golds[1][p - 1] * ysize[1] + golds[1][p];
						int ps = base + preds[1][p - 1] * ysize[1] + preds[1][p];
						adjust(weights, ts, ps);
					}
				}
				if (orders[t] == 2) {
					int ts = base + (golds[1][p - 1] * ysize[0] + golds[0][p])
							* ysize[1] + golds[1][p];
					int ps = base + (preds[1][p - 1] * ysize[0] + preds[0][p])
							* ysize[1] + preds[1][p];
					adjust(weights, ts, ps);
				}
			}
		}
	}
1026230
def noise(matrix):
    

    D = shortest_path_interpolation(matrix, strict=True)
    return np.random.poisson(lam=D)
291673
private void tred2()
    {
        for(int j = 0; j < n; j++)
            d[j] = V.get(n-1, j);

        // Householder reduction to tridiagonal form.

        for (int i = n - 1; i > 0; i--)
        {

            // Scale to avoid under/overflow.

            double scale = 0.0;
            double h = 0.0;
            for (int k = 0; k < i; k++)
            {
                scale = scale + abs(d[k]);
            }
            if (scale == 0.0)
            {
                e[i] = d[i - 1];
                for (int j = 0; j < i; j++)
                {
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                    V.set(j, i, 0.0);
                }
            }
            else
            {

                // Generate Householder vector.

                for (int k = 0; k < i; k++)
                {
                    d[k] /= scale;
                    h += d[k] * d[k];
                }
                double f = d[i - 1];
                double g = sqrt(h);
                if (f > 0)
                    g = -g;
                e[i] = scale * g;
                h -= f * g;
                d[i - 1] = f - g;
                Arrays.fill(e, 0, i, 0.0);

                // Apply similarity transformation to remaining columns.

                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    V.set(j, i, f);
                    g = e[j] + V.get(j, j) * f;
                    for (int k = j + 1; k <= i - 1; k++)
                    {
                        g += V.get(k,j) * d[k];
                        e[k] += V.get(k, j) * f;
                    }
                    e[j] = g;
                }
                f = 0.0;
                for (int j = 0; j < i; j++)
                {
                    e[j] /= h;
                    f += e[j] * d[j];
                }
                double hh = f / (h + h);
                for (int j = 0; j < i; j++)
                {
                    e[j] -= hh * d[j];
                }
                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    g = e[j];
                    
                    for (int k = j; k <= i - 1; k++)
                    {
                        V.increment(k, j, -(f * e[k] + g * d[k]));
                    }
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                }
            }
            d[i] = h;
        }


        // Accumulate transformations.

        for (int i = 0; i < n - 1; i++)
        {
            V.set(n-1, i, V.get(i, i));
            V.set(i, i, 1.0);
            double h = d[i + 1];
            if (h != 0.0)
            {
                for (int k = 0; k <= i; k++)
                {
                    d[k] = V.get(k, i+1) / h;
                }
                for (int j = 0; j <= i; j++)
                {
                    double g = 0.0;
                    for (int k = 0; k <= i; k++)
                    {
                        g += V.get(k, i+1) * V.get(k, j);
                    }
                    
                    RowColumnOps.addMultCol(V, j, 0, i+1, -g, d);
                }
            }

            RowColumnOps.fillCol(V, i+1, 0, i+1, 0.0);
            
        }
        for (int j = 0; j < n; j++)
        {
            d[j] = V.get(n-1, j);
            V.set(n-1, j, 0.0);
        }
        V.set(n-1, n-1, 1.0);
        e[0] = 0.0;
    }
413808
function closestPointnd(c, positions, x, result) {
  var D = numeric.rep([c.length, c.length], 0.0);
  var dvec = numeric.rep([c.length], 0.0);
  for(var i=0; i<c.length; ++i) {
    var pi = positions[c[i]];
    dvec[i] = numeric.dot(pi, x);
    for(var j=0; j<c.length; ++j) {
      var pj = positions[c[j]];
      D[i][j] = D[j][i] = numeric.dot(pi, pj);
    }
  }
  var A = numeric.rep([c.length, c.length+2], 0.0);
  var b = numeric.rep([c.length+2], 0.0);
  b[0] = 1.0-EPSILON;
  b[1] = -(1.0+EPSILON);
  for(var i=0; i<c.length; ++i) {
    A[i][0]   = 1;
    A[i][1]   = -1
    A[i][i+2] = 1;
  }
  for(var attempts=0; attempts<15; ++attempts) {
    var fortran_poop = numeric.solveQP(D, dvec, A, b);
    if(fortran_poop.message.length > 0) {
      //Quadratic form may be singular, perturb and resolve
      for(var i=0; i<c.length; ++i) {
        D[i][i] += 1e-8;
      }
      continue;
    } else if(isNaN(fortran_poop.value[0])) {
      break;
    } else {
      //Success!
      var solution = fortran_poop.solution;
      for(var i=0; i<x.length; ++i) {
        result[i] = 0.0;
        for(var j=0; j<solution.length; ++j) {
          result[i] += solution[j] * positions[c[j]][i];
        }
      }
      return 2.0 * fortran_poop.value[0] + numeric.dot(x,x);
    }
  }
  for(var i=0; i<x.length; ++i) {
    result[i] = Number.NaN;
  }
  return Number.NaN;
}
388266
def indexes_to_fix_for_low_rank(rank, size):
    
    #first we'll work out what to keep, and the do the set difference.

    #here are the indexes of the first column, which are the triangular numbers
    n = np.arange(size)
    triangulars = (n**2 + n) / 2
    keep = []
    for i in range(rank):
        keep.append(triangulars[i:] + i)
    #add the diagonal
    keep.append(triangulars[1:]-1)
    keep.append((size**2 + size)/2 -1)# the very last element
    keep = np.hstack(keep)

    return np.setdiff1d(np.arange((size**2+size)/2), keep)
488268
def polyfit2d(x, y, z, order=3 #bounds=None
              ):
    '''
     
    '''
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(list(range(order+1)), list(range(order+1)))
    for k, (i,j) in enumerate(ij):
        G[:,k] = x**i * y**j
    m = np.linalg.lstsq(G, z)[0]
    return m