isoakernel1.java

a molecular graph kernel based on iterative graph similarity and optimal assignments

    
    /** 
     * Evaluates the iterative optimal assignment kernel without normalization.
     * 
     * @param molA first argument molecuar graph.
     * @param molB second argument molecular graph.
     * @param alpha the weight of the recursive part of the update equation.
     * @param epsilon the desired precision of the result. The result will not deviate by more than epsilon from the true value, as measured by the maximum-norm based metric.
     * @param assignment  the indices of the optimal assignment. Pass null if not needed.
     * @param similarity  the computed similarity matrix. Pass null if not needed.
     * @return the value of the optimal assignment on the iterative similarity matrix.
     */
    
    public float evalRaw(MolecularGraph molA, MolecularGraph molB, float alpha, float epsilon, int[] assignment, Array2Float similarity) {
        // Set up equations.
        computeInternals(molA, molB);
        int iteration = 0;          // Current number of completed iterations.
        int numIterations = 1;      // Number of iterations that guarantee desired precision. Will be established later on.
        
        if (similarity != null && ( similarity.rows() != kk || similarity.cols() != nn )) throw new RuntimeException(String.format("Size of passed similarity matrix does not match size of passed molecules (expected %d x %d but got %d x %d).", kk, nn, similarity.rows(), similarity.cols()));

        // Check special cases.
        if (kknn == 1) {    // Comparison of two atoms. 
            if (assignment != null) assignment[0] = 0;
            final float result = (kvec[0] + 1)/2; 
            if (similarity != null) similarity.set(0, 0, result);
            return result; 
        }  
        
        // Initialize x to vertex kernel.
        System.arraycopy(kvec, 0, x, 0, kknn);
        
        // Save first value of x for computation of termination criterion.
        System.arraycopy(x, 0, xzero, 0, kknn);
        
        // Iteratively compute new values for x.
        do {
            for (int index = kknn; --index >= 0; ) {
                // Compute the contribution due to neighbour similarity.
                float contrib = 0.f;  
                for (int term = pIndices[index][pIndicesRow-1]; --term >= 0; ) {  // Iterate over max terms.
                    float termContrib = 0.f;
                    for (int summand = pIndices[index][(term+1)*tauOne-1]; --summand >= 0; ) 
                        termContrib += pFactors[index][term*tau+summand] * x[pIndices[index][term*tauOne+summand]];                 
                    contrib = Math.max(contrib, termContrib);
                }
                
                // Update the component.
                x[index] = (1.f-alpha)*kvec[index] + alpha*contrib; // Was: x[index] = (kvec[index] + contrib)/2.f;
            }
            
            // Check termination condition.
            if (++iteration == 1) {  // First iteration, establish number of iterations to do.
                final float distance = maxMetric(x, xzero);
                if (distance <= epsilon*(1.f-alpha)/alpha) numIterations = 1;  // d <= e(1-a)/a <--> e(1-a)/d >= a --> log_{a} (e(1-a)/d) <= 1 --> k >= 1.   
                else numIterations = (int) Math.ceil( Math.log( epsilon*(1.f-alpha)/distance ) / Math.log(alpha) );
            }
        } while (iteration < numIterations);
        
        // Compute optimal assignment.
        for (int row = kk; --row >= 0; ) for(int col = nn; --col >= 0; ) assignTemp.set(row, col, x[col*kk+row]);
        final float result = hungAlg.hungarianAlgorithm(assignTemp, OptimizationType.MAX, assignment);
        
        if (similarity != null) for (int row = kk; --row >= 0; ) for (int col = nn; --col >= 0; ) similarity.set(row, col, x[col*kk+row]);
        
        return result;
    }
    
    /** 
     * Evalutes the iterative optimal assignment kernel with normalization.
     * 
     * The result will be in the range [0,1].
     */
    public float evalNorm(MolecularGraph molA, MolecularGraph molB, float alpha, float epsilon, int[] assignment, Array2Float similarity) {
        float result = evalRaw(molA, molB, alpha, epsilon, assignment, similarity);
        //result /= Math.sqrt(evalRaw(molA, molA, epsilon, null, null) * evalRaw(molB, molB, epsilon, null, null));
        result /= Math.sqrt(molA.numAtoms()*molB.numAtoms());   // Make use of k(G,G) = |V|.
        return result;
    }
    
    /** The maximum-norm based metric is used to compute the number of iterations. */
    private static float maxMetric(float[] x, float[] y) {
        float result = 0.f;
        for (int i = x.length; --i >= 0; ) result = Math.max(result, Math.abs(x[i]-y[i]));
        return result;
    }

    /** Allocates memory and precomputes neighbour assignment matrices. */ 
    private void computeInternals(MolecularGraph molA, MolecularGraph molB) {
        // Shortcuts.
        kk = molA.numAtoms();       // Number of atoms in first molecule.   
        nn = molB.numAtoms();       // Number of atoms in second molecule.
        kknn = kk*nn;               // Number of linearized indices.

        // Allocate memory if necessary.
        if (kvec.length < kknn) kvec = new float[kknn];
        if (pFactors.length < kknn) pFactors = new float[kknn][pFactorsRow];
        if (pIndices.length < kknn) pIndices = new int[kknn][pIndicesRow];
        if (x.length < kknn) x = new float[kknn];           // Only allocated, not initialized.
        if (xzero.length < kknn) xzero = new float[kknn];   // Only allocated, not initialized.
        assignTemp.redim(kk, nn);   // Only allocated, not initialized. Redim is necessary because the Hungarian algorithm expects the right matrix dimensions.
        
        // Compute vertex matchings and P matrices.
        for (int index = 0; index < kknn; ++index) {  // Iterate over linearized entries of similarity vector/matrix.
            // Determine row and column for current index.
            int row = index % kk;
            int col = index / kk;

            // Vertex kernel.
            kvec[index] = vertexKernel.eval(molA, row, molB, col);
            
            // P matrices.
            //
            // To compute the assignment matrices, for each vertex pair, all possible neighbour combinations 
            // are enumerated and stored as terms (at most tau!, each one a neighbourhood assignment) with at most tau summands each.
            // The indices give the involved vertex pairs and the factors include division by degree and the edge kernel.
            final int numNeighbA = molA.numNeighbours(row);
            final int numNeighbB = molB.numNeighbours(col);
            
            if (numNeighbA < numNeighbB) {  // Vertex row has fewer neighbours than vertex col.
                final int[][] permPref = Permutations.prefixesPrecomputed(numNeighbA, numNeighbB);  // Possible assignments.
                pIndices[index][pIndicesRow-1] = permPref.length;                                   // Store number of terms in last entry.
                for (int term = permPref.length; --term >= 0; ) {                                   // Determine each term.
                    pIndices[index][(term+1)*tauOne-1] = numNeighbA;                                // Store number of summands.
                    for (int summand = numNeighbA; --summand >= 0; ) {                              // Determine each summand.
                        final int neighbA = molA.neighbour(row, summand);
                        final int neighbB = molB.neighbour(col, permPref[term][summand]);
                        final int ea = molA.bondIndex(row, neighbA);
                        final int eb = molB.bondIndex(col, neighbB);
                        pFactors[index][term*tau + summand] = (1.f/numNeighbB) * edgeKernel.eval(molA, ea, molB, eb);
                        pIndices[index][term*tauOne + summand] = neighbB*kk + neighbA;
                    }
                }
            } else {  // Vertex col has fewer or equally many neighbours than vertex row. This case is symmetric to the first one.
                final int[][] permPref = Permutations.prefixesPrecomputed(numNeighbB, numNeighbA);
                pIndices[index][pIndicesRow-1] = permPref.length;
                for (int term = permPref.length; --term >= 0; ) {
                    pIndices[index][(term+1)*tauOne-1] = numNeighbB;
                    for (int summand = numNeighbB; --summand >= 0; ) {
                        final int neighbA = molA.neighbour(row, permPref[term][summand]);
                        final int neighbB = molB.neighbour(col, summand);
                        final int ea = molA.bondIndex(row, neighbA);
                        final int eb = molB.bondIndex(col, neighbB);
                        pFactors[index][term*tau + summand] = (1.f/numNeighbA) * edgeKernel.eval(molA, ea, molB, eb);
                        pIndices[index][term*tauOne + summand] = neighbB*kk + neighbA;
                    }                   
                }
            }
        }
    }
}