source code of /small-scripts/foo.py

Last modified
Lines 259

Parent directory Download CGIread sitemap Main page

Quick links: __init__ dummy_progressbar extend gcd integral interactive_progressbar next_stage pascals_row test_precission update

  1. #!/usr/bin/python
  2. import os
  3. import sys
  4. import time
  5. class interactive_progressbar():
  6.     def __init__(self, *stages):
  7.         '''
  8.         '''
  9.         self.stage_counter = 0
  10.         self.n_stages = len(stages)
  11.         self.stage_title = stages
  12.         self.i = 0
  13.         self.n = 0
  14.         self.last_check = 0.0
  15.         self.p_past = (0, 0.0)
  16.         self.p_dpast = (0, 0.0)
  17.         self.tf_a = (0, 0.0)
  18.         self.tf_b = (0, 0.0)
  19.         self.tf_c = (0, 0.0)
  20.         self.start_time = 0.0
  21.         self.width = int(os.getenv('COLUMNS', '80')) - 1
  22.         self.prev_length = 0
  23.         # Print initial progress bar.
  24.         self.next_stage()
  25.     
  26.     def next_stage(self):
  27.         '''
  28.         '''
  29.         self.last_check = 0.0
  30.         self.p_past = (0, 0.0)
  31.         self.p_dpast = (0, 0.0)
  32.         self.tf_a = (0, 0.0)
  33.         self.tf_b = (0, 0.0)
  34.         self.tf_c = (0, 0.0)
  35.         self.start_time = time.time()
  36.         self.stage_counter += 1
  37.         sys.stderr.write('\b' * self.prev_length)
  38.         if self.stage_counter > self.n_stages:
  39.             line = '{}/{}: {}'.format(
  40.                 self.n_stages, self.n_stages,
  41.                 self.stage_title[self.n_stages - 1]
  42.             )
  43.         else:
  44.             line = '{}/{}: {}'.format(
  45.                 self.stage_counter, self.n_stages,
  46.                 self.stage_title[self.stage_counter - 1]
  47.             )
  48.         self.prev_length = len(line)
  49.         sys.stderr.write(line)
  50.         sys.stderr.flush()
  51.         if self.stage_counter > self.n_stages:
  52.             blank = self.width - len(line) - len(' [ Done! ]')
  53.             sys.stderr.write(' [' + ' '*blank + ' Done! ]\n')
  54.         
  55.     
  56.     def update(self, i, n):
  57.         if time.time() - self.last_check > 1:
  58.             self.last_check = time.time()
  59.             # Calculate time left
  60.             # uses a pessimistic approximation for O(n)=n^2
  61.             time_past = time.time() - self.start_time
  62.             d_time_past = (time_past - self.p_past[1]) / (i - self.p_past[0])
  63.             dd_time_past = (d_time_past - self.p_dpast[1]) / (i - self.p_dpast[0])
  64.             self.p_past = (i, time_past)
  65.             self.p_dpast = (i, d_time_past)
  66.             # t = f(i)
  67.             # f(x) = a*x^3 + b*x^2  + c*x
  68.             # f'(x) = 3*a*x^2 + 2*b*x + c
  69.             # f''(x) = 6*a*x + 2*b
  70.             a = (dd_time_past*i**2 - 2*d_time_past*i + 2*time_past) / (2*i**3)
  71.             b = -(dd_time_past*i**2 - 3*d_time_past*i + 3*time_past) / (i**2)
  72.             c = (dd_time_past*i**2 - 4*d_time_past*i + 9*time_past) / (2*i)
  73.             if a < 0:
  74.                 a = 0
  75.                 b = d_time_past
  76.                 c = time_past - d_time_past*i
  77.             if b < 0:
  78.                 b = 0
  79.                 c = time_past/i*n
  80.             a = 0
  81.             b = dd_time_past
  82.             c = d_time_past - dd_time_past*i
  83.             if i > 10:
  84.                 self.tf_a = self.tf_a[0] + 1, (self.tf_a[0]*self.tf_a[1] + a) / (self.tf_a[0] + 1)
  85.                 self.tf_b = self.tf_b[0] + 1, (self.tf_b[0]*self.tf_b[1] + b) / (self.tf_b[0] + 1)
  86.                 self.tf_c = self.tf_c[0] + 1, (self.tf_c[0]*self.tf_c[1] + c) / (self.tf_c[0] + 1)
  87.             a, b, c = self.tf_a[1], self.tf_b[1], self.tf_c[1]
  88.             time_left = int(round(a*n**3 + b*n**2 + c*n - time_past))
  89.             progress = time_past / (a*n**3 + b*n**2 + c*n + 1)
  90.             
  91.             # Format time left
  92.             time_left_t = time.gmtime(time_left)
  93.             if time_left < 600:
  94.                 time_left_s = time.strftime('%M:%S', time_left_t)[1:]
  95.             elif time_left < 3600:
  96.                 time_left_s = time.strftime('%M:%S', time_left_t)
  97.             elif time_left < 86400:
  98.                 time_left_s = time.strftime('%H:%M:%S', time_left_t)
  99.             else:
  100.                 time_left_s = str(time_left//86400) + '+'
  101.                 time_left_s += time.strftime('%H:%M', time_left_t)
  102.             # Format progress bar
  103.             stage_s = '{}/{}: {}'.format(
  104.                 self.stage_counter, self.n_stages,
  105.                 self.stage_title[self.stage_counter - 1]
  106.             )
  107.             crud = len(stage_s) + len(' [] ') + len(time_left_s)
  108.             bar = '=' * int(round(progress * (self.width - crud)))
  109.             blank = ' ' * int(round((1.0 - progress) * (self.width - crud)))
  110.             # Print bar
  111.             sys.stderr.write('\b' * self.prev_length)
  112.             line = '{} [{}{}] {}'.format(stage_s, bar, blank, time_left_s)
  113.             self.prev_length = len(line)
  114.             sys.stderr.write(line)
  115.             sys.stderr.flush()
  116.             
  117.            
  118. class dummy_progressbar():
  119.     def __init__(self, *args, **kwargs):
  120.         pass
  121.     def next_stage(self):
  122.         pass
  123.     def update(self, i, n):
  124.         pass
  125. def gcd(a, b):
  126.     while b:
  127.         a, b = b, a % b
  128.     return a
  129. def pascals_row(n, pb=None):
  130.     '''    
  131.     nCr(n, k) = pascals_row(n)[k]
  132.     
  133.     Save pascals_row(n) and use it as an optimized nCr function when
  134.     `n` does not vary.
  135.     '''
  136.     prev = curr = [1]
  137.     i = 0
  138.     while i < n:
  139.         i += 1
  140.         curr = [1]
  141.         j = 1
  142.         while j < i:
  143.             curr.append(prev[j-1] + prev[j])
  144.             j += 1
  145.         curr.append(1)
  146.         prev = curr
  147.         if pb is not None:
  148.             pb.update(i + 1, n)
  149.     return curr
  150. def extend(numerator_a, denominator_a, numerator_b, denominator_b):
  151.     '''
  152.     na, nb, d = extend(na, da, nb, db)
  153.     
  154.     Eg.
  155.     (3/4) + (2/3) = (9-8)/12
  156.     (3/4) - (2/3) = (9-8)/12
  157.     extend(3, 4, 2, 3) -> (9, 8, 12)
  158.     '''
  159.     # na1 = na0 * lcm(da0, db0)/da0
  160.     # nb1 = nb0 * lcm(da0, db0)/db0
  161.     # d1 = lcm(da0, db0)
  162.     # lcm(x, y) = (x*y)/gcd(x, y)
  163.     gcd_cache = gcd(denominator_a, denominator_b)
  164.     new_numerator_a = numerator_a * (denominator_b / gcd_cache)
  165.     new_numerator_b = numerator_b * (denominator_a / gcd_cache)
  166.     new_denominator = denominator_a * (denominator_b / gcd_cache)
  167.     # Dividing with gcd_cache before multiplying is a significant optimization.
  168.     # 4 times faster for test_precission(3000, 2000, 100)
  169.     return new_numerator_a, new_numerator_b, new_denominator
  170.     
  171. def integral(n, k, x, interactive):
  172.     if interactive:
  173.         pb = interactive_progressbar(
  174.             "Pascal's triangle", "Numerator", "Denominator", "Finalization"
  175.         )
  176.     else:
  177.         pb = dummy_progressbar()
  178.     
  179.     a, b, p_d = extend(k, n, 1, 2*x)
  180.     p_low_n = max(0, a - b)
  181.     if a + b > 1:
  182.         p_high_n = p_d
  183.     else:
  184.         p_high_n = a + b
  185.     ncr_row = pascals_row(n - k, pb)
  186.     pb.next_stage()
  187.     # Calculate the main numerator.
  188.     # sum(i=0, n-k, nCr(n-k, i) * (-1)^i * (high^(n-i+1)-low^(n-i+1))/(n-i+1))
  189.     numerator_sum_n, numerator_sum_d = 0, 1
  190.     p_low_exp_n = p_low_n**(k+1)
  191.     p_high_exp_n = p_high_n**(k+1)
  192.     p_exp_d =  p_d**(k+1)
  193.     i = n - k
  194.     while i >= 0:
  195.         # Calculate i'th term in the numerator sum.
  196.         temp_n = ncr_row[i] * (p_high_exp_n - p_low_exp_n)
  197.         temp_d = (n - i + 1) * p_exp_d
  198.         # Add to main numerator
  199.         numerator_sum_n, diff, numerator_sum_d = extend(
  200.             numerator_sum_n, numerator_sum_d,
  201.             temp_n, temp_d
  202.         )
  203.         if i % 2: # (-1)^i
  204.             numerator_sum_n -= diff
  205.         else:
  206.             numerator_sum_n += diff
  207.         # Next p_low^(n-i+1) nad p_high^(n-i+1)
  208.         p_low_exp_n *= p_low_n
  209.         p_high_exp_n *= p_high_n
  210.         p_exp_d *= p_d
  211.         # Loop
  212.         pb.update(n-k-i+1, n-k+1)
  213.         i -= 1
  214.     pb.next_stage()
  215.     # Calculate the main denominator.
  216.     # sum(i=0, n-k, nCr(n-k, i) * (-1)^i / (n-i+1))
  217.     denominator_sum_n, denominator_sum_d = 0, 1
  218.     i = 0
  219.     while i <= n - k:
  220.         temp_n = ncr_row[i]
  221.         temp_d = n - i + 1
  222.         # Add to main denominator
  223.         denominator_sum_n, diff, denominator_sum_d = extend(
  224.             denominator_sum_n, denominator_sum_d,
  225.             temp_n, temp_d
  226.         )
  227.         if i % 2: # (-1)^i
  228.             denominator_sum_n -= diff
  229.         else:
  230.             denominator_sum_n += diff
  231.         # Loop
  232.         pb.update(i+1, n-k+1)
  233.         i += 1
  234.     pb.next_stage()
  235.     # Compute integral and target.
  236.     integral_n = numerator_sum_n * denominator_sum_d
  237.     integral_d = denominator_sum_n * numerator_sum_d
  238.     # Compare
  239.     # NOTE:
  240.     # a/b > c/d  | *b *d
  241.     # ad > bc is wrong if d xor b is negative
  242.     # If d xor b is negative, ad < bc
  243.     if integral_d < 0:
  244.         integral_n, integral_d = -integral_n, -integral_d
  245.     # target_d is x and x is not negative, no test needed
  246.     pb.next_stage()
  247.     return integral_n, integral_d
  248. def test_precission(n, k, x, interactive=True):
  249.     integral_n, integral_d = integral(n, k, x, interactive)
  250.     target_n = x - 1
  251.     target_d = x
  252.     return (integral_n * target_d) > (target_n * integral_d)