001package kernels;
002
003import java.util.Arrays;
004
005import org.apache.commons.math3.linear.BlockRealMatrix;
006import org.apache.commons.math3.linear.EigenDecomposition;
007import org.apache.commons.math3.linear.RealMatrix;
008import org.apache.commons.math3.linear.RealVector;
009
010/**
011 * An integration kernel \( K(\theta,\eta) \) corresponding to an operator on the one-particle space
012 * of an integrable theory. 
013 * 
014 * <p>This class is abstract and doesn't provide the actual kernel function;
015 * this is left to subclasses. Rather, this class provides conversion routines
016 * from the kernel function to matrix elements of the operator, and
017 * for approximately finding its spectral data.</p>
018 * 
019 * <p>To that end, the following orthonormal system of piecewise constant bumps is used as a basis:
020 * \[  \varphi_j (x) = \begin{cases}
021 *     h^{-1/2} \quad & \text{if }  -R + j h < x < -R+(j+1)h , \\
022 *     0 & \text{otherwise} ,
023 *  \end{cases} \]
024 *  
025 *  with \(j=0,\ldots,N-1\), where \(R>0\) is a given "cutoff", \(N\) the number of basis elements, and \(h=2R/N\).
026 *  All approximation results for the eigenvectors, the spectrum, etc. of the operator will depend on these
027 *  approximation parameters \(R\) and \(N\). </p>
028 * 
029 */
030public abstract class OneParticleKernel
031{
032    
033        private boolean evCheckEnabled = true;
034        
035    /**
036     * Evaluates the kernel \( K(\theta,\eta) \) at a given position \( \theta,\eta \).
037     *  
038     * @param theta first argument of the kernel
039     * @param eta second argument of the kernel
040     * @return value of the kernel at that position.
041     */
042    public abstract double kernelValue(double theta, double eta);
043
044    /**
045     * Disable the plausibility check of decaying eigenfunctions for this kernel.
046     * This should be called once, when for this specific kernel it is not expected that the eigenfunctions decay.
047     * The check performed by {@link #checkPlausibility} will then be disabled for all method calls. 
048     * 
049     * @see #findLowestEigenvector
050     * @see #findSpectrum
051     */
052    public void disableEVCheck() {
053        this.evCheckEnabled = false;
054    }
055    
056    /**
057     * Retrieve the matrix element of this kernel between two basis vectors. 
058     * 
059     *  The matrix elements \( (\varphi_i , A \varphi_j) \) are approximated using the midpoint rule.
060     *   
061     * @param i index of the first basis vector
062     * @param j index of the second basis vector
063     * @param n number of basis vectors
064     * @param cutoff the "cutoff" \(R\)
065     * @return the matrix element \( (\varphi_i , A \varphi_j) \)
066     */
067    protected double matrixElement(int i, int j, int n, double cutoff) {
068        double step = 2.0*cutoff/n;
069        double theta = -cutoff + step*(i + 0.5);
070        double eta = -cutoff + step*(j + 0.5);
071        double mtxe = step * kernelValue(theta, eta); 
072        return mtxe;
073    }
074    
075    
076    /**
077     * Compute a matrix approximation of this operator, in the basis described above.
078     *   
079     * @param n number of basis vectors
080     * @param cutoff the "cutoff" \(R\)
081     * @return a matrix approximation of this operator in the given basis
082     */
083    public RealMatrix asRealMatrix(int n, double cutoff) {
084        RealMatrix matrix = new BlockRealMatrix(n, n);
085        for (int i = 0; i < n; i++) {
086            for (int j = 0; j <= i; j++) {
087                double mtxe = matrixElement(i, j, n, cutoff);
088                matrix.setEntry(i, j, mtxe);
089                matrix.setEntry(j, i, mtxe);
090            }
091        }
092        return matrix;            
093    }
094
095    /**
096     * Compute an approximate eigen-decomposition of this operator, in the basis described above.
097     *   
098     * @param n number of basis vectors
099     * @param cutoff the "cutoff" \(R\)
100     * @return eigen-decomposition of this operator in the given basis
101     */
102    public EigenDecomposition findEigendec(int n, double cutoff) {
103        RealMatrix a = this.asRealMatrix(n, cutoff);
104        return new EigenDecomposition(a);
105    }
106
107    /**
108     * Check whether an eigen-decomposition of this operator is plausible, 
109     * namely, whether its lowest eigenfunction decays at the end of the interval \([-R,R]\).
110     * 
111     * This is useful for our main use case, when we want to find the eigenfunction for the lowest eigenvector
112     * independent of the cutoff \(R\).
113     * The method also checks whether the eigenvalues are all real.
114     * (We expect to deal only with symmetric operators.) 
115     * 
116     * @see #disableEVCheck()
117     * @see #findEigendec(int, double)
118     *   
119     * @param ed the eigen-decomposition to verify
120     * @param n number of basis vectors
121     * @throws RuntimeException if the matrix has complex eigenvalues
122     * @throws RuntimeException if the lowest eigenfunction does not decay
123     */
124    protected void checkPlausibility(EigenDecomposition ed, int n) {
125        if (ed.hasComplexEigenvalues()) {
126                throw new RuntimeException("complex eigenvalues detected");
127        }
128        if (evCheckEnabled) {
129                RealVector v = ed.getEigenvector(n-1);
130                double eps = 0.0;
131                for (int i = 0; i < 5; i++) {
132                        eps += Math.abs(v.getEntry(i));
133                        eps += Math.abs(v.getEntry(n-1-i));
134                }
135                if (eps > 1e-4) {
136                        throw new RuntimeException("lowest eigenfunction does not decay");
137                }
138        }
139    }
140    
141    /**
142     * Approximate the eigenfunction for the lowest eigenvalue of this operator, in the basis described above.
143     * 
144     * @param n number of basis vectors
145     * @param cutoff the "cutoff" \(R\)
146     * @return the lowest eigenvector in the given basis
147     */
148    public RealVector findLowestEigenvector(int n, double cutoff) {
149        RealMatrix a = this.asRealMatrix(n, cutoff);
150        EigenDecomposition e = new EigenDecomposition(a);
151        checkPlausibility(e, n);
152        return e.getEigenvector(a.getRowDimension()-1);
153    }
154
155    /**
156     * Approximate the spectrum of this operator, in the basis described above.
157     * The eigenvalues are returned in ascending order. 
158     *   
159     * @param n number of basis vectors
160     * @param cutoff the "cutoff" \(R\)
161     * @return the eigenvalues of the operator in this approximation
162     */
163    public double[] findSpectrum(int n, double cutoff) {
164        RealMatrix a = this.asRealMatrix(n, cutoff);
165        EigenDecomposition e = new EigenDecomposition(a);
166        checkPlausibility(e, n);
167        double[] evs =  e.getRealEigenvalues();
168        Arrays.sort(evs);
169        return evs;
170    }
171
172    /**
173     * Approximate the lowest eigenvalue of this operator, using the basis described above.
174     * 
175     * @param n number of basis vectors
176     * @param cutoff the "cutoff" \(R\)
177     * @return the lowest eigenvalue of the operator in this approximation
178     */
179    public double findLowestEigenvalue(int n, double cutoff) {
180        return findSpectrum(n, cutoff)[0];
181    }
182
183}