contrib/network_model/regression2.py

   1 #!/usr/bin/python
   2 # This script takes the following command line parameters
   3 # 1) an input file containing 2 columns: message size and 1-way trip time
   4 # 2) the maximum relative error for a line segment
   5 # 3) the minimum number of points needed to justify adding a line segment
   6 # 4) the number of links
   7 # 5) the latency
   8 # 6) the bandwidth
   9
  10 import sys
  11
  12 def compute_regression(points):
  13     N = len(points)
  14
  15     if N < 1:
  16         return None
  17
  18     if N < 2:
  19         return (0, points[0][1])
  20
  21     Sx = Sy = Sxx = Syy = Sxy = 0.0
  22
  23     for x, y in points:
  24         Sx  += x
  25         Sy  += y
  26         Sxx += x*x
  27         Syy += y*y
  28         Sxy += x*y
  29     denom = Sxx * N - Sx * Sx
  30     # don't return 0 or negative values as a matter of principle...
  31     m = max(sys.float_info.min, (Sxy * N - Sy * Sx) / denom)
  32     b = max(sys.float_info.min, (Sxx * Sy - Sx * Sxy) / denom)
  33     return (m, b)
  34
  35 def compute_error(m, b, x, y):
  36     yp = m*x+b
  37     return abs(yp - y) / max(min(yp, y), sys.float_info.min)
  38
  39 def compute_max_error(m, b, points):
  40     max_error = 0.0
  41     for x, y in points:
  42         max_error = max(max_error, compute_error(m, b, x, y))
  43     return max_error
  44
  45 def get_max_error_point(m, b, points):
  46     max_error_index = -1
  47     max_error = 0.0
  48
  49     i = 0
  50     while i < len(points):
  51         x, y = points[i]
  52         error = compute_error(m, b, x, y)
  53         if error > max_error:
  54             max_error_index = i
  55             max_error = error
  56         i += 1
  57
  58     return (max_error_index, max_error)
  59
  60 infile_name = sys.argv[1]
  61 error_bound = float(sys.argv[2])
  62 min_seg_points = int(sys.argv[3])
  63 links = int(sys.argv[4])
  64 latency = float(sys.argv[5])
  65 bandwidth = float(sys.argv[6])
  66
  67 infile = open(infile_name, 'r')
  68
  69 # read datafile
  70 points = []
  71 for line in infile:
  72     fields = line.split()
  73     points.append((int(fields[0]), int(fields[1])))
  74 infile.close()
  75
  76 # should sort points by x values
  77 points.sort()
  78
  79 # break points up into segments
  80 pointsets = []
  81 lbi = 0
  82 while lbi < len(points):
  83     min_ubi = lbi
  84     max_ubi = len(points) - 1
  85     while max_ubi - min_ubi > 1:
  86         ubi = (min_ubi + max_ubi) / 2
  87         m, b = compute_regression(points[lbi:ubi+1])
  88         max_error = compute_max_error(m, b, points[lbi:ubi+1])
  89         if max_error > error_bound:
  90             max_ubi = ubi - 1
  91         else:
  92             min_ubi = ubi
  93     ubi = max_ubi
  94     if min_ubi < max_ubi:
  95         m, b = compute_regression(points[lbi:max_ubi+1])
  96         max_error = compute_max_error(m, b, points[lbi:max_ubi+1])
  97         if max_error > error_bound:
  98             ubi = min_ubi
  99     pointsets.append(points[lbi:ubi+1])
 100     lbi = ubi+1
 101
 102 # try to merge larger segments if possible and compute piecewise regression
 103 i = 0
 104 segments = []
 105 notoutliers = 0
 106 while i < len(pointsets):
 107     currpointset = []
 108     j = i
 109     while j < len(pointsets):
 110         newpointset = currpointset + pointsets[j]
 111         # if joining a small segment, we can delete bad points
 112         if len(pointsets[j]) < min_seg_points:
 113             k = 0
 114             while k < len(pointsets[j]):
 115                 m, b = compute_regression(newpointset)
 116                 max_error_index, max_error = get_max_error_point(m, b, newpointset)
 117                 if max_error <= error_bound:
 118                     break
 119                 del newpointset[max_error_index]
 120                 k += 1
 121             # only add new pointset if we had to delete fewer than its length
 122             # points
 123             if k < len(pointsets[j]):
 124                 i = j
 125                 currpointset = newpointset
 126         # otherwise, we just see if it works...
 127         else:
 128             m, b = compute_regression(newpointset)
 129             max_error = compute_max_error(m, b, newpointset)
 130             if max_error > error_bound:
 131                 break
 132             i = j
 133             currpointset = newpointset
 134         j += 1
 135     i += 1
 136     # outliers are ignored when constructing the piecewise funciton
 137     if len(currpointset) < min_seg_points:
 138         continue
 139     notoutliers += len(currpointset)
 140     m, b = compute_regression(currpointset)
 141     lb = min(x for x, y in currpointset)
 142     lat_factor = b / (1.0e6 * links * latency)
 143     bw_factor = 1.0e6 / (m * bandwidth)
 144     segments.append((lb, m, b, lat_factor, bw_factor))
 145
 146 outliers = len(points) - notoutliers
 147 segments.sort()
 148 segments.reverse()
 149
 150 print "/**--------- <copy/paste C code snippet in surf/network.c> -------------"
 151 print "  * produced by:"
 152 print "  *", " ".join(sys.argv)
 153 print "  * outliers:", outliers
 154 print "  * gnuplot: "
 155 print "    plot \"%s\" using 1:2 with lines title \"data\", \\" % (infile_name)
 156 for lb, m, b, lat_factor, bw_factor in segments:
 157     print "        (x >= %d) ? %g*x+%g : \\" % (lb, m, b)
 158 print "        1.0 with lines title \"piecewise function\""
 159 print "  *-------------------------------------------------------------------*/"
 160 print
 161 print "static double smpi_bandwidth_factor(double size)\n{\n"
 162 for lb, m, b, lat_factor, bw_factor in segments:
 163     print "    if (size >= %d) return %g;" % (lb, bw_factor)
 164 print "    return 1.0;\n}\n"
 165 print "static double smpi_latency_factor(double size)\n{\n"
 166 for lb, m, b, lat_factor, bw_factor in segments:
 167     print "    if (size >= %d) return %g;" % (lb, lat_factor)
 168 print "    return 1.0;\n}\n"
 169 print "/**--------- <copy/paste C code snippet in surf/network.c> -----------*/"