quantiztion of BSSEval values
501986
def _q(self, number, precision='.00001'):
        """"""
        if np.isinf(number):
            return np.nan
        else:
            return D(D(number).quantize(D(precision)))
947283
function bnpBitwiseTo(a,op,r) {
var i, f, m = Math.min(a.t,this.t);
for(i = 0; i < m; ++i) r.data[i] = op(this.data[i],a.data[i]);
if(a.t < this.t) {
 f = a.s&this.DM;
 for(i = m; i < this.t; ++i) r.data[i] = op(this.data[i],f);
 r.t = this.t;
} else {
 f = this.s&this.DM;
 for(i = m; i < a.t; ++i) r.data[i] = op(f,a.data[i]);
 r.t = a.t;
}
r.s = op(this.s,a.s);
r.clamp();
}
590572
float[] unquantize(){

    if(maptype==1||maptype==2){
      int quantvals;
      float mindel=float32_unpack(q_min);
      float delta=float32_unpack(q_delta);
      float[] r=new float[entries*dim];

      // maptype 1 and 2 both use a quantized value vector, but
      // different sizes
      switch(maptype){
        case 1:
          // most of the time, entries%dimensions == 0, but we need to be
          // well defined.  We define that the possible vales at each
          // scalar is values == entries/dim.  If entries%dim != 0, we'll
          // have 'too few' values (values*dim<entries), which means that
          // we'll have 'left over' entries; left over entries use zeroed
          // values (and are wasted).  So don't generate codebooks like that
          quantvals=maptype1_quantvals();
          for(int j=0; j<entries; j++){
            float last=0.f;
            int indexdiv=1;
            for(int k=0; k<dim; k++){
              int index=(j/indexdiv)%quantvals;
              float val=quantlist[index];
              val=Math.abs(val)*delta+mindel+last;
              if(q_sequencep!=0)
                last=val;
              r[j*dim+k]=val;
              indexdiv*=quantvals;
            }
          }
          break;
        case 2:
          for(int j=0; j<entries; j++){
            float last=0.f;
            for(int k=0; k<dim; k++){
              float val=quantlist[j*dim+k];
              //if((j*dim+k)==0){System.err.println(" | 0 -> "+val+" | ");}
              val=Math.abs(val)*delta+mindel+last;
              if(q_sequencep!=0)
                last=val;
              r[j*dim+k]=val;
              //if((j*dim+k)==0){System.err.println(" $ r[0] -> "+r[0]+" | ");}
            }
          }
          //System.err.println("\nr[0]="+r[0]);
      }
      return (r);
    }
    return (null);
  }
1208996
function(q, lowerTail, logp) {
            var ind, ret;

            logp = logp === true;
            lowerTail = lowerTail !== false;

            if (utils.hasNaN(q)) { return NaN; }
            if (q < xs[0]) {
               ret = lowerTail ? 0 : 1;
            } else if (q >= xs[xs.length - 1]) {
               ret = lowerTail ? 1 : 0;
            } else {
               ind = binSearch(xs, function(y) { return y <= q; });
               ret = lowerTail ? cumsLeft[ind] : cumsRight[ind];
            }
            ret = Math.min(ret, 1);
            return logp ? Math.log(ret) : ret;
         }
769366
function bnpRShiftTo(n, r) {
  r.s = this.s;
  var ds = Math.floor(n / this.DB);
  if (ds >= this.t) {
    r.t = 0;
    return;
  }
  var bs = n % this.DB;
  var cbs = this.DB - bs;
  var bm = (1 << bs) - 1;
  r[0] = this[ds] >> bs;
  for (var i = ds + 1; i < this.t; ++i) {
    r[i - ds - 1] |= (this[i] & bm) << cbs;
    r[i - ds] = this[i] >> bs;
  }
  if (bs > 0) r[this.t - ds - 1] |= (this.s & bm) << cbs;
  r.t = this.t - ds;
  r.clamp();
}
757602
function getSqSegDist(p, p1, p2) {

    var x = p1.x,
        y = p1.y,
        dx = p2.x - x,
        dy = p2.y - y;

    if (dx !== 0 || dy !== 0) {

        var t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy);

        if (t > 1) {
            x = p2.x;
            y = p2.y;

        } else if (t > 0) {
            x += dx * t;
            y += dy * t;
        }
    }

    dx = p.x - x;
    dy = p.y - y;

    return dx * dx + dy * dy;
}
415596
function ptsub(q, n) {
        var eps;

        if (n === 1 && q > 0.001 && q < 0.01) {
            eps = 1.0e-4;
        } else if (n === 2 && q < 0.0001) {
            eps = 1.0e-4;
        } else if (n === 1 && q < 0.001) {
            eps = 1.0e-2;
        } else {
            eps = 1.0e-5;
        }

        s = 10000.0;
        w = 0.0;
        while (true) {
            w = w + s;

            if (s <= eps) { return w; }
            if ((qe = 2.0 - p_t(n, w) * 2.0 - q) === 0.0) { return w; }
            if (qe < 0.0) {
                w = w - s;
                s = s / 10.0;
            }
        }
    }
105778
function (x, y) {
    var i, j, xi, Zi, Vi, lsb_Vi, w=sjcl.bitArray, xor=w._xor4;

    Zi = [0,0,0,0];
    Vi = y.slice(0);

    // Block size is 128 bits, run 128 times to get Z_128
    for (i=0; i<128; i++) {
      xi = (x[Math.floor(i/32)] & (1 << (31-i%32))) !== 0;
      if (xi) {
        // Z_i+1 = Z_i ^ V_i
        Zi = xor(Zi, Vi);
      }

      // Store the value of LSB(V_i)
      lsb_Vi = (Vi[3] & 1) !== 0;

      // V_i+1 = V_i >> 1
      for (j=3; j>0; j--) {
        Vi[j] = (Vi[j] >>> 1) | ((Vi[j-1]&1) << 31);
      }
      Vi[0] = Vi[0] >>> 1;

      // If LSB(V_i) is 1, V_i+1 = (V_i >> 1) ^ R
      if (lsb_Vi) {
        Vi[0] = Vi[0] ^ (0xe1 << 24);
      }
    }
    return Zi;
  }
1250871
function basym(a, b, x, y) {
      var z, p, q, bg, bgn, Ls, as, bs, cs, es, r, ab, sgn, abnp1,
         root2z, root2znm1, ser;

      p = a / (a + b);
      q = b / (a + b);
      bg = a < b ? Math.sqrt(q / a) : Math.sqrt(p / b); // betaGamma (11.3)
      bgn = 1 / bg;                        // betaGamma^n
      z = Math.sqrt(a * phi(x / p) + b * phi(y / q));   // (11.12)
      Ls = [];
      as = [];
      cs = [];
      es = [];
      ab = a <= b ? a / b : b / a;
      abnp1 = 1;  // (-1)^n * (a/b)^(n+1) or (b/a)^(n+1) for (11.4 - 11.5)
      root2z = Math.sqrt(2) * z;
      root2znm1 = 1; // (sqrt(2) * z) ^ (n - 1) for (11.21)
      sgn = -1; // (-1)^n
      ser = series(function(n) {
               var i, seriesTerms;

               bgn *= bg;
               abnp1 = -abnp1 * ab;
               sgn = -sgn;
               if (n === 0) {
                  es[0] = 1;
                  Ls[0] = Math.sqrt(Math.PI) / 4 * Math.exp(z * z) *
                           erfc(z);
               } else if (n === 1) {
                  as[1] = 2 / 3 * (1 - abnp1) * (a <= b ? q : sgn * p);
                  cs[1] = -0.5 * as[1];
                  es[1] = -cs[1];
                  Ls[1] = Math.pow(2, -1.5);
               } else {
                  root2znm1 *= root2z;
                  Ls[n] = Math.pow(2, -1.5) * root2znm1 +
                           (n - 1) * Ls[n - 2];     // (11.21)
                  as[n] = 2 / (n + 2) * (1 - abnp1) * (a <= b ? q : sgn * p);
                  r = (-n - 1) / 2;
                  seriesTerms = function(j) {
                        if (j === 0) { return 0; }
                        return ((i - j) * r - j) * bs[j] * as[i - j];
                     };
                  bs = [1, r * as[1]]; // first two bs for current r
                  for (i = 2; i <= n; i += 1) {
                     // add the series for b_n^r
                     bs[i] = r * as[i] + 1 / i * series(seriesTerms, i);
                  }
                  cs[n] = bs[n] / (n + 1); // (11.6)
                  es[n] = -series(function(k) {
                     return es[k] * cs[n - k];
                  }, n);
               }
               return es[n] * (Ls[n] * bgn);
             }, 100);
      return Math.log(2 / Math.sqrt(Math.PI)) -
             bcorr(a, b) -   // (11.20)
             z * z + Math.log(ser);
   }
1229789
def _se_all(self):
        """"""
        x = np.atleast_2d(self.x)
        err = np.atleast_1d(self.ms_err)
        se = np.sqrt(np.diagonal(np.linalg.inv(x.T @ x)) * err[:, None])
        return np.squeeze(se)
983745
function hslToRgb(h, s, l, a, out) {
		function hueToRgb(m1, m2, h) {
			h = h % 1;
			if (h < 0) {
				h += 1;
			}
			if (h < 1 / 6) {
				return m1 + (m2 - m1) * h * 6;
			}
			if (h < 1 / 2) {
				return m2;
			}
			if (h < 2 / 3) {
				return m1 + (m2 - m1) * (2/3 - h) * 6;
			}
			return m1;
		}

		var m1, m2;
		if (l < 0.5) {
			m2 = l * (s + 1);
		} else {
			m2 = l + s - l * s;
		}
		m1 = l * 2 - m2;

		if (!out) {
			out = [];
		}

		out[0] = hueToRgb(m1, m2, h + 1/3);
		out[1] = hueToRgb(m1, m2, h);
		out[2] = hueToRgb(m1, m2, h - 1/3);
		out[3] = a;

		return out;
	}
313958
def bootstrap_ts(y, func, B=1000, b=3):
    
    beta_star = np.empty(B)
    z = y
    z_star = np.empty(len(z))
    for boot_i in range(B):
        for block_i, start in enumerate(np.random.randint(len(z) - b + 1, size=len(z) / b)):
            z_star[block_i * b:(block_i + 1) * b] = z[start:start + b]
            beta_star[boot_i] = func(z_star)
    return beta_star
483798
function sqSegSegDist(x0, y0, x1, y1, x2, y2, x3, y3) {
    var ux = x1 - x0;
    var uy = y1 - y0;
    var vx = x3 - x2;
    var vy = y3 - y2;
    var wx = x0 - x2;
    var wy = y0 - y2;
    var a = ux * ux + uy * uy;
    var b = ux * vx + uy * vy;
    var c = vx * vx + vy * vy;
    var d = ux * wx + uy * wy;
    var e = vx * wx + vy * wy;
    var D = a * c - b * b;

    var sc, sN, tc, tN;
    var sD = D;
    var tD = D;

    if (D === 0) {
        sN = 0;
        sD = 1;
        tN = e;
        tD = c;
    } else {
        sN = b * e - c * d;
        tN = a * e - b * d;
        if (sN < 0) {
            sN = 0;
            tN = e;
            tD = c;
        } else if (sN > sD) {
            sN = sD;
            tN = e + b;
            tD = c;
        }
    }

    if (tN < 0.0) {
        tN = 0.0;
        if (-d < 0.0) sN = 0.0;
        else if (-d > a) sN = sD;
        else {
            sN = -d;
            sD = a;
        }
    } else if (tN > tD) {
        tN = tD;
        if ((-d + b) < 0.0) sN = 0;
        else if (-d + b > a) sN = sD;
        else {
            sN = -d + b;
            sD = a;
        }
    }

    sc = sN === 0 ? 0 : sN / sD;
    tc = tN === 0 ? 0 : tN / tD;

    var cx = (1 - sc) * x0 + sc * x1;
    var cy = (1 - sc) * y0 + sc * y1;
    var cx2 = (1 - tc) * x2 + tc * x3;
    var cy2 = (1 - tc) * y2 + tc * y3;
    var dx = cx2 - cx;
    var dy = cy2 - cy;

    return dx * dx + dy * dy;
}
738702
def sde(self):
        
        variance = float(self.variance.values)
        lengthscale = float(self.lengthscale)

        F  = np.array(((-1.0/lengthscale,),))
        L  = np.array(((1.0,),))
        Qc = np.array( ((2.0*variance/lengthscale,),) )
        H = np.array(((1.0,),))
        Pinf = np.array(((variance,),))
        P0 = Pinf.copy()

        dF = np.zeros((1,1,2));
        dQc = np.zeros((1,1,2));
        dPinf = np.zeros((1,1,2));

        dF[:,:,0] = 0.0
        dF[:,:,1] = 1.0/lengthscale**2

        dQc[:,:,0] = 2.0/lengthscale
        dQc[:,:,1] = -2.0*variance/lengthscale**2

        dPinf[:,:,0] = 1.0
        dPinf[:,:,1] = 0.0

        dP0 = dPinf.copy()

        return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
102025
def ρ(self, e):
        
        d = self._dichalcogenide
        at, λ = d.at, d.λ
        return abs(2 * e - λ) * (2 * at**2)**(-1)
1040421
double qobs(int dummy, double[] xt) {
        int x = (int) (xt[0]);
        double p = Math.random();
        double[] distr = accumP[x];
        int i, left = 0, right = distr.length - 1;
        // Find and return least i s.t. p <= distr[i]
        while (left <= right) {
            // Here distr[left-1] <= p <= distr[right+1]
            i = (left + right) / 2;
            if (p < distr[i]) {
                right = i - 1;
            } else if (distr[i] < p) {
                left = i + 1;
            } else {
                return i;
            }
        }
        // Now p <= distr[left]
        return left;
    }
408877
def get_dorefa(bitW, bitA, bitG):
    
    def quantize(x, k):
        n = float(2 ** k - 1)

        @tf.custom_gradient
        def _quantize(x):
            return tf.round(x * n) / n, lambda dy: dy

        return _quantize(x)

    def fw(x):
        if bitW == 32:
            return x

        if bitW == 1:   # BWN
            E = tf.stop_gradient(tf.reduce_mean(tf.abs(x)))

            @tf.custom_gradient
            def _sign(x):
                return tf.where(tf.equal(x, 0), tf.ones_like(x), tf.sign(x / E)) * E, lambda dy: dy

            return _sign(x)

        x = tf.tanh(x)
        x = x / tf.reduce_max(tf.abs(x)) * 0.5 + 0.5
        return 2 * quantize(x, bitW) - 1

    def fa(x):
        if bitA == 32:
            return x
        return quantize(x, bitA)

    def fg(x):
        if bitG == 32:
            return x

        @tf.custom_gradient
        def _identity(input):
            def grad_fg(x):
                rank = x.get_shape().ndims
                assert rank is not None
                maxx = tf.reduce_max(tf.abs(x), list(range(1, rank)), keep_dims=True)
                x = x / maxx
                n = float(2**bitG - 1)
                x = x * 0.5 + 0.5 + tf.random_uniform(
                    tf.shape(x), minval=-0.5 / n, maxval=0.5 / n)
                x = tf.clip_by_value(x, 0.0, 1.0)
                x = quantize(x, bitG) - 0.5
                return x * maxx * 2

            return input, grad_fg

        return _identity(x)
    return fw, fa, fg