Chebyshev basis functions T_n.
961743
def basis(self, n):
        
        if n == 0:
            return self(np.array([1.]))
        vals = np.ones(n+1)
        vals[1::2] = -1
        return self(vals)
1242341
def sandwich_decompositions(coords0, coords1, samples=SAMPLES):
    """"""
    decomps = []
    for _ in range(samples):
        circ = qf.Circuit()
        circ += qf.CANONICAL(*coords0, 0, 1)
        circ += qf.random_gate([0])
        circ += qf.random_gate([1])
        circ += qf.CANONICAL(*coords1, 0, 1)
        gate = circ.asgate()

        coords = qf.canonical_coords(gate)
        decomps.append(coords)

    return decomps
1282771
def xyz2lab(x__, y__, z__):
    

    def f__(arr):
        """f__ function
        """

        return np.where(arr > 216.0 / 24389.0,
                        arr ** (1.0/3.0),
                        (1.0 / 3.0) * (29.0 / 6.0) * (29.0 / 6.0) * arr
                        + 4.0 / 29.0)
    fy_ = f__(y__ / 100.0)
    l__ = 116 * fy_ - 16
    a__ = 500.0 * (f__(x__ / 95.047) - fy_)
    b__ = 200.0 * (fy_ - f__(z__ / 108.883))
    return l__, a__, b__
348905
def stSpectralEntropy(X, n_short_blocks=10):
    """"""
    L = len(X)                         # number of frame samples
    Eol = numpy.sum(X ** 2)            # total spectral energy

    sub_win_len = int(numpy.floor(L / n_short_blocks))   # length of sub-frame
    if L != sub_win_len * n_short_blocks:
        X = X[0:sub_win_len * n_short_blocks]

    sub_wins = X.reshape(sub_win_len, n_short_blocks, order='F').copy()  # define sub-frames (using matrix reshape)
    s = numpy.sum(sub_wins ** 2, axis=0) / (Eol + eps)                      # compute spectral sub-energies
    En = -numpy.sum(s*numpy.log2(s + eps))                                    # compute spectral entropy

    return En
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
85829
def function(self, x, y, amp, R_sersic, n_sersic, center_x=0, center_y=0):
        
        #if n_sersic < 0.2:
        #    n_sersic = 0.2
        #if R_sersic < 10.**(-6):
        #    R_sersic = 10.**(-6)
        R_sersic = np.maximum(0, R_sersic)
        x_shift = x - center_x
        y_shift = y - center_y
        R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
        if isinstance(R, int) or isinstance(R, float):
            R = max(self._smoothing, R)
        else:
            R[R < self._smoothing] = self._smoothing
        _, bn = self.k_bn(n_sersic, R_sersic)
        R_frac = R/R_sersic
        #R_frac = R_frac.astype(np.float32)
        if isinstance(R, int) or isinstance(R, float):
            if R_frac > 100:
                result = 0
            else:
                exponent = -bn*(R_frac**(1./n_sersic)-1.)
                result = amp * np.exp(exponent)
        else:
            R_frac_real = R_frac[R_frac <= 100]
            exponent = -bn*(R_frac_real**(1./n_sersic)-1.)
            result = np.zeros_like(R)
            result[R_frac <= 100] = amp * np.exp(exponent)
        return np.nan_to_num(result)
1092191
private function hqr2()
    {
        //  Initialize
        $nn = $this->n;
        $n  = $nn - 1;
        $low = 0;
        $high = $nn - 1;
        $eps = pow(2.0, -52.0);
        $exshift = 0.0;
        $p = $q = $r = $s = $z = 0;
        // Store roots isolated by balanc and compute matrix norm
        $norm = 0.0;

        for ($i = 0; $i < $nn; ++$i) {
            if (($i < $low) or ($i > $high)) {
                $this->d[$i] = $this->H[$i][$i];
                $this->e[$i] = 0.0;
            }
            for ($j = max($i-1, 0); $j < $nn; ++$j) {
                $norm = $norm + abs($this->H[$i][$j]);
            }
        }

        // Outer loop over eigenvalue index
        $iter = 0;
        while ($n >= $low) {
            // Look for single small sub-diagonal element
            $l = $n;
            while ($l > $low) {
                $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
                if ($s == 0.0) {
                    $s = $norm;
                }
                if (abs($this->H[$l][$l-1]) < $eps * $s) {
                    break;
                }
                --$l;
            }
            // Check for convergence
            // One root found
            if ($l == $n) {
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->d[$n] = $this->H[$n][$n];
                $this->e[$n] = 0.0;
                --$n;
                $iter = 0;
            // Two roots found
            } elseif ($l == $n-1) {
                $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
                $q = $p * $p + $w;
                $z = sqrt(abs($q));
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
                $x = $this->H[$n][$n];
                // Real pair
                if ($q >= 0) {
                    if ($p >= 0) {
                        $z = $p + $z;
                    } else {
                        $z = $p - $z;
                    }
                    $this->d[$n-1] = $x + $z;
                    $this->d[$n] = $this->d[$n-1];
                    if ($z != 0.0) {
                        $this->d[$n] = $x - $w / $z;
                    }
                    $this->e[$n-1] = 0.0;
                    $this->e[$n] = 0.0;
                    $x = $this->H[$n][$n-1];
                    $s = abs($x) + abs($z);
                    $p = $x / $s;
                    $q = $z / $s;
                    $r = sqrt($p * $p + $q * $q);
                    $p = $p / $r;
                    $q = $q / $r;
                    // Row modification
                    for ($j = $n-1; $j < $nn; ++$j) {
                        $z = $this->H[$n-1][$j];
                        $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
                        $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
                    }
                    // Column modification
                    for ($i = 0; $i <= $n; ++$i) {
                        $z = $this->H[$i][$n-1];
                        $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
                        $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
                    }
                    // Accumulate transformations
                    for ($i = $low; $i <= $high; ++$i) {
                        $z = $this->V[$i][$n-1];
                        $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
                        $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
                    }
                // Complex pair
                } else {
                    $this->d[$n-1] = $x + $p;
                    $this->d[$n] = $x + $p;
                    $this->e[$n-1] = $z;
                    $this->e[$n] = -$z;
                }
                $n = $n - 2;
                $iter = 0;
            // No convergence yet
            } else {
                // Form shift
                $x = $this->H[$n][$n];
                $y = 0.0;
                $w = 0.0;
                if ($l < $n) {
                    $y = $this->H[$n-1][$n-1];
                    $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                }
                // Wilkinson's original ad hoc shift
                if ($iter == 10) {
                    $exshift += $x;
                    for ($i = $low; $i <= $n; ++$i) {
                        $this->H[$i][$i] -= $x;
                    }
                    $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
                    $x = $y = 0.75 * $s;
                    $w = -0.4375 * $s * $s;
                }
                // MATLAB's new ad hoc shift
                if ($iter == 30) {
                    $s = ($y - $x) / 2.0;
                    $s = $s * $s + $w;
                    if ($s > 0) {
                        $s = sqrt($s);
                        if ($y < $x) {
                            $s = -$s;
                        }
                        $s = $x - $w / (($y - $x) / 2.0 + $s);
                        for ($i = $low; $i <= $n; ++$i) {
                            $this->H[$i][$i] -= $s;
                        }
                        $exshift += $s;
                        $x = $y = $w = 0.964;
                    }
                }
                // Could check iteration count here.
                $iter = $iter + 1;
                // Look for two consecutive small sub-diagonal elements
                $m = $n - 2;
                while ($m >= $l) {
                    $z = $this->H[$m][$m];
                    $r = $x - $z;
                    $s = $y - $z;
                    $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
                    $q = $this->H[$m+1][$m+1] - $z - $r - $s;
                    $r = $this->H[$m+2][$m+1];
                    $s = abs($p) + abs($q) + abs($r);
                    $p = $p / $s;
                    $q = $q / $s;
                    $r = $r / $s;
                    if ($m == $l) {
                        break;
                    }
                    if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
                        $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
                        break;
                    }
                    --$m;
                }
                for ($i = $m + 2; $i <= $n; ++$i) {
                    $this->H[$i][$i-2] = 0.0;
                    if ($i > $m+2) {
                        $this->H[$i][$i-3] = 0.0;
                    }
                }
                // Double QR step involving rows l:n and columns m:n
                for ($k = $m; $k <= $n-1; ++$k) {
                    $notlast = ($k != $n-1);
                    if ($k != $m) {
                        $p = $this->H[$k][$k-1];
                        $q = $this->H[$k+1][$k-1];
                        $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
                        $x = abs($p) + abs($q) + abs($r);
                        if ($x != 0.0) {
                            $p = $p / $x;
                            $q = $q / $x;
                            $r = $r / $x;
                        }
                    }
                    if ($x == 0.0) {
                        break;
                    }
                    $s = sqrt($p * $p + $q * $q + $r * $r);
                    if ($p < 0) {
                        $s = -$s;
                    }
                    if ($s != 0) {
                        if ($k != $m) {
                            $this->H[$k][$k-1] = -$s * $x;
                        } elseif ($l != $m) {
                            $this->H[$k][$k-1] = -$this->H[$k][$k-1];
                        }
                        $p = $p + $s;
                        $x = $p / $s;
                        $y = $q / $s;
                        $z = $r / $s;
                        $q = $q / $p;
                        $r = $r / $p;
                        // Row modification
                        for ($j = $k; $j < $nn; ++$j) {
                            $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
                            if ($notlast) {
                                $p = $p + $r * $this->H[$k+2][$j];
                                $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
                            }
                            $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
                            $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
                        }
                        // Column modification
                        for ($i = 0; $i <= min($n, $k+3); ++$i) {
                            $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->H[$i][$k+2];
                                $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
                            }
                            $this->H[$i][$k] = $this->H[$i][$k] - $p;
                            $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
                        }
                        // Accumulate transformations
                        for ($i = $low; $i <= $high; ++$i) {
                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->V[$i][$k+2];
                                $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
                            }
                            $this->V[$i][$k] = $this->V[$i][$k] - $p;
                            $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
                        }
                    }  // ($s != 0)
                }  // k loop
            }  // check convergence
        }  // while ($n >= $low)

        // Backsubstitute to find vectors of upper triangular form
        if ($norm == 0.0) {
            return;
        }

        for ($n = $nn-1; $n >= 0; --$n) {
            $p = $this->d[$n];
            $q = $this->e[$n];
            // Real vector
            if ($q == 0) {
                $l = $n;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-1; $i >= 0; --$i) {
                    $w = $this->H[$i][$i] - $p;
                    $r = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $s = $r;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0.0) {
                            if ($w != 0.0) {
                                $this->H[$i][$n] = -$r / $w;
                            } else {
                                $this->H[$i][$n] = -$r / ($eps * $norm);
                            }
                        // Solve real equations
                        } else {
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
                            $t = ($x * $s - $z * $r) / $q;
                            $this->H[$i][$n] = $t;
                            if (abs($x) > abs($z)) {
                                $this->H[$i+1][$n] = (-$r - $w * $t) / $x;
                            } else {
                                $this->H[$i+1][$n] = (-$s - $y * $t) / $z;
                            }
                        }
                        // Overflow control
                        $t = abs($this->H[$i][$n]);
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n] = $this->H[$j][$n] / $t;
                            }
                        }
                    }
                }
            // Complex vector
            } elseif ($q < 0) {
                $l = $n-1;
                // Last vector component imaginary so matrix is triangular
                if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
                    $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
                    $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
                } else {
                    $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
                    $this->H[$n-1][$n-1] = $this->cdivr;
                    $this->H[$n-1][$n]   = $this->cdivi;
                }
                $this->H[$n][$n-1] = 0.0;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-2; $i >= 0; --$i) {
                    // double ra,sa,vr,vi;
                    $ra = 0.0;
                    $sa = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
                        $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    $w = $this->H[$i][$i] - $p;
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $r = $ra;
                        $s = $sa;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0) {
                            $this->cdiv(-$ra, -$sa, $w, $q);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                        } else {
                            // Solve complex equations
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
                            $vi = ($this->d[$i] - $p) * 2.0 * $q;
                            if ($vr == 0.0 & $vi == 0.0) {
                                $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
                            }
                            $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                            if (abs($x) > (abs($z) + abs($q))) {
                                $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
                                $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
                            } else {
                                $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
                                $this->H[$i+1][$n-1] = $this->cdivr;
                                $this->H[$i+1][$n]   = $this->cdivi;
                            }
                        }
                        // Overflow control
                        $t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n]));
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
                                $this->H[$j][$n]   = $this->H[$j][$n] / $t;
                            }
                        }
                    } // end else
                } // end for
            } // end else for complex case
        } // end for

        // Vectors of isolated roots
        for ($i = 0; $i < $nn; ++$i) {
            if ($i < $low | $i > $high) {
                for ($j = $i; $j < $nn; ++$j) {
                    $this->V[$i][$j] = $this->H[$i][$j];
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix
        for ($j = $nn-1; $j >= $low; --$j) {
            for ($i = $low; $i <= $high; ++$i) {
                $z = 0.0;
                for ($k = $low; $k <= min($j, $high); ++$k) {
                    $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
                }
                $this->V[$i][$j] = $z;
            }
        }
    }
805465
function (m, u) {
        return [
            m[0][0] * u[0] + m[0][1] * u[1] + m[0][2] * u[2],
            m[1][0] * u[0] + m[1][1] * u[1] + m[1][2] * u[2],
            m[2][0] * u[0] + m[2][1] * u[1] + m[2][2] * u[2]
        ];
    }
75233
function interpolatePairs( _ps, t1, t2 ) {
  var pn = [];
  for (var i=0; i<_ps.length; i++) {
    var next = (i==_ps.length-1) ? 0 : i+1;
    pn.push( new Pair( _ps[i].interpolate( t1 ) ).to( _ps[next].interpolate( 1-t2 ) ) );
  }
  return pn;
}
1271192
def __basis0(self, xi):
        """"""
        return np.where(np.all([self.knot_vector[:-1] <=  xi,
                                xi < self.knot_vector[1:]],axis=0), 1.0, 0.0)
1371562
function fft2d(re, im) {
    var tre = [],
        tim = [],
        i = 0;
    // x-axis
    for (var y = 0; y < _n; y++) {
        i = y * _n;
        for (var x1 = 0; x1 < _n; x1++) {
            tre[x1] = re[x1 + i];
            tim[x1] = im[x1 + i];
        }
        fft1d(tre, tim);
        for (var x2 = 0; x2 < _n; x2++) {
            re[x2 + i] = tre[x2];
            im[x2 + i] = tim[x2];
        }
    }
    // y-axis
    for (var x = 0; x < _n; x++) {
        for (var y1 = 0; y1 < _n; y1++) {
            i = x + y1 * _n;
            tre[y1] = re[i];
            tim[y1] = im[i];
        }
        fft1d(tre, tim);
        for (var y2 = 0; y2 < _n; y2++) {
            i = x + y2 * _n;
            re[i] = tre[y2];
            im[i] = tim[y2];
        }
    }
}
1096874
def bernstein(n, t):
    """"""
    t1 = 1-t
    return [n_choose_k(n, k) * t1**(n-k) * t**k for k in range(n+1)]
459067
def _Aterm(N, h, m, k, dW):
    """"""
    sqrt2h = np.sqrt(2.0/h)
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    term1 = _dot(Xk, _t(Yk + sqrt2h*dW))
    term2 = _dot(Yk + sqrt2h*dW, _t(Xk))
    return (term1 - term2)/k
40356
def fourier(x, N):
    """"""

    term = 0.
    for n in range(1, N, 2):
        term += (1. / n) * math.sin(n * math.pi * x / L)
    return (4. / (math.pi)) * term
915966
def sbesselj(x, N):
    

    nmax = N - 1;
    out = np.zeros(N, dtype=np.float64)
    z = x ** 2

    out[0] = np.sin(x) / x
    j1 = np.sin(x) / z - np.cos(x) / x

    u = 1
    v = x / (2.0 * nmax + 1.0)
    w = v
    n = nmax

    while(np.abs(v / w) > 1e-20):
        n = n + 1
        u = 1 / (1 - z * u / (4.0 * n ** 2 - 1.0))
        v *= u - 1
        w += v

    out[nmax] = w

    for n in xrange(nmax - 1, 0, -1):
        out[n] = 1.0 / ((2.0 * n + 1.0) / x - out[n + 1])

    if(np.abs(out[0]) >= np.abs(j1)):
        out[1] *= out[0]
    else:
        out[1] = j1

    for n in xrange(1, nmax):
        out[n + 1] *= out[n]

    return out
1351328
def bb2hw(a:Collection[int])->np.ndarray:
    ""
    return np.array([a[1],a[0],a[3]-a[1],a[2]-a[0]])
1314882
def fpbspl(t, n, k, x, l, h, hh):
    
    h[0] = 1.0
    for j in range(1, k + 1):
        hh[0:j] = h[0:j]
        h[0] = 0.0
        for i in range(j):
            li = l+i
            f = hh[i]/(t[li] - t[li - j])
            h[i] = h[i] + f*(t[li] - x)
            h[i + 1] = f*(x - t[li - j])