Maximum profit perfect matching

    for minimum cost perfect matching just inverse the weights

    :param G: squared weight matrix of a complete bipartite graph
    :complexity: :math:`O(n^4)`
775712
def kuhn_munkres(G):      # maximum profit bipartite matching in O(n^4)
    
    assert len(G) == len(G[0])
    n = len(G)
    mu = [None] * n                 # Empty matching
    mv = [None] * n
    lu = [max(row) for row in G]    # Trivial labels
    lv = [0] * n
    for u0 in range(n):
        if mu[u0] is None:          # Free node
            while True:
                au = [False] * n    # Empty alternating tree
                av = [False] * n
                if improve_matching(G, u0, mu, mv, au, av, lu, lv):
                    break
                improve_labels(G, au, av, lu, lv)
    return (mu,  sum(lu) + sum(lv))
833669
def prim(G, start, weight='weight'):
    
    if len(connected_components(G)) != 1:
        raise GraphInsertError("Prim algorithm work with connected graph only")
    if start not in G.vertices:
        raise GraphInsertError("Vertex %s doesn't exist." % (start,))
    pred = {}
    key = {}
    pqueue = {}
    lowest = 0
    for edge in G.edges:
        if G.edges[edge][weight] > lowest:
            lowest = G.edges[edge][weight]
    for vertex in G.vertices:
        pred[vertex] = None
        key[vertex] = 2 * lowest
    key[start] = 0
    for vertex in G.vertices:
        pqueue[vertex] = key[vertex]
    while pqueue:
        current = popmin(pqueue, lowest)
        for neighbor in G.vertices[current]:
            if (neighbor in pqueue and
                    G.edges[(current, neighbor)][weight] < key[neighbor]):
                pred[neighbor] = current
                key[neighbor] = G.edges[(current, neighbor)][weight]
                pqueue[neighbor] = G.edges[(current, neighbor)][weight]
    return pred
823928
def from_edge_pairs(pairs, num_vertices=None, symmetric=False, weights=None):
  '''
  '''
  if not symmetric:
    if weights is None:
      return EdgePairGraph(pairs, num_vertices=num_vertices)
    row, col = np.asarray(pairs).T
    row, weights = np.broadcast_arrays(row, weights)
    shape = None if num_vertices is None else (num_vertices, num_vertices)
    adj = ss.coo_matrix((weights, (row, col)), shape=shape)
    return SparseAdjacencyMatrixGraph(adj)
  # symmetric case
  G = SymmEdgePairGraph(pairs, num_vertices=num_vertices)
  if weights is None:
    return G
  # Convert to sparse adj graph with provided edge weights
  s = G.matrix('coo').astype(float)
  # shenanigans to assign edge weights in the right order
  flat_idx = np.ravel_multi_index(s.nonzero(), s.shape)
  r, c = np.transpose(pairs)
  rc_idx = np.ravel_multi_index((r,c), s.shape)
  cr_idx = np.ravel_multi_index((c,r), s.shape)
  order = np.argsort(flat_idx)
  flat_idx = flat_idx[order]
  s.data[order[np.searchsorted(flat_idx, rc_idx)]] = weights
  s.data[order[np.searchsorted(flat_idx, cr_idx)]] = weights
  return SparseAdjacencyMatrixGraph(s)
1071560
def tsp_gurobi(edges):
    
    from gurobipy import Model, GRB, quicksum

    edges = populate_edge_weights(edges)
    incoming, outgoing, nodes = node_to_edge(edges)
    idx = dict((n, i) for i, n in enumerate(nodes))
    nedges = len(edges)
    n = len(nodes)

    m = Model()

    def step(x): return "u_{0}".format(x)
    # Create variables
    vars = {}
    for i, (a, b, w) in enumerate(edges):
        vars[i] = m.addVar(obj=w, vtype=GRB.BINARY, name=str(i))
    for u in nodes[1:]:
        u = step(u)
        vars[u] = m.addVar(obj=0, vtype=GRB.INTEGER, name=u)
    m.update()

    # Bounds for step variables
    for u in nodes[1:]:
        u = step(u)
        vars[u].lb = 1
        vars[u].ub = n - 1

    # Add degree constraint
    for v in nodes:
        incoming_edges = incoming[v]
        outgoing_edges = outgoing[v]
        m.addConstr(quicksum(vars[x] for x in incoming_edges) == 1)
        m.addConstr(quicksum(vars[x] for x in outgoing_edges) == 1)

    # Subtour elimination
    edge_store = dict(((idx[a], idx[b]), i)
                      for i, (a, b, w) in enumerate(edges))

    # Given a list of edges, finds the shortest subtour
    def subtour(s_edges):
        visited = [False] * n
        cycles = []
        lengths = []
        selected = [[] for i in range(n)]
        for x, y in s_edges:
            selected[x].append(y)
        while True:
            current = visited.index(False)
            thiscycle = [current]
            while True:
                visited[current] = True
                neighbors = [x for x in selected[current] if not visited[x]]
                if len(neighbors) == 0:
                    break
                current = neighbors[0]
                thiscycle.append(current)
            cycles.append(thiscycle)
            lengths.append(len(thiscycle))
            if sum(lengths) == n:
                break
        return cycles[lengths.index(min(lengths))]

    def subtourelim(model, where):
        if where != GRB.callback.MIPSOL:
            return
        selected = []
        # make a list of edges selected in the solution
        sol = model.cbGetSolution([model._vars[i] for i in range(nedges)])
        selected = [edges[i] for i, x in enumerate(sol) if x > .5]
        selected = [(idx[a], idx[b]) for a, b, w in selected]
        # find the shortest cycle in the selected edge list
        tour = subtour(selected)
        if len(tour) == n:
            return
        # add a subtour elimination constraint
        c = tour
        incident = [edge_store[a, b] for a, b in pairwise(c + [c[0]])]
        model.cbLazy(quicksum(model._vars[x]
                              for x in incident) <= len(tour) - 1)

    m.update()

    m._vars = vars
    m.params.LazyConstraints = 1
    m.optimize(subtourelim)

    selected = [v.varName for v in m.getVars() if v.x > .5]
    selected = [int(x) for x in selected if x[:2] != "u_"]
    results = sorted(x for i, x in enumerate(edges) if i in selected) \
        if selected else None
    return results
259489
def leastsq_NxN(x, y, fit_offset=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
        x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()
    else:
        weights = None

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        xx_ = prod_sum_obs(x, x)
        xy_ = prod_sum_obs(x, y)

        if fit_offset:
            n_obs = x.shape[0] if weights is None else sum_obs(weights)
            x_ = sum_obs(x) / n_obs
            y_ = sum_obs(y) / n_obs
            gamma = (xy_ / n_obs - x_ * y_) / (xx_ / n_obs - x_ ** 2)
            offset = y_ - gamma * x_

            # fix negative offsets:
            idx = offset < 0
            gamma[idx] = xy_[idx] / xx_[idx]
            offset = np.clip(offset, 0, None)
        else:
            gamma = xy_ / xx_
            offset = np.zeros(x.shape[1])
    offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
    return offset, gamma
1156004
def _get_upper_bound(self):
        """

        if self.lap_type == 'normalized':
            return 2  # Equal iff the graph is bipartite.
        elif self.lap_type == 'combinatorial':
            bounds = []
            # Equal for full graphs.
            bounds += [self.n_vertices * np.max(self.W)]
            # Gershgorin circle theorem. Equal for regular bipartite graphs.
            # Special case of the below bound.
            bounds += [2 * np.max(self.dw)]
            # Anderson, Morley, Eigenvalues of the Laplacian of a graph.
            # Equal for regular bipartite graphs.
            if self.n_edges > 0:
                sources, targets, _ = self.get_edge_list()
                bounds += [np.max(self.dw[sources] + self.dw[targets])]
            # Merris, A note on Laplacian graph eigenvalues.
            if not self.is_directed():
                W = self.W
            else:
                W = utils.symmetrize(self.W, method='average')
            m = W.dot(self.dw) / self.dw  # Mean degree of adjacent vertices.
            bounds += [np.max(self.dw + m)]
            # Good review: On upper bounds for Laplacian graph eigenvalues.
            return min(bounds)
        else:
            raise ValueError('Unknown Laplacian type '
                             '{}'.format(self.lap_type))
220108
def g(number_of_links, size_x, size_y)
      a = (size_x + size_y) ** p
      b = size_x ** p
      c = size_x ** p

      number_of_links / (a - b - c)
    end
106207
function normalizeRanks(g) {
  var min = _.min(_.map(g.nodes(), function(v) { return g.node(v).rank; }));
  _.forEach(g.nodes(), function(v) {
    var node = g.node(v);
    if (_.has(node, "rank")) {
      node.rank -= min;
    }
  });
}
1129710
def optimize(self, G, params0=None, n_times=10, verbose=False, vmax=5, perturb=1e-3, factr=1e7):
        
        # set params0 from null if params0 is None
        if params0 is None:
            if self.null is None:
                if verbose:     print(".. fitting null model")
                self.fitNull()
            if self.bgRE:
                params0 = sp.concatenate([self.null['params0_g'], self.null['params0_n']])
            else:
                params0 = self.null['params0_n']
            params_was_None = True
        else:
            params_was_None = False
        G *= sp.sqrt(self.N/(G**2).sum())
        self._gp.covar.G = G
        start = TIME.time()
        for i in range(n_times):
            if params_was_None:
                n_params = self.Cr.getNumberParams()
                _params0 = {'covar': sp.concatenate([1e-3*sp.randn(n_params), params0])}
            else:
                _params0 = {'covar': params0}
            self._gp.setParams(_params0)
            conv, info = self._gp.optimize(factr=factr, verbose=verbose)
            conv *= self.Cr.K().diagonal().max()<vmax
            conv *= self.getLMLgrad() < 0.1
            if conv or not params_was_None: break
        self.infoOpt = info
        if not conv:
            warnings.warn("not converged")
        # return value
        RV = {}
        if self.P>1:
            RV['Cr']  = self.Cr.K()
            if self.bgRE: RV['Cg']  = self.Cg.K()
            RV['Cn']  = self.Cn.K()
        RV['time']  = sp.array([TIME.time()-start])
        RV['params0'] = _params0
        RV['nit'] = sp.array([info['nit']])
        RV['funcalls'] = sp.array([info['funcalls']])
        RV['var']    = self.getVariances()
        RV['conv']  = sp.array([conv])
        RV['NLLAlt']  = sp.array([self.getNLLAlt()])
        RV['LLR']    = sp.array([self.getLLR()])
        RV['LMLgrad'] = sp.array([self.getLMLgrad()])
        return RV
310864
def disjoint_mst(X, num_spanning_trees=3, metric='euclidean'):
  ''''''
  D = pairwise_distances(X, metric=metric)
  if metric == 'precomputed':
    D = D.copy()
  mst = minimum_spanning_tree(D)
  W = mst.copy()
  for i in range(1, num_spanning_trees):
    ii,jj = mst.nonzero()
    D[ii,jj] = np.inf
    D[jj,ii] = np.inf
    mst = minimum_spanning_tree(D)
    W = W + mst
  # MSTs are all one-sided, so we symmetrize here
  return Graph.from_adj_matrix(W + W.T)
1233595
def solve_full(z, Fval, DPhival, G, A):    
    M, N=G.shape
    P, N=A.shape

    """"""
    m=M    

    """Primal variable"""
    x=z[0:N]

    """Multiplier for equality constraints"""
    nu=z[N:N+P]

    """Multiplier for inequality constraints"""
    l=z[N+P:N+P+M]

    """Slacks"""
    s=z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = np.diag(l/s)

    """Condensed system"""
    if issparse(DPhival):
        if not issparse(A):
            A = csr_matrix(A)        
        H = DPhival + mydot(G.T, mydot(SIG, G))
        J = bmat([[H, A.T], [A, None]])
    else:
        if issparse(A):
            A = A.toarray()
        J = np.zeros((N+P, N+P))
        J[0:N, 0:N] = DPhival + mydot(G.T, mydot(SIG, G))            
        J[0:N, N:] = A.T
        J[N:, 0:N] = A

    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))

    """Prepare iterative solve via MINRES"""
    sign = np.zeros(N+P)
    sign[0:N/2] = 1.0
    sign[N/2:] = -1.0
    S = diags(sign, 0)
    J_new = mydot(S, csr_matrix(J))
    b_new = mydot(S, b)

    dJ_new = np.abs(J_new.diagonal())
    dPc = np.ones(J_new.shape[0])
    ind = (dJ_new > 0.0)
    dPc[ind] = 1.0/dJ_new[ind]
    Pc = diags(dPc, 0)    
    dxnu, info = minres(J_new, b_new, tol=1e-8, M=Pc)
    
    # dxnu = solve(J, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
262101
function simplify(g) {
  var simplified = new Graph().setGraph(g.graph());
  _.forEach(g.nodes(), function(v) { simplified.setNode(v, g.node(v)); });
  _.forEach(g.edges(), function(e) {
    var simpleLabel = simplified.edge(e.v, e.w) || { weight: 0, minlen: 1 };
    var label = g.edge(e);
    simplified.setEdge(e.v, e.w, {
      weight: simpleLabel.weight + label.weight,
      minlen: Math.max(simpleLabel.minlen, label.minlen)
    });
  });
  return simplified;
}
157619
public void getMaxBranching(CompleteGraph graph, int root,
                                int[] invertedMaxBranching) {
        // Number of nodes.
        int numNodes = graph.getNumberOfNodes();

        // Branching.
        SparseGraph maxBranching = new SparseGraph();

        /*
         * Weighted edges used to allow weight modification and avoid any impact
         * in the original graph.
         */
        SimpleWeightedEdge[][] edges = new SimpleWeightedEdge[numNodes][numNodes];

        // Disjoint sets for the strongly connected components (SCC).
        DisjointSets sPartition = new DisjointSets(numNodes);
        // Disjoint sets for the weakly connected components (WCC).
        DisjointSets wPartition = new DisjointSets(numNodes);

        /*
         * Priority queues for each strongly connected component. In the
         * beginning, each SCC is composed by exactly one node.
         */
        ArrayList<PriorityQueue<SimpleWeightedEdge>> incomingEdges = new ArrayList<PriorityQueue<SimpleWeightedEdge>>(
                                                                                                                      numNodes);

        // Unique incoming edges for each SCC. Initially, no SCC has such edge.
        SimpleWeightedEdge[] enter = new SimpleWeightedEdge[numNodes];

        /*
         * List of root components, i.e., SCCs that have no incoming edges
         * (enter[scc] == null). In the beginning, every SCC is a root
         * component.
         */
        LinkedList<Integer> rootComponents = new LinkedList<Integer>();

        // Root node of each root component.
        int[] min = new int[numNodes];

        for (int node = 0; node < numNodes; ++node) {

            // Every SCC is a root component.
            rootComponents.add(node);

            /*
             * The head of its root component is its only node. This array is
             * called min in Tarjan's paper.
             */
            min[node] = node;

            // Create a priority queue for each SCC.
            PriorityQueue<SimpleWeightedEdge> sccPriorityQueue = new PriorityQueue<SimpleWeightedEdge>();
            incomingEdges.add(sccPriorityQueue);

            // No incoming edge is considered (nor created) for the root node.
            if (node != root) {
                /*
                 * Create and add all incoming edges of <code>node</code> to its
                 * SCC priority queue.
                 */
                for (int from = 0; from < numNodes; ++from) {
                    if (from == node)
                        // Skip autocycle edges.
                        continue;
                    // Create an weighted edge and add it to the priority queue.
                    edges[from][node] = new SimpleWeightedEdge(from, node,
                                                               graph.getEdgeWeight(from, node));
                    sccPriorityQueue.add(edges[from][node]);
                }
            }
        }

        // Root component with no available incoming edges.
        LinkedList<Integer> doneRootComponents = new LinkedList<Integer>();

        while (!rootComponents.isEmpty()) {
            // Get some arbitrary root component.
            int sccTo = rootComponents.pop();
            // Maximum edge entering the root component 'sccTo'.
            SimpleWeightedEdge maxInEdge = incomingEdges.get(sccTo).poll();

            if (maxInEdge == null) {
                // No edge left to consider in this component.
                doneRootComponents.add(sccTo);
                continue;
            }

            // SCC component of edge 'e' from node: e = (from, to).
            int sccFrom = sPartition.find(maxInEdge.from);

            if (sccFrom == sccTo) {
                // Skip, for now, this component.
                rootComponents.add(sccTo);
                continue;
            }

            // Include the selected edge in the current branching.
            maxBranching.addEdge(maxInEdge.from, maxInEdge.to);

            // SCC component of edge 'e' from node, where e = (from, to).
            int wssFrom = wPartition.find(maxInEdge.from);
            // SCC component of edge 'e' to node, where e = (from, to).
            int wssTo = wPartition.find(maxInEdge.to);

            // Edge connects two weakly connected components.
            if (wssFrom != wssTo) {
                wPartition.union(wssFrom, wssTo);
                enter[sccTo] = maxInEdge;
                continue;
            }

            /*
             * Edge is within the same WCC, thus it inclusion will create a new
             * SCC by uniting some old SCCs (the ones on the path from e.to to
             * e.from).
             */
            double minEdgeWeight = Double.POSITIVE_INFINITY;
            int minScc = -1;
            SimpleWeightedEdge tmpEdge = maxInEdge;
            while (tmpEdge != null) {
                if (tmpEdge.weight < minEdgeWeight) {
                    minEdgeWeight = tmpEdge.weight;
                    minScc = sPartition.find(tmpEdge.to);
                }

                tmpEdge = enter[sPartition.find(tmpEdge.from)];
            }

            // Increment incoming edges weight.
            double inc = minEdgeWeight - maxInEdge.weight;
            for (SimpleWeightedEdge e : incomingEdges.get(sccTo))
                e.weight += inc;

            // Set the head of the current SCC.
            min[sccTo] = min[minScc];

            // Include all used SCCs in the current SCC.
            tmpEdge = enter[sccFrom];
            while (tmpEdge != null) {
                /*
                 * Increment incoming edges weight and include them in the
                 * current SCC priority queue.
                 */
                int tmpSccTo = sPartition.find(tmpEdge.to);
                inc = minEdgeWeight - tmpEdge.weight;
                for (SimpleWeightedEdge e : incomingEdges.get(tmpSccTo)) {
                    e.weight += inc;
                    incomingEdges.get(sccTo).add(e);
                }
                // Remove the priority queue of this SCC.
                incomingEdges.set(tmpSccTo, null);
                sPartition.union(sccTo, tmpSccTo);

                // Next edge.
                tmpEdge = enter[sPartition.find(tmpEdge.from)];
            }

            // Include the new SCC to be considered in the future.
            rootComponents.add(sccTo);
        }

        // Invert the maximum branching.
        boolean[] visited = new boolean[numNodes];
        for (int scc : doneRootComponents)
            invertMaximumBranching(min[scc], maxBranching, visited,
                                   invertedMaxBranching);
    }

    /**
     * Walk through the given branching from <code>node</code> and store the
     * inverted branching in <code>invertedMaxBranching</code>.
     * 
     * In fact, the given branching can include cycles. But it is only necessary
     * to disconsider the last edges of each cycle to get the real branching.
     * Thus, we use the array <code>visited</code>.
     * 
     * @param node
     * @param branching
     * @param visited
     * @param invertedMaxBranching
     */
    private void invertMaximumBranching(int node, SparseGraph branching,
                                        boolean[] visited, int[] invertedMaxBranching) {
        visited[node] = true;
        Set<Integer> toNodes = branching.getOutEdges(node);
        if (toNodes == null)
            return;
        for (int to : toNodes) {
            if (visited[to])
                continue;
            invertedMaxBranching[to] = node;
            invertMaximumBranching(to, branching, visited, invertedMaxBranching);
        }
    }
741168
public <E extends Edge> double[] compute(Graph<E> g) {
        // Perform a quick test for whether the vertices of g are a contiguous
        // sequence starting at 0, which makes the vertex mapping trivial
        if (!hasContiguousVertices(g))
            throw new IllegalArgumentException(
                "Vertices must be in continugous order");
        
        double[] centralities = new double[g.order()];
        IntIterator vertexIter = g.vertices().iterator();
        while (vertexIter.hasNext()) {
            int s = vertexIter.nextInt();
            Deque<Integer> S = new ArrayDeque<Integer>();

            // Initialize P to an empty list for each vertex
            List<List<Integer>> P = new ArrayList<List<Integer>>(g.order());
            for (int i = 0; i < g.order(); ++i)
                P.add(new ArrayList<Integer>());

            double[] sigma = new double[g.order()];
            sigma[s] = 1;

            double[] d = new double[g.order()];
            Arrays.fill(d, -1);
            d[s] = 0;
            
            Queue<Integer> Q = new ArrayDeque<Integer>();
            Q.add(s);
            while (!Q.isEmpty()) {
                int v = Q.poll();
                S.offer(v);
                IntIterator neighborIter = g.getNeighbors(v).iterator();
                while (neighborIter.hasNext()) { 
                    int w = neighborIter.nextInt();
                    // Check whether this is the first time we've seen vertex w
                    if (d[w] < 0) {
                        Q.offer(w);
                        d[w] = d[v] + 1;
                    }
                    // Check whether the shortest path to w is through v
                    if (d[w] == d[v] + 1) {
                        sigma[w] += sigma[v];
                        P.get(w).add(v);
                    }
                }
            }
            double[] delta = new double[g.order()];
            // S as a stack returns vertices in order of their non-increasing
            // distance from vertex s
            while (!S.isEmpty()) {
                int w = S.pollLast(); // get the top of the stack
                for (int v : P.get(w)) {
                    delta[v] += (sigma[v] / sigma[w]) * (1 + delta[w]);
                }
                if (w != s) {
                    centralities[w] += delta[w];
                }            
            }
        }
        return centralities;
    }
64027
private static MaxMarginalSet cliqueTreeToMaxMarginalSet(CliqueTree cliqueTree,
      FactorGraph originalFactorGraph) {
    for (int i = 0; i < cliqueTree.numFactors(); i++) {
      computeMarginal(cliqueTree, i, false);
    }
    return new FactorMaxMarginalSet(cliqueTree, originalFactorGraph.getConditionedValues());
  }
574812
def comparison_graph_py(arr1, arr2):
    
    # timeout is not implemented
    u1, v1, c1 = zip(*arr1)
    u2, v2, c2 = zip(*arr2)
    c1 = np.array(c1, dtype=int)
    c2 = np.array(c2, dtype=int)
    product = nx.Graph()
    c1 = c1[:, np.newaxis]  # transpose
    if NUMEXPR_AVAILABLE:
        m = ne.evaluate("c2 == c1")
    else:
        m = c2 == c1
    edges = []
    for x, y in zip(*np.nonzero(m)):
        edges.append({"u1": u1[x], "v1": v1[x], "u2": u2[y], "v2": v2[y]})
    # Graph.add_edges is expensive. Add adjacency dict manually.
    node = {}
    for e in edges:
        node[(e["u1"], e["u2"])] = {}
        node[(e["v1"], e["v2"])] = {}
    adj = node.copy()
    for e in edges:
        adj[(e["u1"], e["u2"])][(e["v1"], e["v2"])] = {}
        adj[(e["v1"], e["v2"])][(e["u1"], e["u2"])] = {}
    product = nx.Graph()
    product.node = node
    product.adj = adj
    return product
73272
def bound_bboxes(bboxes):
    
    group_x0 = min(map(lambda l: l[x0], bboxes))
    group_y0 = min(map(lambda l: l[y0], bboxes))
    group_x1 = max(map(lambda l: l[x1], bboxes))
    group_y1 = max(map(lambda l: l[y1], bboxes))
    return (group_x0, group_y0, group_x1, group_y1)