Re-order the eigenvalues according to the similarity of the eigenvectors
169217
def band_reorder(self):
        
        eiv = self.eigendisplacements
        eig = self.bands
    
        nphonons,nqpoints = self.bands.shape
        order = np.zeros([nqpoints,nphonons],dtype=int)
        order[0] = np.array(range(nphonons))
        
        #get the atomic masses
        atomic_masses = [ site.specie.atomic_mass for site in self.structure.sites ]
 
        #get order
        for nq in range(1,nqpoints):
            old_eiv = eigenvectors_from_displacements(eiv[:,nq-1],atomic_masses)
            new_eiv = eigenvectors_from_displacements(eiv[:,nq],  atomic_masses)
            order[nq] = estimate_band_connection(old_eiv.reshape([nphonons,nphonons]).T,
                                                 new_eiv.reshape([nphonons,nphonons]).T,
                                                 order[nq-1])

        #reorder
        for nq in range(1,nqpoints):
            eivq=eiv[:,nq]
            eigq=eig[:,nq]
            eiv[:,nq] = eivq[order[nq]]
            eig[:,nq] = eigq[order[nq]]
563087
def vector_similarity(self, vector, items):
        """"""
        vector = self.normalize(vector)
        items_vec = np.stack([self.norm_vectors[self.items[x]] for x in items])
        return vector.dot(items_vec.T)
1108755
def cosine_similarity_vec(num_tokens, num_removed_vec):
    
    remaining = -np.array(num_removed_vec) + num_tokens
    return remaining / (np.sqrt(num_tokens + 1e-6) * np.sqrt(remaining + 1e-6))
408799
def tf_min_eig_vec(self):
    """"""
    # Full eigen decomposition requires the explicit psd matrix M
    _, matrix_m = self.dual_object.get_full_psd_matrix()
    [eig_vals, eig_vectors] = tf.self_adjoint_eig(matrix_m)
    index = tf.argmin(eig_vals)
    return tf.reshape(
        eig_vectors[:, index], shape=[eig_vectors.shape[0].value, 1])
109285
public void rearrangeTriangles(int[] indices) {
        assert indices.length % 3 == 0;
        // Step 1: Make sure that the first index of each triangle is the smallest
        // one (rotate triangle nodes if necessary)
        for (int off = 0; off < indices.length; off += 3) {
            if ((indices[off + 1] < indices[off]) && (indices[off + 1] < indices[off + 2])) {
                int tmp = indices[off];
                indices[off] = indices[off + 1];
                indices[off + 1] = indices[off + 2];
                indices[off + 2] = tmp;
            } else if ((indices[off + 2] < indices[off]) && (indices[off + 2] < indices[off + 1])) {
                int tmp = indices[off];
                indices[off] = indices[off + 2];
                indices[off + 2] = indices[off + 1];
                indices[off + 1] = tmp;
            }
        }

        // Step 2: Sort the triangles based on the first triangle index
        Triangle[] tris = new Triangle[indices.length / 3];
        for (int i = 0; i < tris.length; i++) {
            int off = i * 3;
            tris[i] = new Triangle(indices, off);
        }

        Arrays.sort(tris);

        for (int i = 0; i < tris.length; i++) {
            int off = i * 3;
            tris[i].copyBack(indices, off);
        }
    }
946169
def eigendecompose(self, normalise=False):
        
        self.eigenvalues, self.eigenvectors = utils.eigendecompose(self.tensor,
                                                                   normalise)
        return self.eigenvalues, self.eigenvectors
1350223
def _symmetrize_correlograms(correlograms):
    """"""

    n_clusters, _, n_bins = correlograms.shape
    assert n_clusters == _

    # We symmetrize c[i, j, 0].
    # This is necessary because the algorithm in correlograms()
    # is sensitive to the order of identical spikes.
    correlograms[..., 0] = np.maximum(correlograms[..., 0],
                                      correlograms[..., 0].T)

    sym = correlograms[..., 1:][..., ::-1]
    sym = np.transpose(sym, (1, 0, 2))

    return np.dstack((sym, correlograms))
104809
def eigenvectors_right_samples(self):
         """
        res = np.empty((self.nsamples, self.nstates, self.nstates), dtype=config.dtype)
        for i in range(self.nsamples):
            res[i, :, :] = self._sampled_hmms[i].eigenvectors_right
        return res
400752
def spectral_embedding(self, n):
        
        coords = spectral_embedding(self._affinity, n)
        return CoordinateMatrix(normalise_rows(coords))
902855
def find_best_rsquared(list_of_fits):
        """"""
        res = sorted(list_of_fits, key=lambda x: x.rsquared)
        return res[-1]
188936
def fitlin(imgarr,refarr):
    
    # Initialize variables
    _mat = np.zeros((3,3),dtype=np.float64)
    _xorg = imgarr[0][0]
    _yorg = imgarr[0][1]
    _xoorg = refarr[0][0]
    _yoorg = refarr[0][1]
    _sigxox = 0.
    _sigxoy = 0.
    _sigxo = 0.
    _sigyox = 0.
    _sigyoy = 0.
    _sigyo = 0.

    _npos = len(imgarr)
    # Populate matrices
    for i in range(_npos):
        _mat[0][0] += np.power((imgarr[i][0] - _xorg),2)
        _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg)
        _mat[0][2] += (imgarr[i][0] - _xorg)
        _mat[1][1] += np.power((imgarr[i][1] - _yorg),2)
        _mat[1][2] += imgarr[i][1] - _yorg

        _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg)
        _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg)
        _sigxo += refarr[i][0] - _xoorg
        _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg)
        _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg)
        _sigyo += refarr[i][1] - _yoorg

    _mat[2][2] = _npos
    _mat[1][0] = _mat[0][1]
    _mat[2][0] = _mat[0][2]
    _mat[2][1] = _mat[1][2]

    # Now invert this matrix
    _mat = linalg.inv(_mat)

    _a  = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2]
    _b  = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2])
    #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2]

    _c  = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2]
    _d  = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2]
    #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2]

    _xt = _xoorg - _a*_xorg+_b*_yorg
    _yt = _yoorg - _d*_xorg-_c*_yorg

    return [_a,_b,_xt],[_c,_d,_yt]
767896
def eigenvalue_band_properties(self):
        
        vbm = -float("inf")
        vbm_kpoint = None
        cbm = float("inf")
        cbm_kpoint = None
        for spin, d in self.eigenvalues.items():
            for k, val in enumerate(d):
                for (eigenval, occu) in val:
                    if occu > self.occu_tol and eigenval > vbm:
                        vbm = eigenval
                        vbm_kpoint = k
                    elif occu <= self.occu_tol and eigenval < cbm:
                        cbm = eigenval
                        cbm_kpoint = k
        return max(cbm - vbm, 0), cbm, vbm, vbm_kpoint == cbm_kpoint
865013
def higher_order_diff(eqs, syms, order=2):
    ''''''

    import numpy

    eqs = list([sympy.sympify(eq) for eq in eqs])
    syms = list([sympy.sympify(s) for s in syms])

    neq = len(eqs)
    p = len(syms)

    D = [numpy.array(eqs)]
    orders = []

    for i in range(1,order+1):

        par = D[i-1]
        mat = numpy.empty([neq] + [p]*i, dtype=object)   #.append( numpy.zeros(orders))
        for ind in non_decreasing_series(p,i):

            ind_parent = ind[:-1]
            k = ind[-1]

            for line in range(neq):

                ii = [line] + ind
                iid = [line] + ind_parent
                eeq = par[ tuple(iid) ]
                mat[tuple(ii)] = eeq.diff(syms[k])

        D.append(mat)

    return D
129561
def eigen_table(self):
        """"""
        idx = ["Eigenvalue", "Variability (%)", "Cumulative (%)"]
        table = pd.DataFrame(
            np.array(
                [self.eigenvalues, self.inertia, self.cumulative_inertia]
            ),
            columns=["F%s" % i for i in range(1, self.keep + 1)],
            index=idx,
        )

        return table
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;
            }
        }
    }
695968
def PCA(Y, components):
	

	N,D = Y.shape
	sv = linalg.svd(Y, full_matrices=0);
	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
	v = s0.std(axis=0)
	s0 /= v;
	w0 *= v;
	return [s0, w0]

	if N>D:
		sv = linalg.svd(Y, full_matrices=0);
		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
		v = s0.std(axis=0)
		s0 /= v;
		w0 *= v;
		return [s0, w0]
	else:
		K=np.cov(Y)
		sv = linalg.eigh(K)
		std_var = np.sqrt(sv[0])
		pc = sv[1]*std_var[np.newaxis(),0]
		#import ipdb
		#ipdb.set_trace()
		return [pc,std_var]
508465
@Override
    public FloatMatrix repelem(final int rowRepeats, final int colRepeats) {
        N.checkArgument(rowRepeats > 0 && colRepeats > 0, "rowRepeats=%s and colRepeats=%s must be bigger than 0", rowRepeats, colRepeats);

        final float[][] c = new float[rows * rowRepeats][cols * colRepeats];

        for (int i = 0; i < rows; i++) {
            final float[] fr = c[i * rowRepeats];

            for (int j = 0; j < cols; j++) {
                N.copy(Array.repeat(a[i][j], colRepeats), 0, fr, j * colRepeats, colRepeats);
            }

            for (int k = 1; k < rowRepeats; k++) {
                N.copy(fr, 0, c[i * rowRepeats + k], 0, fr.length);
            }
        }

        return new FloatMatrix(c);
    }