google-research

1
# coding=utf-8
2
# Copyright 2024 The Google Research Authors.
3
#
4
# Licensed under the Apache License, Version 2.0 (the "License");
5
# you may not use this file except in compliance with the License.
6
# You may obtain a copy of the License at
7
#
8
#     http://www.apache.org/licenses/LICENSE-2.0
9
#
10
# Unless required by applicable law or agreed to in writing, software
11
# distributed under the License is distributed on an "AS IS" BASIS,
12
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13
# See the License for the specific language governing permissions and
14
# limitations under the License.
15

16
"""This file contains the warm-start Ford-Fulkerson algorithm implmentation, the image sequence experiments
17

18
and the subroutines used by the experiment.
19
"""
20

21
import numpy as np
22
from queue import Queue
23
from collections import defaultdict
24
from copy import deepcopy
25
import time
26
from imagesegmentation import parseArgs, imageSegmentation
27
from augmentingPath import augmentingPath
28
import os
29
"""
30
One iteration of feasibility projection.
31
Algorithm: given a node u, BFS until find an excess-deficit path in residual graph
32
Prioritize paths that do not contain s or t, if no such path exists, allow using s or t.
33
Such a path must be found.
34
"""
35

36

37
def FeasRestoreIter(p, rGraph, excesses, s, t):
38

39
  def FindPath(chain, node, excesses, s, t, out_dir=True):
40
    # here into_dir = True denotes the path should carry flow OUT OF v.
41
    v, path = node, [node]
42
    path_flow = abs(excesses[v]) if v not in [s, t] else float('inf')
43
    u = chain[v]
44
    while u >= 0:
45
      path.append(u)
46
      (head, tail) = (u, v) if out_dir else (v, u)
47
      path_flow = min(rGraph[head][tail], path_flow)
48
      v = u
49
      u = chain[v]
50
    if v not in [s, t]:
51
      path_flow = min(path_flow, abs(excesses[v]))
52
    if out_dir:
53
      path.reverse()
54
    return path, path_flow
55

56
  if excesses[p] == 0:
57
    return [p], 0
58
  direc = (excesses[p] > 0)
59
  chain = defaultdict(int)
60
  q = Queue()
61
  q.put(p)
62
  visited = {p}
63
  chain[p] = -1
64

65
  while not q.empty():
66
    u = q.get()
67

68
    for v in list(rGraph[u].keys()) + [s, t]:
69
      (head, tail) = (u, v) if direc else (v, u)
70
      if v not in visited and rGraph[head][tail] > 0:
71
        q.put(v)
72
        chain[v] = u
73
        visited.add(v)
74
        if excesses[p] * excesses[v] < 0:
75
          return FindPath(chain, v, excesses, s, t, direc)
76

77
  if t in visited:
78
    return FindPath(chain, t, excesses, s, t, direc)
79
  if s in visited:
80
    return FindPath(chain, s, excesses, s, t, direc)
81

82
  return None, None
83

84

85
"""
86
Given a flow and the original graph, build the residual graph in the same dictionary format.
87
"""
88

89

90
def BuildRdGraph(flows, graph):
91
  rGraph = deepcopy(graph)
92
  for i in graph:
93
    for j in graph[i]:
94
      assert flows[i][j] <= graph[i][j]
95
      rGraph[i][j] -= flows[i][j]
96
      rGraph[j][i] += flows[i][j]
97
  return rGraph
98

99

100
"""
101
Given an infeasible flow, project it to a feasible one satisfying capacity constraints on the new graph
102
Iteratively calling the FeasRestoreIter() subroutine.
103
"""
104

105

106
def FeasProj(flows, graph, excesses, V, s, t):
107
  assert len(excesses) == V
108
  num_paths = 0
109
  total_path_len = 0
110
  rGraph = BuildRdGraph(flows, graph)
111
  proj_flows = deepcopy(flows)
112

113
  begin = time.time()
114
  for i in range(V):
115
    if i in [s, t] or excesses[i] == 0:
116
      continue
117
    while excesses[i] != 0:
118
      path, pathFlow = FeasRestoreIter(i, rGraph, excesses, s, t)
119
      if pathFlow is None:
120
        raise Exception(
121
            "Code has a bug. Didn't find a path during feasibility projection.")
122

123
      num_paths += 1
124
      total_path_len += (len(path) - 1)
125

126
      for j in range(len(path) - 1):
127
        u, v = path[j], path[j + 1]
128
        rGraph[u][v] -= pathFlow
129
        rGraph[v][u] += pathFlow
130
        proj_flows[u][v] += max(0, pathFlow - proj_flows[v][u])
131
        proj_flows[v][u] = max(0, proj_flows[v][u] - pathFlow)
132

133
      if path[0] not in [s, t]:
134
        excesses[path[0]] -= pathFlow
135
      if path[-1] not in [s, t]:
136
        excesses[path[-1]] += pathFlow
137
  end = time.time()
138

139
  print('feasiblity projection time:', end - begin)
140

141
  if num_paths > 0:
142
    feas_proj_avg_len = float(total_path_len) / num_paths
143
  else:
144
    feas_proj_avg_len = 0
145
  print('# of paths: %d, average length: %f' % (num_paths, feas_proj_avg_len))
146
  return proj_flows, rGraph, end - begin, num_paths, float(
147
      total_path_len) / num_paths
148

149

150
"""
151
The main function for warm-starting Ford-Fulkerson with a potentially infeasible flow.
152
  - First perform feasibility projection.
153
  - Then, based on the feasible flow, keep finding augmenting path and increasing the flow until termination.
154
"""
155

156

157
def WarmStartFlow(inf_flows, graph, excesses, V, s, t):
158
  begin = time.time()
159
  proj_flows, rGraph, proj_time, num_paths_proj, avg_len_proj = FeasProj(
160
      inf_flows, graph, excesses, V, s, t)
161
  proj_value = sum([proj_flows[s][v] for v in proj_flows[s]])
162
  max_flows, cuts, num_paths_aug, avg_len_aug = augmentingPath(
163
      graph, V, s, t, proj_flows, rGraph)
164
  end = time.time()
165
  print('warmstart total time: ', end - begin)
166
  print(
167
      'finding maximum flow # of augmenting path: %d, average augmenting path length: %f'
168
      % (num_paths_aug, avg_len_aug))
169
  return max_flows, cuts, proj_value, proj_time, end - begin, num_paths_proj, avg_len_proj, num_paths_aug, avg_len_aug
170

171

172
"""
173
Given another flow conserving flow and graph, round down the flow on every arc if it exceeds the capacity.
174
Returns the rounded flow and the excess at every node.
175
This function is not needed for warm-start itself but mostly for the experiment setting.
176
"""
177

178

179
def RoundDown(flows, graph, V, s, t):
180
  excesses = np.zeros(V, dtype=int)
181
  for i in range(V):
182
    for j in graph[i]:
183
      diff = max(0, flows[i][j] - graph[i][j])
184
      if i != s:
185
        excesses[i] += diff
186
      if j != t:
187
        excesses[j] -= diff
188
      flows[i][j] -= diff
189
  return excesses
190

191

192
"""
193
This function prints out the excesses/deficits on the graph so that we can see where they are.
194
"""
195

196

197
def ExcessGraph(excesses, size):
198
  excess_graph = np.zeros((size, size))
199
  for u in range(size * size):
200
    i, j = u // size, u % size
201
    excess_graph[i][j] = excesses[u]
202
  print(excess_graph)
203
  return
204

205

206
"""
207
The main function for running the full experiment with an existing set of seeds. The experiment does the following from image to image:
208
  - Construct the graph on the new image.
209
  - Take the old max flow solution from the last image, and round down according to the new capacities, thus generating excesses/deficits.
210
  - Run warm-start Ford-Fulkerson on the resulting infeasible flow.
211
"""
212

213

214
def Exp(args):
215
  folder, group, size, algo, loadseed = args.folder, args.group, args.size, args.algo, args.loadseed
216
  image_dir = folder + '/' + group + '_cropped'
217
  image_list = os.listdir(image_dir)
218
  V = size * size + 2
219
  SOURCE, SINK = V - 2, V - 1
220

221
  result_dir = folder + '/' + group + '_results' + '/'
222
  if not os.path.exists(result_dir):
223
    os.makedirs(result_dir)
224
  time_dir = result_dir + str(size) + '_time.txt'
225
  path_dir = result_dir + str(size) + '_path.txt'
226

227
  time_f = open(time_dir, 'w')
228
  time_f.write('image_name\tff\twarm_start\tfeas_proj\n')
229
  path_f = open(path_dir, 'w')
230
  path_f.write(
231
      'image_name\tflow_value\texcess\trecoverd_flow\tnum_aug_path\tavg_path_len\tnum_proj_path\tavg_proj_len\tnum_warm_start_path\tavg_warm_start_len\n'
232
  )
233

234
  num_images = len(image_list)
235
  old_flows = None
236
  for i in range(num_images):
237
    new_image = image_list[i]
238
    new_flows, min_cut, path_count, average_path_len, new_graph, ff_time = imageSegmentation(
239
        new_image, folder, group, (size, size), algo, loadseed)
240
    if old_flows is None:
241
      old_flows = deepcopy(new_flows)
242
      continue
243
    excesses = RoundDown(old_flows, new_graph, V, SOURCE, SINK)
244
    total_excess = max(
245
        np.sum(np.maximum(excesses, 0)), np.sum(np.maximum(-excesses, 0)))
246
    print('total excess/deficit to round down:' + str(total_excess))
247
    warmstart_flows, warmstart_cuts, proj_value, proj_time, warmstart_time, num_paths_proj, avg_len_proj, num_paths_aug, avg_len_aug = WarmStartFlow(
248
        old_flows, new_graph, excesses, V, SOURCE, SINK)
249

250
    # check that the warm start algo reaches optimality
251
    opt_flow_value = sum([new_flows[SOURCE][v] for v in new_flows[SOURCE]])
252
    ws_flow_value = sum(
253
        [warmstart_flows[SOURCE][v] for v in warmstart_flows[SOURCE]])
254
    assert opt_flow_value == ws_flow_value
255
    old_flows = deepcopy(new_flows)
256

257
    # write the results
258
    time_f.write(new_image.split('.')[0] + '\t')
259
    time_f.write('{}\t{}\t{}\n'.format(ff_time, warmstart_time, proj_time))
260
    time_f.flush()
261
    path_f.write(new_image.split('.')[0] + '\t')
262
    path_f.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(
263
        opt_flow_value, total_excess, proj_value, path_count,
264
        round(average_path_len, 4), num_paths_proj, round(avg_len_proj, 4),
265
        num_paths_aug, round(avg_len_aug, 4)))
266
    path_f.flush()
267

268
  time_f.close()
269
  path_f.close()
270
  return
271

272

273
if __name__ == '__main__':
274
  args = parseArgs()
275
  Exp(args)
276