private static DoubleMatrix2D getGt(final DenseDoubleMatrix2D p, final DenseDoubleMatrix2D q, double lambda) {
    final int K = p.columns();

    DenseDoubleMatrix2D A1 = new DenseDoubleMatrix2D(K, K);
    q.zMult(q, A1, 1.0, 0.0, true, false);
    for (int k = 0; k < K; k++) {
        A1.setQuick(k, k, lambda + A1.getQuick(k, k));
    }

    EigenvalueDecomposition eig = new EigenvalueDecomposition(A1);
    DoubleMatrix1D d = eig.getRealEigenvalues();
    DoubleMatrix2D gt = eig.getV();
    for (int k = 0; k < K; k++) {
        double a = sqrt(d.get(k));
        gt.viewColumn(k).assign(x -> a * x);
    }

    return gt;
}
 

public double[] getCartesianCoords(double[] polars) {
	if (polars.length==0) { return new double[0]; }
	
	EigenvalueDecomposition evd = new EigenvalueDecomposition(A);
	DoubleMatrix2D Q = evd.getV();
	DoubleMatrix2D L = evd.getD();
	DoubleMatrix2D qT = Q.viewDice().copy();
	Algebra alg = new Algebra();		
	
	int n = polars.length;		
	double radius = polars[0];
	double[] phis = new double[n-1];
	for (int i=1; i<n; i++) { phis[i-1] = polars[i]; }
	
	double[] cartCoords = new double[n];
	for (int i=0; i<n; i++) {
		double prod = 1;
		for (int j=0; j<i; j++) { prod = prod * Math.sin(phis[j]); }
		if (i<n-1) { prod = prod * Math.cos(phis[i]); } 			
		cartCoords[i] = radius*prod;
	}
	
	// transform cartesian coordinates back!
	
	double[] u = new double[cartCoords.length];
	for (int i=0; i<u.length; i++) {
		u[i] = cartCoords[i]*Math.sqrt(L.get(i, i));
	}
	double[] s = alg.mult(Q, DoubleFactory1D.dense.make(u)).toArray();
	double[] x = new double[s.length];
	for (int i=0; i<x.length; i++) { x[i] = s[i] + c.get(i); }
	
	return x;
}
 

public EPolyPC( EPoly template, int fullOrder, int PCOrder, double PCFac) {
    //set up the principal component basis and the DOF information
    //coeffs will be added later
    //since we will be using this EPolyPC object when performing the fitting
    
    super(template.numDOFs,template.DOFmax,template.DOFmin,template.center,
            template.minE,null,fullOrder,template.DOFNames);
    
    //get the coordinate transformation into the eigenbasis of the template's Hessian
    //so we can use the high-eigenvalue directions as our principcal components
    DoubleMatrix2D hess = SeriesFitter.getHessian(template.coeffs,numDOFs,false);
    EigenvalueDecomposition edec = new EigenvalueDecomposition(hess);

    DoubleMatrix1D eigVals = edec.getRealEigenvalues();
    DoubleMatrix2D eigVecs = edec.getV();//The columns of eigVec are the eigenvectors

    invAxisCoeffs = eigVecs;//axisCoeffs' inverse is its transpose, since it's orthogonal
    axisCoeffs = Algebra.DEFAULT.transpose(eigVecs);
    

    this.fullOrder = fullOrder;
    this.PCOrder = PCOrder;


    //Now figure out what PC's to use
    //Let's try the criterion to use all PCs
    //whose eigenvalues' absolute values are within a factor of PCFac
    //of the biggest

    double maxEigVal = 0;
    for(int n=0; n<numDOFs; n++)
        maxEigVal = Math.max( Math.abs(eigVals.get(n)), maxEigVal );

    isPC = new boolean[numDOFs];

    for(int n=0; n<numDOFs; n++){
        if( Math.abs(eigVals.get(n)) >= maxEigVal*PCFac )
            isPC[n] = true;
    }
}
 

public double[] getEllipsoidalCoords(double[] dihedrals) {
    	
        if (dihedrals.length==0) { return new double[0]; }
    	
    	// for now we're just using the unit sphere
    	DoubleMatrix2D A = DoubleFactory2D.dense.identity(dihedrals.length);
    	DoubleMatrix1D c = DoubleFactory1D.dense.make(new double[dihedrals.length]);
    	
        EigenvalueDecomposition evd = new EigenvalueDecomposition(A);
        DoubleMatrix2D Q = evd.getV();
        DoubleMatrix2D L = evd.getD();
        DoubleMatrix2D qT = Q.viewDice().copy();
        Algebra alg = new Algebra();
    	
    	// first transform the cartesian coordinates based on the ellipse
    	double[] s = new double[dihedrals.length];
    	for (int i=0; i<dihedrals.length; i++) { 
    		s[i] = dihedrals[i]-c.get(i);
    	}
    	double[] u = alg.mult(qT, DoubleFactory1D.dense.make(s)).toArray();
    	double[] x = new double[u.length];
    	for (int i=0; i<u.length; i++) {
    		x[i] = u[i]/Math.sqrt(L.get(i, i));
    	}    	
    	dihedrals = x;
    	
    	// now get elliptical coordinates
/*    	double radius = 0;
    	for (double d : dihedrals) { radius += d*d; }
    	radius = Math.sqrt(radius);*/
    	int n = dihedrals.length;
    	double[] phi = new double[n-1];
    	for (int i=0; i<n-1; i++) {
    		double d=0;
    		for (int j=i; j<n; j++) { d += dihedrals[j]*dihedrals[j]; }
    		double quot = dihedrals[i]/Math.sqrt(d);
    		phi[i] = Math.acos(quot);
    	}
    	if (dihedrals[n-1] < 0) { phi[n-2] = 2*Math.PI - phi[n-2]; }
    	double[] ellCoords = new double[n];
    	ellCoords[0] = 0; //radius;
    	for (int i=1; i<n; i++) { ellCoords[i] = phi[i-1]; }
    	return ellCoords;
    }
 

static DoubleMatrix2D invertX(DoubleMatrix2D X){
        
        if(X.rows()!=X.columns())
            throw new RuntimeException("ERROR: Can't invert square matrix");
        
        if(X.rows()==1){//easy case of 1-D
            if(Math.abs(X.get(0,0))<1e-8)
                return DoubleFactory2D.dense.make(1,1,0);
            else
                return DoubleFactory2D.dense.make(1,1,1./X.get(0,0));
        }
    
        /*DoubleMatrix2D invX = null;
    
        try{
            invX = Algebra.DEFAULT.inverse(X);
        }
        catch(Exception e){
            System.err.println("ERROR: Singular matrix in MinVolEllipse calculation");
        }*/
        
        //Invert X, but if singular, this probably indicates a lack of points
        //but we can still make a meaningful lower-dimensional ellipsoid for the dimensions we have data
        //so we construct a pseudoinverse by setting those eigencomponents to 0
        //X is assumed to be symmetric (if not would need SVD instead...)
    
        EigenvalueDecomposition edec = new EigenvalueDecomposition(X);
        
        DoubleMatrix2D newD = edec.getD().copy();
        for(int m=0; m<X.rows(); m++){
            if( Math.abs(newD.get(m,m)) < 1e-8 ){
                //System.out.println("Warning: singular X in MinVolEllipse calculation");
                //this may come up somewhat routinely, especially with discrete DOFs in place
                newD.set(m,m,0);
            }
            else
                newD.set(m,m,1./newD.get(m,m));
        }

        DoubleMatrix2D invX = Algebra.DEFAULT.mult(edec.getV(),newD).zMult(edec.getV(), null, 1, 0, false, true);
        
        return invX;
}