Column reduction and reduction transfer steps from LAPJV algorithm
13407
def _column_reduction(self):
        
        #assign each column to its lowest cost row, ensuring that only row
        #or column is assigned once
        i1, j = np.unique(np.argmin(self.c, axis=0), return_index=True)
        self._x[i1] = j

        #if problem is solved, return
        if len(i1) == self.n:
            return False

        self._y[j] = i1

        #reduction_transfer
        #tempc is array with previously assigned matchings masked
        self._v = np.min(self.c, axis=0)
        tempc = self.c.copy()
        tempc[i1, j] = np.inf
        mu = np.min(tempc[i1, :] - self._v[None, :], axis=1)
        self._v[j] -= mu
        return True
241475
def d2Ibr_dV2(Ybr, V, lam):
    
    nb = len(V)
    diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = spdiag(mul(-(Ybr.T * lam), V))
    Hva = -1j * Haa * diaginvVm
    Hav = Hva
    Hvv = spmatrix([], [], [], (nb, nb))

    return Haa, Hav, Hva, Hvv
1246933
public static void tred2(int n, double v[][], double d[], double e[]) {

    //  This is derived from the Algol procedures tred2 by
    //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    //  Fortran subroutine in EISPACK.

    System.arraycopy(v[n - 1], 0, d, 0, n);

    // 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 + Math.abs(d[k]);
      }
      if (scale == 0.0) {
        e[i] = d[i - 1];
        for (int j = 0; j < i; j++) {
          d[j] = v[i - 1][j];
          v[i][j] = 0.0;
          v[j][i] = 0.0;
        }
      } else {

        h = householderIteration(i, scale, v, d, e);

      }
      d[i] = h;
    }

    // Accumulate transformations.
    accumulateTransformations(n, v, d);

    e[0] = 0.0;
  }
773514
private void dfs(Visitor visitor, int v, int[] cc, int id) {
        visitor.visit(v);
        cc[v] = id;
        for (int t = 0; t < n; t++) {
            if (graph[v][t] != 0.0) {
                if (cc[t] == -1) {
                    dfs(visitor, t, cc, id);
                }
            }
        }
    }
681440
def precompute_optimzation_Y(laplacian_matrix, n_samples, relaxation_kwds):
    """"""
    relaxation_kwds.setdefault('presave',False)
    relaxation_kwds.setdefault('presave_name','pre_comp_current.npy')
    relaxation_kwds.setdefault('verbose',False)
    if relaxation_kwds['verbose']:
        print ('Making Lk and nbhds')
    Lk_tensor, nbk, si_map = \
        compute_Lk(laplacian_matrix, n_samples, relaxation_kwds['subset'])
    if relaxation_kwds['presave']:
        raise NotImplementedError('Not yet implemented presave')
    return { 'Lk': Lk_tensor, 'nbk': nbk, 'si_map': si_map }
1097517
def subdivide(face, detail)
      cols = 2.0 ** detail
      a = prepare(@vertices[face.a])
      b = prepare(@vertices[face.b])
      c = prepare(@vertices[face.c])
      v = []

      # Construct all of the vertices for this subdivision.
      for i in 0..cols do
        v[i] = []

        aj = prepare(a.clone.lerp(c, i.to_f / cols.to_f))
        bj = prepare(b.clone.lerp(c, i.to_f / cols.to_f))
        rows = cols - i

        for j in 0..rows do
          v[i][j] = if j.zero? && i == cols
                      aj
                    else
                      prepare(aj.clone.lerp(bj, j.to_f / rows.to_f))
                    end
        end
      end

      # Construct all of the faces
      for i in 0...cols do
        for j in (0...(2 * (cols - i) - 1)) do
          k = j/2

          if j.even?
            make(
              v[i][k + 1],
              v[i + 1][k],
              v[i][k]
            )
          else
            make(
              v[i][k + 1],
              v[i + 1][k + 1],
              v[i + 1][k]
            )
          end
        end
      end
    end
1145532
def multidimensional_discretization(rho, sigma, N=3, method='rouwenhorst', m=2):
    

    # rho is assumed to be a scalar
    # sigma is a positive symmetric matrix
    # N number of points in each non-degenerate dimension
    # m : standard deviations to approximate

    import scipy.linalg
    from itertools import product

    d = sigma.shape[1]

    sigma = sigma.copy()

    zero_columns = np.where(sigma.sum(axis=0)==0)[0]
    for i in zero_columns:
        sigma[i,i] = 1

    L = scipy.linalg.cholesky(sigma)

    N = int(N)

    if method=='tauchen':
        [nodes_1d, probas_1d] = tauchen(N, 0, rho, 1, m=m)
    elif method=='rouwenhorst':
        [nodes_1d, probas_1d] = rouwenhorst(rho, 1, N)

    markov_nodes = np.array( list( product( *([nodes_1d]*d))) ).T
    markov_indices = np.array( list( product( *([range(N)]*d)  )  ) ).T

    markov_nodes = np.dot(L, markov_nodes)


    transition_matrix = 1
    for i in range(d):
        transition_matrix = np.kron(transition_matrix, probas_1d)

    markov_nodes = np.ascontiguousarray(markov_nodes.T)

    for i in zero_columns:
        markov_nodes[:,i] = 0

    return [markov_nodes, transition_matrix]
561948
function spfa(graph, s) {
  const distance = {};
  const previous = {};
  const queue = {};
  const isInQue = {};
  const cnt = {};
  let head = 0;
  let tail = 1;
  // initialize
  distance[s] = 0;
  queue[0] = s;
  isInQue[s] = true;
  cnt[s] = 1;
  graph.vertices.forEach(v => {
    if (v !== s) {
      distance[v] = Infinity;
      isInQue[v] = false;
      cnt[v] = 0;
    }
  });

  let currNode;
  while (head !== tail) {
    currNode = queue[head++];
    isInQue[currNode] = false;
    const neighbors = graph.neighbors(currNode);
    for (let i = 0; i < neighbors.length; i++) {
      const v = neighbors[i];
      // relaxation
      const newDistance = distance[currNode] + graph.edge(currNode, v);
      if (newDistance < distance[v]) {
        distance[v] = newDistance;
        previous[v] = currNode;
        if (!isInQue[v]) {
          queue[tail++] = v;
          isInQue[v] = true;
          cnt[v]++;
          if (cnt[v] > graph.vertices.size) {
            // indicates negative-weighted cycle
            return {
              distance: {}
            };
          }
        }
      }
    }
  }

  return {
    distance,
    previous
  };
}
523198
public static void jacobian_Control4( DMatrixRMaj L_full ,
										  double beta[] , DMatrixRMaj A )
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2]; double b3 = beta[3];

		final double ld[] = L_full.data;

		for( int i = 0; i < 6; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2 +   ld[li+3]*b3;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+4]*b1 +   ld[li+5]*b2 +   ld[li+6]*b3;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+5]*b1 + 2*ld[li+7]*b2 +   ld[li+8]*b3;
			A.data[indexA++] =   ld[li+3]*b0 +   ld[li+6]*b1 +   ld[li+8]*b2 + 2*ld[li+9]*b3;
		}
	}
500618
function svdJs() {
                var A = this;
                var V = Matrix$b.I(A.rows());
                var S = A.transpose();
                var U = Matrix$b.I(A.cols());
                var err = Number.MAX_VALUE;
                var i = 0;
                var maxLoop = 100;

                while(err > 2.2737e-13 && i < maxLoop) {
                    var qr = S.transpose().qrJs();
                    S = qr.R;
                    V = V.x(qr.Q);
                    qr = S.transpose().qrJs();
                    U = U.x(qr.Q);
                    S = qr.R;

                    var e = S.triu(1).unroll().norm();
                    var f = S.diagonal().norm();

                    if(f == 0)
                        f = 1;

                    err = e / f;

                    i++;
                }

                var ss = S.diagonal();
                var s = [];

                for(var i = 1; i <= ss.cols(); i++) {
                    var ssn = ss.e(i);
                    s.push(Math.abs(ssn));

                    if(ssn < 0) {
                        for(var j = 0; j < U.rows(); j++) {
                            V.elements[j][i - 1] = -(V.elements[j][i - 1]);
                        }
                    }
                }

                return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
            }
403717
function d3_geom_hullCCW(i1, i2, i3, v) {
  var t, a, b, c, d, e, f;
  t = v[i1]; a = t[0]; b = t[1];
  t = v[i2]; c = t[0]; d = t[1];
  t = v[i3]; e = t[0]; f = t[1];
  return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
}
562076
private static int getDataRows(Chunk[] chunks, Frame f, int[] chunksIds, int cols) {
        double totalRows = 0;
        if(null != chunks) {
            for (Chunk ch : chunks) {
                totalRows += ch.len();
            }
        } else {
            for(int chunkId : chunksIds) {
                totalRows += f.anyVec().chunkLen(chunkId);
            }
        }
        return (int) Math.ceil(totalRows * cols / SPARSE_MATRIX_DIM);
    }
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }
891587
def __computeDecomp(self):
        
        primeNumbers = [getPrimeFactors(d) for d in self.globalDims]

        ns = [len(pns) for pns in primeNumbers]
        validDecomps = []
        self.validProcs = []
        for it in MultiArrayIter(ns):
            inds = it.getIndices()
            decomp = [primeNumbers[d][inds[d]] for d in range(self.ndims)]
            self.validProcs.append(reduce(operator.mul, decomp, 1))
            if reduce(operator.mul, decomp, 1) == self.nprocs:
                validDecomps.append(decomp)

        # sort and remove duplicates
        self.validProcs.sort()
        vprocs = []
        for vp in self.validProcs:
            if len(vprocs) == 0 or (len(vprocs) >= 1 and vp != vprocs[-1]):
                vprocs.append(vp)
        self.validProcs = vprocs

        if len(validDecomps) == 0:
            # no solution
            return

        # find the optimal decomp among all valid decomps
        minCost = float('inf')
        bestDecomp = validDecomps[0]
        for decomp in validDecomps:
            sizes = [self.globalDims[d]//decomp[d] for d in range(self.ndims)]
            volume = reduce(operator.mul, sizes, 1)
            surface = 0
            for d in range(self.ndims):
                surface += 2*reduce(operator.mul, sizes[:d], 1) * \
                    reduce(operator.mul, sizes[d+1:], 1)
            cost = surface / float(volume)
            if cost < minCost:
                bestDecomp = decomp
                minCost = cost
        self.decomp = bestDecomp

        # ok, we have a valid decomp, now build the sub-domain iterator
        self.mit = MultiArrayIter(self.decomp, rowMajor=self.rowMajor)

        # fill in the proc to index set map
        procId = 0
        self.proc2IndexSet = {}
        numCellsPerProc = [self.globalDims[d]//self.decomp[d]
                           for d in range(self.ndims)]
        for it in self.mit:
            nps = it.getIndices()
            self.proc2IndexSet[procId] = []
            for d in range(self.ndims):
                sbeg = nps[d]*numCellsPerProc[d]
                send = (nps[d] + 1)*numCellsPerProc[d]
                self.proc2IndexSet[procId].append(slice(sbeg, send))
            procId += 1
750377
function constrainedMDSLayout(areas, params) {
        params = params || {};
        var restarts = params.restarts || 10;

        // bidirectionally map sets to a rowid  (so we can create a matrix)
        var sets = [],
            setids = {},
            i;
        for (i = 0; i < areas.length; ++i) {
            var area = areas[i];
            if (area.sets.length == 1) {
                setids[area.sets[0]] = sets.length;
                sets.push(area);
            }
        }

        var matrices = getDistanceMatrices(areas, sets, setids),
            distances = matrices.distances,
            constraints = matrices.constraints;

        // keep distances bounded, things get messed up otherwise.
        // TODO: proper preconditioner?
        var norm = norm2(distances.map(norm2)) / distances.length;
        distances = distances.map(function (row) {
            return row.map(function (value) {
                return value / norm;
            });
        });

        var obj = function obj(x, fxprime) {
            return constrainedMDSGradient(x, fxprime, distances, constraints);
        };

        var best, current;
        for (i = 0; i < restarts; ++i) {
            var initial = zeros(distances.length * 2).map(Math.random);

            current = conjugateGradient(obj, initial, params);
            if (!best || current.fx < best.fx) {
                best = current;
            }
        }
        var positions = best.x;

        // translate rows back to (x,y,radius) coordinates
        var circles = {};
        for (i = 0; i < sets.length; ++i) {
            var set = sets[i];
            circles[set.sets[0]] = {
                x: positions[2 * i] * norm,
                y: positions[2 * i + 1] * norm,
                radius: Math.sqrt(set.size / Math.PI)
            };
        }

        if (params.history) {
            for (i = 0; i < params.history.length; ++i) {
                scale(params.history[i].x, norm);
            }
        }
        return circles;
    }
1009588
private void dfs(int v, int[] cc, int id) {
        cc[v] = id;
        for (int t = 0; t < n; t++) {
            if (graph[v][t] != 0.0) {
                if (cc[t] == -1) {
                    dfs(t, cc, id);
                }
            }
        }
    }
449886
private void intersectionMatrix(int[][] res, ArrayList<int[]> partsx, ArrayList<int[]> partsy, int gridsize) {
    for(int x = 0; x < gridsize; x++) {
      final int[] px = partsx.get(x);
      final int[] rowx = res[x];
      for(int y = 0; y < gridsize; y++) {
        int[] py = partsy.get(y);
        rowx[y] = intersectionSize(px, py);
      }
    }
  }