Gamma function from picomath.org
see http://picomath.org/php/gamma.php.html for commented
implementation

TODO: These should probably move outside of NlpTools together
with the Random namespace and form a nice php math library
765130
private function gamma($x)
    {
        $gamma = 0.577215664901532860606512090; # Euler's gamma constant
        if ($x < 0.001) {
            return 1.0/($x*(1.0 + $gamma*$x));
        }
        if ($x < 12.0) {
            # The algorithm directly approximates gamma over (1,2) and uses
            # reduction identities to reduce other arguments to this interval.
            $y = $x;
            $n = 0;
            $arg_was_less_than_one = ($y < 1.0);
            # Add or subtract integers as necessary to bring y into (1,2)
            # Will correct for this below
            if ($arg_was_less_than_one) {
                $y += 1.0;
            } else {
                $n = floor($y) - 1;  # will use n later
                $y -= $n;
            }
            # numerator coefficients for approximation over the interval (1,2)
            $p =
            array(
                -1.71618513886549492533811E+0,
                 2.47656508055759199108314E+1,
                -3.79804256470945635097577E+2,
                 6.29331155312818442661052E+2,
                 8.66966202790413211295064E+2,
                -3.14512729688483675254357E+4,
                -3.61444134186911729807069E+4,
                 6.64561438202405440627855E+4
            );

            # denominator coefficients for approximation over the interval (1,2)
            $q =
            array(
                -3.08402300119738975254353E+1,
                 3.15350626979604161529144E+2,
                -1.01515636749021914166146E+3,
                -3.10777167157231109440444E+3,
                 2.25381184209801510330112E+4,
                 4.75584627752788110767815E+3,
                -1.34659959864969306392456E+5,
                -1.15132259675553483497211E+5
            );

            $num = 0.0;
            $den = 1.0;

            $z = $y - 1;
            for ($i = 0; $i < 8; $i++) {
                $num = ($num + $p[$i])*$z;
                $den = $den*$z + $q[$i];
            }
            $result = $num/$den + 1.0;

            # Apply correction if argument was not initially in (1,2)
            if ($arg_was_less_than_one) {
                # Use identity gamma(z) = gamma(z+1)/z
                # The variable "result" now holds gamma of the original y + 1
                # Thus we use y-1 to get back the orginal y.
                $result /= ($y-1.0);
            } else {
                # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
                for ($i = 0; $i < $n; $i++) {
                    $result *= $y++;
                }
            }

            return $result;
        }

        ###########################################################################
        # Third interval: [12, infinity)

        if ($x > 171.624) {
            # Correct answer too large to display.

            return Double.POSITIVE_INFINITY;
        }

        return exp($this->log_gamma($x));
    }
922125
function toLatLng()
	{
		$wgs84 = new ReferenceEllipsoid(ReferenceEllipsoid::WGS84_MAJ, ReferenceEllipsoid::WGS84_MIN);

		$UTM_F0	 = 0.9996;
		$a = $wgs84->maj;
		$eSquared = $wgs84->ecc;
		$ePrimeSquared = $eSquared / (1.0 - $eSquared);
		$e1 = (1 - sqrt(1 - $eSquared)) / (1 + sqrt(1 - $eSquared));

		$x = $this->easting - 500000.0;;
		$y = $this->northing;

		$zoneNumber = $this->lngZone;
		$zoneLetter = $this->latZone;

		$longitudeOrigin = ($zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;

		// Correct y for southern hemisphere
		if ((ord($zoneLetter) - ord("N")) < 0) {
			$y -= 10000000.0;
		}

		$m = $y / $UTM_F0;

		$mu = $m
			/ ($a
			* (1.0
				- $eSquared / 4.0
				- 3.0 * $eSquared * $eSquared / 64.0
				- 5.0
				* pow($eSquared, 3.0)
				/ 256.0));

		$phi1Rad = $mu
			+ (3.0 * $e1 / 2.0 - 27.0 * pow($e1, 3.0) / 32.0) * sin(2.0 * $mu)
			+ (21.0 * $e1 * $e1 / 16.0 - 55.0 * pow($e1, 4.0) / 32.0)
			* sin(4.0 * $mu)
			+ (151.0 * pow($e1, 3.0) / 96.0) * sin(6.0 * $mu);

		$n = $a / sqrt(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad));

		$t = tan($phi1Rad) * tan($phi1Rad);
		$c = $ePrimeSquared * cos($phi1Rad) * cos($phi1Rad);

		$r = $a
			* (1.0 - $eSquared)
			/ pow(
			1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad),
			1.5);

		$d = $x / ($n * $UTM_F0);

		$latitude = ($phi1Rad
			- ($n * tan($phi1Rad) / $r)
			* ($d * $d / 2.0
				- (5.0
				+ (3.0 * $t)
				+ (10.0 * $c)
				- (4.0 * $c * $c)
				- (9.0 * $ePrimeSquared))
				* pow($d, 4.0)
				/ 24.0
				+ (61.0
				+ (90.0 * $t)
				+ (298.0 * $c)
				+ (45.0 * $t * $t)
				- (252.0 * $ePrimeSquared)
				- (3.0 * $c * $c))
				* pow($d, 6.0)
				/ 720.0)) * (180.0 / pi());

		$longitude = $longitudeOrigin + (($d
			- (1.0 + 2.0 * $t + $c) * pow($d, 3.0) / 6.0
			+ (5.0
			- (2.0 * $c)
			+ (28.0 * $t)
			- (3.0 * $c * $c)
			+ (8.0 * $ePrimeSquared)
			+ (24.0 * $t * $t))
			* pow($d, 5.0)
			/ 120.0)
			/ cos($phi1Rad)) * (180.0 / pi());

		return new LatLng($latitude, $longitude);
	}
648775
function WriteBarcode2($code, $x = '', $y = '', $size = 1, $height = 1, $bgcol = false, $col = false, $btype = 'IMB', $print_ratio = '', $k = 1)
	{
		if (empty($code)) {
			return;
		}

		$this->barcode = new Barcode();
		$arrcode = $this->barcode->getBarcodeArray($code, $btype, $print_ratio);

		if (empty($x)) {
			$x = $this->x;
		}

		if (empty($y)) {
			$y = $this->y;
		}

		$prevDrawColor = $this->DrawColor;
		$prevTextColor = $this->TextColor;
		$prevFillColor = $this->FillColor;
		$lw = $this->LineWidth;
		$this->SetLineWidth(0.01);
		$size /= $k; // in case resized in a table
		$xres = $arrcode['nom-X'] * $size;

		if ($btype === 'IMB' || $btype === 'RM4SCC' || $btype === 'KIX' || $btype === 'POSTNET' || $btype === 'PLANET') {
			$llm = $arrcode['quietL'] / $k; // Left Quiet margin
			$rlm = $arrcode['quietR'] / $k; // Right Quiet margin
			$tlm = $blm = $arrcode['quietTB'] / $k;
			$height = 1;  // Overrides
		} elseif (in_array($btype, ['C128A', 'C128B', 'C128C', 'EAN128A', 'EAN128B', 'EAN128C', 'C39', 'C39+', 'C39E', 'C39E+', 'S25', 'S25+', 'I25', 'I25+', 'I25B', 'I25B+', 'C93', 'MSI', 'MSI+', 'CODABAR', 'CODE11'])) {
			$llm = $arrcode['lightmL'] * $xres; // Left Quiet margin
			$rlm = $arrcode['lightmR'] * $xres; // Right Quiet margin
			$tlm = $blm = $arrcode['lightTB'] * $xres * $height;
		}

		$bcw = ($arrcode["maxw"] * $xres);
		$fbw = $bcw + $llm + $rlm;  // Full barcode width incl. light margins

		$bch = ($arrcode["nom-H"] * $size * $height);
		$fbh = $bch + $tlm + $blm;  // Full barcode height

		// PRINT border background color
		$xpos = $x;
		$ypos = $y;

		if ($col) {
			$this->SetDColor($col);
			$this->SetTColor($col);
		} else {
			$this->SetDColor($this->colorConverter->convert(0, $this->PDFAXwarnings));
			$this->SetTColor($this->colorConverter->convert(0, $this->PDFAXwarnings));
		}

		if ($bgcol) {
			$this->SetFColor($bgcol);
		} else {
			$this->SetFColor($this->colorConverter->convert(255, $this->PDFAXwarnings));
		}

		// PRINT BARS
		if ($col) {
			$this->SetFColor($col);
		} else {
			$this->SetFColor($this->colorConverter->convert(0, $this->PDFAXwarnings));
		}
		$xpos = $x + $llm;

		if ($arrcode !== false) {
			foreach ($arrcode["bcode"] as $v) {
				$bw = ($v["w"] * $xres);
				if ($v["t"]) {
					$ypos = $y + $tlm + ($bch * $v['p'] / $arrcode['maxh']);
					$this->Rect($xpos, $ypos, $bw, ($v['h'] * $bch / $arrcode['maxh']), 'F');
				}
				$xpos += $bw;
			}
		}

		// PRINT BEARER BARS
		if ($btype == 'I25B' || $btype == 'I25B+') {
			$this->Rect($x, $y, $fbw, ($arrcode['lightTB'] * $xres * $height), 'F');
			$this->Rect($x, $y + $tlm + $bch, $fbw, ($arrcode['lightTB'] * $xres * $height), 'F');
		}

		// Restore **************
		$this->DrawColor = $prevDrawColor;
		$this->TextColor = $prevTextColor;
		$this->FillColor = $prevFillColor;
		$this->SetLineWidth($lw);
		$this->SetY($y);
	}
1141866
def gammalnStirling(z):
    
    return (0.5 * (np.log(2. * np.pi) - np.log(z))) \
           + (z * (np.log(z + (1. / ((12. * z) - (1. / (10. * z))))) - 1.))
875400
function NeuQuant(pixels, samplefac) {
  var network; // int[netsize][4]
  var netindex; // for network lookup - really 256

  // bias and freq arrays for learning
  var bias;
  var freq;
  var radpower;

  /*
    Private Method: init

    sets up arrays
  */
  function init() {
    network = [];
    netindex = new Int32Array(256);
    bias = new Int32Array(netsize);
    freq = new Int32Array(netsize);
    radpower = new Int32Array(netsize >> 3);

    var i, v;
    for (i = 0; i < netsize; i++) {
      v = (i << (netbiasshift + 8)) / netsize;
      network[i] = new Float64Array([v, v, v, 0]);
      //network[i] = [v, v, v, 0]
      freq[i] = intbias / netsize;
      bias[i] = 0;
    }
  }

  /*
    Private Method: unbiasnet

    unbiases network to give byte values 0..255 and record position i to prepare for sort
  */
  function unbiasnet() {
    for (var i = 0; i < netsize; i++) {
      network[i][0] >>= netbiasshift;
      network[i][1] >>= netbiasshift;
      network[i][2] >>= netbiasshift;
      network[i][3] = i; // record color number
    }
  }

  /*
    Private Method: altersingle

    moves neuron *i* towards biased (b,g,r) by factor *alpha*
  */
  function altersingle(alpha, i, b, g, r) {
    network[i][0] -= (alpha * (network[i][0] - b)) / initalpha;
    network[i][1] -= (alpha * (network[i][1] - g)) / initalpha;
    network[i][2] -= (alpha * (network[i][2] - r)) / initalpha;
  }

  /*
    Private Method: alterneigh

    moves neurons in *radius* around index *i* towards biased (b,g,r) by factor *alpha*
  */
  function alterneigh(radius, i, b, g, r) {
    var lo = Math.abs(i - radius);
    var hi = Math.min(i + radius, netsize);

    var j = i + 1;
    var k = i - 1;
    var m = 1;

    var p, a;
    while ((j < hi) || (k > lo)) {
      a = radpower[m++];

      if (j < hi) {
        p = network[j++];
        p[0] -= (a * (p[0] - b)) / alpharadbias;
        p[1] -= (a * (p[1] - g)) / alpharadbias;
        p[2] -= (a * (p[2] - r)) / alpharadbias;
      }

      if (k > lo) {
        p = network[k--];
        p[0] -= (a * (p[0] - b)) / alpharadbias;
        p[1] -= (a * (p[1] - g)) / alpharadbias;
        p[2] -= (a * (p[2] - r)) / alpharadbias;
      }
    }
  }

  /*
    Private Method: contest

    searches for biased BGR values
  */
  function contest(b, g, r) {
    /*
      finds closest neuron (min dist) and updates freq
      finds best neuron (min dist-bias) and returns position
      for frequently chosen neurons, freq[i] is high and bias[i] is negative
      bias[i] = gamma * ((1 / netsize) - freq[i])
    */

    var bestd = ~(1 << 31);
    var bestbiasd = bestd;
    var bestpos = -1;
    var bestbiaspos = bestpos;

    var i, n, dist, biasdist, betafreq;
    for (i = 0; i < netsize; i++) {
      n = network[i];

      dist = Math.abs(n[0] - b) + Math.abs(n[1] - g) + Math.abs(n[2] - r);
      if (dist < bestd) {
        bestd = dist;
        bestpos = i;
      }

      biasdist = dist - ((bias[i]) >> (intbiasshift - netbiasshift));
      if (biasdist < bestbiasd) {
        bestbiasd = biasdist;
        bestbiaspos = i;
      }

      betafreq = (freq[i] >> betashift);
      freq[i] -= betafreq;
      bias[i] += (betafreq << gammashift);
    }

    freq[bestpos] += beta;
    bias[bestpos] -= betagamma;

    return bestbiaspos;
  }

  /*
    Private Method: inxbuild

    sorts network and builds netindex[0..255]
  */
  function inxbuild() {
    var i, j, p, q, smallpos, smallval, previouscol = 0, startpos = 0;
    for (i = 0; i < netsize; i++) {
      p = network[i];
      smallpos = i;
      smallval = p[1]; // index on g
      // find smallest in i..netsize-1
      for (j = i + 1; j < netsize; j++) {
        q = network[j];
        if (q[1] < smallval) { // index on g
          smallpos = j;
          smallval = q[1]; // index on g
        }
      }
      q = network[smallpos];
      // swap p (i) and q (smallpos) entries
      if (i != smallpos) {
        j = q[0];   q[0] = p[0];   p[0] = j;
        j = q[1];   q[1] = p[1];   p[1] = j;
        j = q[2];   q[2] = p[2];   p[2] = j;
        j = q[3];   q[3] = p[3];   p[3] = j;
      }
      // smallval entry is now in position i

      if (smallval != previouscol) {
        netindex[previouscol] = (startpos + i) >> 1;
        for (j = previouscol + 1; j < smallval; j++)
          netindex[j] = i;
        previouscol = smallval;
        startpos = i;
      }
    }
    netindex[previouscol] = (startpos + maxnetpos) >> 1;
    for (j = previouscol + 1; j < 256; j++)
      netindex[j] = maxnetpos; // really 256
  }

  /*
    Private Method: inxsearch

    searches for BGR values 0..255 and returns a color index
  */
  function inxsearch(b, g, r) {
    var a, p, dist;

    var bestd = 1000; // biggest possible dist is 256*3
    var best = -1;

    var i = netindex[g]; // index on g
    var j = i - 1; // start at netindex[g] and work outwards

    while ((i < netsize) || (j >= 0)) {
      if (i < netsize) {
        p = network[i];
        dist = p[1] - g; // inx key
        if (dist >= bestd) i = netsize; // stop iter
        else {
          i++;
          if (dist < 0) dist = -dist;
          a = p[0] - b; if (a < 0) a = -a;
          dist += a;
          if (dist < bestd) {
            a = p[2] - r; if (a < 0) a = -a;
            dist += a;
            if (dist < bestd) {
              bestd = dist;
              best = p[3];
            }
          }
        }
      }
      if (j >= 0) {
        p = network[j];
        dist = g - p[1]; // inx key - reverse dif
        if (dist >= bestd) j = -1; // stop iter
        else {
          j--;
          if (dist < 0) dist = -dist;
          a = p[0] - b; if (a < 0) a = -a;
          dist += a;
          if (dist < bestd) {
            a = p[2] - r; if (a < 0) a = -a;
            dist += a;
            if (dist < bestd) {
              bestd = dist;
              best = p[3];
            }
          }
        }
      }
    }

    return best;
  }

  /*
    Private Method: learn

    "Main Learning Loop"
  */
  function learn() {
    var i;

    var lengthcount = pixels.length;
    var alphadec = 30 + ((samplefac - 1) / 3);
    var samplepixels = lengthcount / (3 * samplefac);
    var delta = ~~(samplepixels / ncycles);
    var alpha = initalpha;
    var radius = initradius;

    var rad = radius >> radiusbiasshift;

    if (rad <= 1) rad = 0;
    for (i = 0; i < rad; i++)
      radpower[i] = alpha * (((rad * rad - i * i) * radbias) / (rad * rad));

    var step;
    if (lengthcount < minpicturebytes) {
      samplefac = 1;
      step = 3;
    } else if ((lengthcount % prime1) !== 0) {
      step = 3 * prime1;
    } else if ((lengthcount % prime2) !== 0) {
      step = 3 * prime2;
    } else if ((lengthcount % prime3) !== 0)  {
      step = 3 * prime3;
    } else {
      step = 3 * prime4;
    }

    var b, g, r, j;
    var pix = 0; // current pixel

    i = 0;
    while (i < samplepixels) {
      b = (pixels[pix] & 0xff) << netbiasshift;
      g = (pixels[pix + 1] & 0xff) << netbiasshift;
      r = (pixels[pix + 2] & 0xff) << netbiasshift;

      j = contest(b, g, r);

      altersingle(alpha, j, b, g, r);
      if (rad !== 0) alterneigh(rad, j, b, g, r); // alter neighbours

      pix += step;
      if (pix >= lengthcount) pix -= lengthcount;

      i++;

      if (delta === 0) delta = 1;
      if (i % delta === 0) {
        alpha -= alpha / alphadec;
        radius -= radius / radiusdec;
        rad = radius >> radiusbiasshift;

        if (rad <= 1) rad = 0;
        for (j = 0; j < rad; j++)
          radpower[j] = alpha * (((rad * rad - j * j) * radbias) / (rad * rad));
      }
    }
  }

  /*
    Method: buildColormap

    1. initializes network
    2. trains it
    3. removes misconceptions
    4. builds colorindex
  */
  function buildColormap() {
    init();
    learn();
    unbiasnet();
    inxbuild();
  }
  this.buildColormap = buildColormap;

  /*
    Method: getColormap

    builds colormap from the index

    returns array in the format:

    >
    > [r, g, b, r, g, b, r, g, b, ..]
    >
  */
  function getColormap() {
    var map = [];
    var index = [];

    for (var i = 0; i < netsize; i++)
      index[network[i][3]] = i;

    var k = 0;
    for (var l = 0; l < netsize; l++) {
      var j = index[l];
      map[k++] = (network[j][0]);
      map[k++] = (network[j][1]);
      map[k++] = (network[j][2]);
    }
    return map;
  }
  this.getColormap = getColormap;

  /*
    Method: lookupRGB

    looks for the closest *r*, *g*, *b* color in the map and
    returns its index
  */
  this.lookupRGB = inxsearch;
}
292189
function gamma1_64(x)
	{
		var rotr19 = rotr_64(x, 19), rotr61 = rotr_64(x, 61),
			shr6 = shr_64(x, 6);

		return new Int_64(
				rotr19.highOrder ^ rotr61.highOrder ^ shr6.highOrder,
				rotr19.lowOrder ^ rotr61.lowOrder ^ shr6.lowOrder
			);
	}
795457
def rdirichlet(theta, size=1):
    
    gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
    if size > 1 and np.size(theta) > 1:
        return (gammas.T / gammas.sum(1))[:-1].T
    elif np.size(theta) > 1:
        return (gammas[0] / gammas[0].sum())[:-1]
    else:
        return 1.
337577
def afterglow(self, src=None, event=None, dst=None, **kargs):
        """"""
        if src is None:
            src = lambda x: x['IP'].src
        if event is None:
            event = lambda x: x['IP'].dport
        if dst is None:
            dst = lambda x: x['IP'].dst
        sl = {}
        el = {}
        dl = {}
        for i in self.res:
            try:
                s, e, d = src(i), event(i), dst(i)
                if s in sl:
                    n, lst = sl[s]
                    n += 1
                    if e not in lst:
                        lst.append(e)
                    sl[s] = (n, lst)
                else:
                    sl[s] = (1, [e])
                if e in el:
                    n, lst = el[e]
                    n += 1
                    if d not in lst:
                        lst.append(d)
                    el[e] = (n, lst)
                else:
                    el[e] = (1, [d])
                dl[d] = dl.get(d, 0) + 1
            except Exception:
                continue

        import math

        def normalize(n):
            return 2 + math.log(n) / 4.0

        def minmax(x):
            m, M = reduce(lambda a, b: (min(a[0], b[0]), max(a[1], b[1])),
                          ((a, a) for a in x))
            if m == M:
                m = 0
            if M == 0:
                M = 1
            return m, M

        mins, maxs = minmax(x for x, _ in six.itervalues(sl))
        mine, maxe = minmax(x for x, _ in six.itervalues(el))
        mind, maxd = minmax(six.itervalues(dl))

        gr = 'digraph "afterglow" {\n\tedge [len=2.5];\n'

        gr += "# src nodes\n"
        for s in sl:
            n, _ = sl[s]
            n = 1 + float(n - mins) / (maxs - mins)
            gr += '"src.%s" [label = "%s", shape=box, fillcolor="#FF0000", style=filled, fixedsize=1, height=%.2f,width=%.2f];\n' % (repr(s), repr(s), n, n)  # noqa: E501
        gr += "# event nodes\n"
        for e in el:
            n, _ = el[e]
            n = n = 1 + float(n - mine) / (maxe - mine)
            gr += '"evt.%s" [label = "%s", shape=circle, fillcolor="#00FFFF", style=filled, fixedsize=1, height=%.2f, width=%.2f];\n' % (repr(e), repr(e), n, n)  # noqa: E501
        for d in dl:
            n = dl[d]
            n = n = 1 + float(n - mind) / (maxd - mind)
            gr += '"dst.%s" [label = "%s", shape=triangle, fillcolor="#0000ff", style=filled, fixedsize=1, height=%.2f, width=%.2f];\n' % (repr(d), repr(d), n, n)  # noqa: E501

        gr += "###\n"
        for s in sl:
            n, lst = sl[s]
            for e in lst:
                gr += ' "src.%s" -> "evt.%s";\n' % (repr(s), repr(e))
        for e in el:
            n, lst = el[e]
            for d in lst:
                gr += ' "evt.%s" -> "dst.%s";\n' % (repr(e), repr(d))

        gr += "}"
        return do_graph(gr, **kargs)
565089
private static float[][] generateGaussianBlurKernels(int level) {
		float[][] pascalsTriangle = generatePascalsTriangle(level);
		float[][] gaussianTriangle = new float[pascalsTriangle.length][];
		for (int i = 0; i < gaussianTriangle.length; i++) {
			float total = 0.0f;
			gaussianTriangle[i] = new float[pascalsTriangle[i].length];
			for (int j = 0; j < pascalsTriangle[i].length; j++)
				total += pascalsTriangle[i][j];
			float coefficient = 1 / total;
			for (int j = 0; j < pascalsTriangle[i].length; j++)
				gaussianTriangle[i][j] = coefficient * pascalsTriangle[i][j];
		}
		return gaussianTriangle;
	}

	/**
	 * Generate Pascal's triangle
	 * 
	 * @param level The level of the triangle to generate
	 * @return The Pascal's triangle kernel
	 */
	private static float[][] generatePascalsTriangle(int level) {
		if (level < 2) level = 2;
		float[][] triangle = new float[level][];
		triangle[0] = new float[1];
		triangle[1] = new float[2];
		triangle[0][0] = 1.0f;
		triangle[1][0] = 1.0f;
		triangle[1][1] = 1.0f;
		for (int i = 2; i < level; i++) {
			triangle[i] = new float[i + 1];
			triangle[i][0] = 1.0f;
			triangle[i][i] = 1.0f;
			for (int j = 1; j < triangle[i].length - 1; j++)
				triangle[i][j] = triangle[i - 1][j - 1] + triangle[i - 1][j];
		}
		return triangle;
	}
}
1235223
def transQuad(length=0.0, k1=0.0, gamma=None):
    
    m = np.eye(6, 6, dtype=np.float64)
    if length == 0.0:
        print("warning: 'length' should be a positive float number.")
    elif gamma is not None and gamma != 0.0:
        if k1 == 0:
            print("warning: 'k1' should be a positive float number.")
            m[0, 1] = m[2, 3] = 1.0
            m[4, 5] = float(length) / gamma / gamma
        else:
            sqrtk = np.sqrt(complex(k1))
            sqrtkl = sqrtk * length
            m[0, 0] = m[1, 1] = (np.cos(sqrtkl)).real
            m[0, 1] = (np.sin(sqrtkl) / sqrtk).real
            m[1, 0] = (-np.sin(sqrtkl) * sqrtk).real
            m[2, 2] = m[3, 3] = (np.cosh(sqrtkl)).real
            m[2, 3] = (np.sinh(sqrtkl) / sqrtk).real
            m[3, 2] = (-np.sinh(sqrtkl) * sqrtk).real
            m[4, 5] = float(length) / gamma / gamma
    else:
        print("warning: 'gamma' should be a positive float number.")
    return m
91393
function _setupInlineCrypt()
    {
        $lambda_functions =& self::_getLambdaFunctions();

        // We create max. 10 hi-optimized code for memory reason. Means: For each $key one ultra fast inline-crypt function.
        // (Currently, for Blowfish, one generated $lambda_function cost on php5.5@32bit ~100kb unfreeable mem and ~180kb on php5.5@64bit)
        // After that, we'll still create very fast optimized code but not the hi-ultimative code, for each $mode one.
        $gen_hi_opt_code = (bool)(count($lambda_functions) < 10);

        // Generation of a unique hash for our generated code
        $code_hash = "Crypt_Blowfish, {$this->mode}";
        if ($gen_hi_opt_code) {
            $code_hash = str_pad($code_hash, 32) . $this->_hashInlineCryptFunction($this->key);
        }

        $safeint = $this->safe_intval_inline();

        if (!isset($lambda_functions[$code_hash])) {
            switch (true) {
                case $gen_hi_opt_code:
                    $p = $this->bctx['p'];
                    $init_crypt = '
                        static $sb_0, $sb_1, $sb_2, $sb_3;
                        if (!$sb_0) {
                            $sb_0 = $self->bctx["sb"][0];
                            $sb_1 = $self->bctx["sb"][1];
                            $sb_2 = $self->bctx["sb"][2];
                            $sb_3 = $self->bctx["sb"][3];
                        }
                    ';
                    break;
                default:
                    $p   = array();
                    for ($i = 0; $i < 18; ++$i) {
                        $p[] = '$p_' . $i;
                    }
                    $init_crypt = '
                        list($sb_0, $sb_1, $sb_2, $sb_3) = $self->bctx["sb"];
                        list(' . implode(',', $p) . ') = $self->bctx["p"];

                    ';
            }

            // Generating encrypt code:
            $encrypt_block = '
                $in = unpack("N*", $in);
                $l = $in[1];
                $r = $in[2];
            ';
            for ($i = 0; $i < 16; $i+= 2) {
                $encrypt_block.= '
                    $l^= ' . $p[$i] . ';
                    $r^= ' . sprintf($safeint, '(' . sprintf($safeint, '$sb_0[$l >> 24 & 0xff] + $sb_1[$l >> 16 & 0xff]') . ' ^
                          $sb_2[$l >>  8 & 0xff]) +
                          $sb_3[$l       & 0xff]') . ';

                    $r^= ' . $p[$i + 1] . ';
                    $l^= ' . sprintf($safeint, '(' . sprintf($safeint, '$sb_0[$r >> 24 & 0xff] + $sb_1[$r >> 16 & 0xff]') . '  ^
                          $sb_2[$r >>  8 & 0xff]) +
                          $sb_3[$r       & 0xff]') . ';
                ';
            }
            $encrypt_block.= '
                $in = pack("N*",
                    $r ^ ' . $p[17] . ',
                    $l ^ ' . $p[16] . '
                );
            ';

            // Generating decrypt code:
            $decrypt_block = '
                $in = unpack("N*", $in);
                $l = $in[1];
                $r = $in[2];
            ';

            for ($i = 17; $i > 2; $i-= 2) {
                $decrypt_block.= '
                    $l^= ' . $p[$i] . ';
                    $r^= ' . sprintf($safeint, '(' . sprintf($safeint, '$sb_0[$l >> 24 & 0xff] + $sb_1[$l >> 16 & 0xff]') . ' ^
                          $sb_2[$l >>  8 & 0xff]) +
                          $sb_3[$l       & 0xff]') . ';

                    $r^= ' . $p[$i - 1] . ';
                    $l^= ' . sprintf($safeint, '(' . sprintf($safeint, '$sb_0[$r >> 24 & 0xff] + $sb_1[$r >> 16 & 0xff]') . ' ^
                          $sb_2[$r >>  8 & 0xff]) +
                          $sb_3[$r       & 0xff]') . ';
                ';
            }

            $decrypt_block.= '
                $in = pack("N*",
                    $r ^ ' . $p[0] . ',
                    $l ^ ' . $p[1] . '
                );
            ';

            $lambda_functions[$code_hash] = $this->_createInlineCryptFunction(
                array(
                   'init_crypt'    => $init_crypt,
                   'init_encrypt'  => '',
                   'init_decrypt'  => '',
                   'encrypt_block' => $encrypt_block,
                   'decrypt_block' => $decrypt_block
                )
            );
        }
        $this->inline_crypt = $lambda_functions[$code_hash];
    }
1303980
function ldexp( frac, exp ) {
	var high;
	var tmp;
	var w;
	var m;
	if (
		frac === 0 || // handles +-0
		frac !== frac || // handles NaN
		frac === PINF ||
		frac === NINF
	) {
		return frac;
	}
	// Normalize the input fraction:
	tmp = normalize( frac );
	frac = tmp[ 0 ];
	exp += tmp[ 1 ];

	// Extract the exponent from `frac` and add it to `exp`:
	exp += floatExp( frac );

	// Check for underflow/overflow...
	if ( exp < MIN_SUBNORMAL_EXPONENT ) {
		return copysign( 0, frac );
	}
	if ( exp > MAX_EXPONENT ) {
		if ( frac < 0 ) {
			return NINF;
		}
		return PINF;
	}
	// Check for a subnormal and scale accordingly to retain precision...
	if ( exp <= MAX_SUBNORMAL_EXPONENT ) {
		exp += 52;
		m = TWO52_INV;
	} else {
		m = 1.0;
	}
	// Split the fraction into higher and lower order words:
	w = toWords( frac );
	high = w[ 0 ];

	// Clear the exponent bits within the higher order word:
	high &= CLEAR_EXP_MASK;

	// Set the exponent bits to the new exponent:
	high |= ((exp+BIAS) << 20);

	// Create a new floating-point number:
	return m * fromWords( high, w[ 1 ] );
}
1363852
def Eg(self, **kwargs):
        '''
        
        '''
        return min(self.Eg_Gamma(**kwargs),
                   self.Eg_L(**kwargs),
                   self.Eg_X(**kwargs))
358759
def generate(self, x, **kwargs):
    
    # Parse and save attack-specific parameters
    assert self.parse_params(**kwargs)

    if self.symbolic_impl:
      # Create random targets if y_target not provided
      if self.y_target is None:
        from random import randint

        def random_targets(gt):
          result = gt.copy()
          nb_s = gt.shape[0]
          nb_classes = gt.shape[1]

          for i in range(nb_s):
            result[i, :] = np.roll(result[i, :],
                                   randint(1, nb_classes - 1))

          return result

        labels, nb_classes = self.get_or_guess_labels(x, kwargs)
        self.y_target = tf.py_func(random_targets, [labels],
                                   self.tf_dtype)
        self.y_target.set_shape([None, nb_classes])

      x_adv = jsma_symbolic(
          x,
          model=self.model,
          y_target=self.y_target,
          theta=self.theta,
          gamma=self.gamma,
          clip_min=self.clip_min,
          clip_max=self.clip_max)
    else:
      raise NotImplementedError("The jsma_batch function has been removed."
                                " The symbolic_impl argument to SaliencyMapMethod will be removed"
                                " on 2019-07-18 or after. Any code that depends on the non-symbolic"
                                " implementation of the JSMA should be revised. Consider using"
                                " SaliencyMapMethod.generate_np() instead.")

    return x_adv
352148
function qrPolynomial(num, shift) {

		if (typeof num.length == 'undefined') {
			throw new Error(num.length + '/' + shift);
		}

		var _num = function() {
			var offset = 0;
			while (offset < num.length && num[offset] == 0) {
				offset += 1;
			}
			var _num = new Array(num.length - offset + shift);
			for (var i = 0; i < num.length - offset; i += 1) {
				_num[i] = num[i + offset];
			}
			return _num;
		}();

		var _this = {};

		_this.getAt = function(index) {
			return _num[index];
		};

		_this.getLength = function() {
			return _num.length;
		};

		_this.multiply = function(e) {

			var num = new Array(_this.getLength() + e.getLength() - 1);

			for (var i = 0; i < _this.getLength(); i += 1) {
				for (var j = 0; j < e.getLength(); j += 1) {
					num[i + j] ^= QRMath.gexp(QRMath.glog(_this.getAt(i) ) + QRMath.glog(e.getAt(j) ) );
				}
			}

			return qrPolynomial(num, 0);
		};

		_this.mod = function(e) {

			if (_this.getLength() - e.getLength() < 0) {
				return _this;
			}

			var ratio = QRMath.glog(_this.getAt(0) ) - QRMath.glog(e.getAt(0) );

			var num = new Array(_this.getLength() );
			for (var i = 0; i < _this.getLength(); i += 1) {
				num[i] = _this.getAt(i);
			}

			for (var i = 0; i < e.getLength(); i += 1) {
				num[i] ^= QRMath.gexp(QRMath.glog(e.getAt(i) ) + ratio);
			}

			// recursive call
			return qrPolynomial(num, 0).mod(e);
		};

		return _this;
	}
1370914
def lpgamma(p, a):
    
    res = p * (p - 1.0) / 4.0 * log(pi)
    for ii in range(1, p+1):
        res += gammaln(a + (1.0 - ii) / 2.0)

    return res
1233893
def afterglow(self, src=None, event=None, dst=None, **kargs):
        """"""
        if src is None:
            src = lambda x: x['IP'].src
        if event is None:
            event = lambda x: x['IP'].dport
        if dst is None:
            dst = lambda x: x['IP'].dst
        sl = {}
        el = {}
        dl = {}
        for i in self.res:
            try:
                s,e,d = src(i),event(i),dst(i)
                if s in sl:
                    n,l = sl[s]
                    n += 1
                    if e not in l:
                        l.append(e)
                    sl[s] = (n,l)
                else:
                    sl[s] = (1,[e])
                if e in el:
                    n,l = el[e]
                    n+=1
                    if d not in l:
                        l.append(d)
                    el[e] = (n,l)
                else:
                    el[e] = (1,[d])
                dl[d] = dl.get(d,0)+1
            except:
                continue

        import math
        def normalize(n):
            return 2+math.log(n)/4.0

        def minmax(x):
            m,M = min(x),max(x)
            if m == M:
                m = 0
            if M == 0:
                M = 1
            return m,M

        #mins,maxs = minmax(map(lambda (x,y): x, sl.values()))
        mins,maxs = minmax([ a[0] for a in sl.values()])
        #mine,maxe = minmax(map(lambda (x,y): x, el.values()))
        mine,maxe = minmax([ a[0] for a in el.values()])
        mind,maxd = minmax(dl.values())
    
        gr = 'digraph "afterglow" {\n\tedge [len=2.5];\n'

        gr += "# src nodes\n"
        for s in sl:
            n,l = sl[s]; n = 1+(n-mins)/(maxs-mins)
            gr += '"src.%s" [label = "%s", shape=box, fillcolor="#FF0000", style=filled, fixedsize=1, height=%.2f,width=%.2f];\n' % (repr(s),repr(s),n,n)
        gr += "# event nodes\n"
        for e in el:
            n,l = el[e]; n = n = 1+(n-mine)/(maxe-mine)
            gr += '"evt.%s" [label = "%s", shape=circle, fillcolor="#00FFFF", style=filled, fixedsize=1, height=%.2f, width=%.2f];\n' % (repr(e),repr(e),n,n)
        for d in dl:
            n = dl[d]; n = n = 1+(n-mind)/(maxd-mind)
            gr += '"dst.%s" [label = "%s", shape=triangle, fillcolor="#0000ff", style=filled, fixedsize=1, height=%.2f, width=%.2f];\n' % (repr(d),repr(d),n,n)

        gr += "###\n"
        for s in sl:
            n,l = sl[s]
            for e in l:
                gr += ' "src.%s" -> "evt.%s";\n' % (repr(s),repr(e)) 
        for e in el:
            n,l = el[e]
            for d in l:
                gr += ' "evt.%s" -> "dst.%s";\n' % (repr(e),repr(d)) 
            
        gr += "}"
        return do_graph(gr, **kargs)