/**
* pmfg-kernel.cl - Planar Maximally Filtered Graph for Linguistics
*
* OpenCL kernels for GPU-accelerated PMFG construction.
* Integrates with atomese-simd via OpenclNode.
*
* Key insight: In linguistics, "planarity" means "no crossing dependencies"
* when words are laid out in sentence order. This is MUCH simpler than
* general graph planarity (no need for Boyer-Myrvold).
*
* Crossing test: edges (i,j) and (k,l) cross iff:
* (i < k < j < l) OR (k < i < l < j)
* where i,j,k,l are word positions in the sentence.
*/
#if defined(cl_khr_fp64)
# pragma OPENCL EXTENSION cl_khr_fp64: enable
#elif defined(cl_amd_fp64)
# pragma OPENCL EXTENSION cl_amd_fp64: enable
#else
# error double precision is not supported
#endif
// ═══════════════════════════════════════════════════════════════════
// HELPER: Check if two edges cross in linear word order
// ═══════════════════════════════════════════════════════════════════
inline bool edges_cross(int i, int j, int k, int l) {
// Normalize so i < j and k < l
if (i > j) { int t = i; i = j; j = t; }
if (k > l) { int t = k; k = l; l = t; }
// Check crossing: i < k < j < l OR k < i < l < j
return (i < k && k < j && j < l) || (k < i && i < l && l < j);
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 1: Compute all-pairs MI for a sentence
// ═══════════════════════════════════════════════════════════════════
/**
* compute_sentence_mi - Get MI values for all word pairs in sentences
*
* @param mi_matrix Global MI matrix [vocab_size × vocab_size]
* @param sentences Word indices [max_len × n_sentences]
* @param sentence_lens Length of each sentence
* @param pair_mi Output: MI for each pair [max_pairs × n_sentences]
* @param pair_indices Output: (pos1, pos2) for each pair
* @param pair_counts Output: number of pairs per sentence
* @param vocab_size Vocabulary size
* @param max_len Maximum sentence length
* @param sz Number of sentences
*/
__kernel void compute_sentence_mi(
__global const double* mi_matrix,
__global const int* sentences,
__global const int* sentence_lens,
__global double* pair_mi,
__global int* pair_indices,
__global int* pair_counts,
const int vocab_size,
const int max_len,
const unsigned long sz)
{
size_t sent_idx = get_global_id(0);
if (sent_idx >= sz) return;
int n = sentence_lens[sent_idx];
__global const int* words = sentences + sent_idx * max_len;
int max_pairs = max_len * (max_len - 1) / 2;
__global double* mi_out = pair_mi + sent_idx * max_pairs;
__global int* idx_out = pair_indices + sent_idx * max_pairs * 2;
int pair_count = 0;
for (int i = 0; i < n; i++) {
int w1 = words[i];
if (w1 < 0) continue;
for (int j = i + 1; j < n; j++) {
int w2 = words[j];
if (w2 < 0) continue;
// Store MI and positions
double mi = mi_matrix[w1 * vocab_size + w2];
mi_out[pair_count] = mi;
idx_out[pair_count * 2] = i;
idx_out[pair_count * 2 + 1] = j;
pair_count++;
}
}
pair_counts[sent_idx] = pair_count;
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 2: Sort pairs by MI (parallel bitonic sort)
// ═══════════════════════════════════════════════════════════════════
/**
* bitonic_sort_step - One step of bitonic sort
*
* Call this kernel log²(n) times with different j values.
*/
__kernel void bitonic_sort_step(
__global double* mi_values,
__global int* pair_indices,
const int j,
const int k,
const unsigned long sz)
{
size_t i = get_global_id(0);
if (i >= sz) return;
size_t ixj = i ^ j;
if (ixj > i) {
if ((i & k) == 0) {
// Sort ascending
if (mi_values[i] < mi_values[ixj]) {
// Swap MI
double tmp_mi = mi_values[i];
mi_values[i] = mi_values[ixj];
mi_values[ixj] = tmp_mi;
// Swap indices
int tmp_i1 = pair_indices[i * 2];
int tmp_i2 = pair_indices[i * 2 + 1];
pair_indices[i * 2] = pair_indices[ixj * 2];
pair_indices[i * 2 + 1] = pair_indices[ixj * 2 + 1];
pair_indices[ixj * 2] = tmp_i1;
pair_indices[ixj * 2 + 1] = tmp_i2;
}
} else {
// Sort descending
if (mi_values[i] > mi_values[ixj]) {
double tmp_mi = mi_values[i];
mi_values[i] = mi_values[ixj];
mi_values[ixj] = tmp_mi;
int tmp_i1 = pair_indices[i * 2];
int tmp_i2 = pair_indices[i * 2 + 1];
pair_indices[i * 2] = pair_indices[ixj * 2];
pair_indices[i * 2 + 1] = pair_indices[ixj * 2 + 1];
pair_indices[ixj * 2] = tmp_i1;
pair_indices[ixj * 2 + 1] = tmp_i2;
}
}
}
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 3: Greedy PMFG construction (main algorithm)
// ═══════════════════════════════════════════════════════════════════
/**
* pmfg_greedy_batch - Build PMFG for multiple sentences
*
* Greedy algorithm:
* 1. Sort pairs by MI (descending)
* 2. For each pair (highest MI first):
* - Check if adding it would cross existing edges
* - If not, add it
* - Stop when we have 3N-6 edges (maximally planar)
*
* @param sorted_mi MI values sorted descending [max_pairs × n_sent]
* @param sorted_indices Pair positions [max_pairs × 2 × n_sent]
* @param pair_counts Number of pairs per sentence
* @param sentence_lens Length of each sentence
* @param pmfg_edges Output: PMFG edges [max_edges × 2 × n_sent]
* @param pmfg_edge_counts Output: edges per sentence
* @param max_len Max sentence length
* @param sz Number of sentences
*/
__kernel void pmfg_greedy_batch(
__global const double* sorted_mi,
__global const int* sorted_indices,
__global const int* pair_counts,
__global const int* sentence_lens,
__global int* pmfg_edges,
__global int* pmfg_edge_counts,
const int max_len,
const unsigned long sz)
{
size_t sent_idx = get_global_id(0);
if (sent_idx >= sz) return;
int n = sentence_lens[sent_idx];
int n_pairs = pair_counts[sent_idx];
int max_pmfg_edges = 3 * n - 6; // Maximally planar
if (max_pmfg_edges < 0) max_pmfg_edges = 0;
int max_pairs = max_len * (max_len - 1) / 2;
__global const int* pairs = sorted_indices + sent_idx * max_pairs * 2;
__global int* edges = pmfg_edges + sent_idx * (3 * max_len - 6) * 2;
// Local storage for added edges
int added_edges[64]; // Max 30 edges × 2 = 60, round up
int n_edges = 0;
// Greedy: try to add each pair in MI order
for (int p = 0; p < n_pairs && n_edges < max_pmfg_edges; p++) {
int pos1 = pairs[p * 2];
int pos2 = pairs[p * 2 + 1];
// Check if this edge crosses any existing edge
bool crosses = false;
for (int e = 0; e < n_edges && !crosses; e++) {
int e1 = added_edges[e * 2];
int e2 = added_edges[e * 2 + 1];
crosses = edges_cross(pos1, pos2, e1, e2);
}
if (!crosses) {
// Add edge
added_edges[n_edges * 2] = pos1;
added_edges[n_edges * 2 + 1] = pos2;
edges[n_edges * 2] = pos1;
edges[n_edges * 2 + 1] = pos2;
n_edges++;
}
}
pmfg_edge_counts[sent_idx] = n_edges;
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 4: Extract Sections from PMFG
// ═══════════════════════════════════════════════════════════════════
/**
* extract_sections_batch - Convert PMFG edges to Sections
*
* Section = word + (left_connectors, right_connectors)
*
* @param sentences Word indices [max_len × n_sent]
* @param sentence_lens Sentence lengths
* @param pmfg_edges PMFG edges [max_edges × 2 × n_sent]
* @param pmfg_edge_counts Edges per sentence
* @param section_left Output: left connector count [max_len × n_sent]
* @param section_right Output: right connector count [max_len × n_sent]
* @param max_len Max sentence length
* @param sz Number of sentences
*/
__kernel void extract_sections_batch(
__global const int* sentences,
__global const int* sentence_lens,
__global const int* pmfg_edges,
__global const int* pmfg_edge_counts,
__global int* section_left,
__global int* section_right,
const int max_len,
const unsigned long sz)
{
size_t sent_idx = get_global_id(0);
if (sent_idx >= sz) return;
int n = sentence_lens[sent_idx];
int n_edges = pmfg_edge_counts[sent_idx];
__global const int* edges = pmfg_edges + sent_idx * (3 * max_len - 6) * 2;
__global int* left = section_left + sent_idx * max_len;
__global int* right = section_right + sent_idx * max_len;
// Initialize
for (int i = 0; i < n; i++) {
left[i] = 0;
right[i] = 0;
}
// Count connectors
for (int e = 0; e < n_edges; e++) {
int pos1 = edges[e * 2];
int pos2 = edges[e * 2 + 1];
if (pos1 < pos2) {
right[pos1]++; // pos1 connects right
left[pos2]++; // pos2 connects left
} else {
right[pos2]++;
left[pos1]++;
}
}
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 5: Aggregate Section counts across corpus
// ═══════════════════════════════════════════════════════════════════
/**
* aggregate_sections - Sum Section counts using atomic operations
*
* @param sentences Word indices
* @param sentence_lens Sentence lengths
* @param section_left Left connector counts per sentence
* @param section_right Right connector counts per sentence
* @param global_counts Output: global Section counts [vocab × max_conn × 2]
* @param vocab_size Vocabulary size
* @param max_len Max sentence length
* @param max_conn Max connectors per word
* @param sz Number of sentences
*/
__kernel void aggregate_sections(
__global const int* sentences,
__global const int* sentence_lens,
__global const int* section_left,
__global const int* section_right,
__global int* global_counts,
const int vocab_size,
const int max_len,
const int max_conn,
const unsigned long sz)
{
size_t sent_idx = get_global_id(0);
if (sent_idx >= sz) return;
int n = sentence_lens[sent_idx];
__global const int* words = sentences + sent_idx * max_len;
__global const int* left = section_left + sent_idx * max_len;
__global const int* right = section_right + sent_idx * max_len;
for (int i = 0; i < n; i++) {
int word = words[i];
if (word < 0 || word >= vocab_size) continue;
int left_count = left[i];
int right_count = right[i];
// Clamp to max_conn
if (left_count > max_conn) left_count = max_conn;
if (right_count > max_conn) right_count = max_conn;
// Atomic increment for this Section type
// Index: word * max_conn * 2 + left_count * 2 + 0 (left)
// word * max_conn * 2 + right_count * 2 + 1 (right)
int idx_l = word * max_conn * 2 + left_count * 2;
int idx_r = word * max_conn * 2 + right_count * 2 + 1;
atomic_add(&global_counts[idx_l], 1);
atomic_add(&global_counts[idx_r], 1);
}
}
// ═══════════════════════════════════════════════════════════════════
// KERNEL 6: Single-sentence PMFG (for interactive use)
// ═══════════════════════════════════════════════════════════════════
/**
* pmfg_single - Complete PMFG for one sentence
*
* All-in-one kernel for interactive parsing.
* Uses double arrays for atomese-simd FloatValue compatibility.
*
* @param mi_matrix Global MI matrix [vocab × vocab]
* @param words Word indices for this sentence (as doubles)
* @param n_words_arr Number of words (as double, single element)
* @param pmfg_edges Output: PMFG edges (as doubles)
* @param n_edges_out Output: number of edges (as double)
* @param vocab_size Vocabulary size
* @param sz Array length (for atomese-simd)
*/
__kernel void pmfg_single(
__global const double* mi_matrix,
__global const double* words,
__global const double* n_words_arr,
__global const double* vocab_size_arr,
__global double* pmfg_edges,
__global double* n_edges_out,
const unsigned long sz)
{
if (get_global_id(0) != 0) return;
int n = (int)n_words_arr[0];
int vocab_size = (int)vocab_size_arr[0];
if (n < 2) {
n_edges_out[0] = 0.0;
return;
}
int max_edges = 3 * n - 6;
if (max_edges < 0) max_edges = 0;
// Compute MI for all pairs and sort
double pair_mi[128]; // Max 66 pairs for 12 words
int pair_pos[256]; // pos1, pos2 for each pair
int n_pairs = 0;
for (int i = 0; i < n; i++) {
int w1 = (int)words[i];
for (int j = i + 1; j < n; j++) {
int w2 = (int)words[j];
pair_mi[n_pairs] = mi_matrix[w1 * vocab_size + w2];
pair_pos[n_pairs * 2] = i;
pair_pos[n_pairs * 2 + 1] = j;
n_pairs++;
}
}
// Simple insertion sort (small n, descending by MI)
for (int i = 1; i < n_pairs; i++) {
double key_mi = pair_mi[i];
int key_p1 = pair_pos[i * 2];
int key_p2 = pair_pos[i * 2 + 1];
int j = i - 1;
while (j >= 0 && pair_mi[j] < key_mi) {
pair_mi[j + 1] = pair_mi[j];
pair_pos[(j + 1) * 2] = pair_pos[j * 2];
pair_pos[(j + 1) * 2 + 1] = pair_pos[j * 2 + 1];
j--;
}
pair_mi[j + 1] = key_mi;
pair_pos[(j + 1) * 2] = key_p1;
pair_pos[(j + 1) * 2 + 1] = key_p2;
}
// Greedy PMFG construction
int added_edges[64];
int n_edges = 0;
for (int p = 0; p < n_pairs && n_edges < max_edges; p++) {
int pos1 = pair_pos[p * 2];
int pos2 = pair_pos[p * 2 + 1];
// Check crossings
bool crosses = false;
for (int e = 0; e < n_edges && !crosses; e++) {
crosses = edges_cross(pos1, pos2,
added_edges[e * 2],
added_edges[e * 2 + 1]);
}
if (!crosses) {
added_edges[n_edges * 2] = pos1;
added_edges[n_edges * 2 + 1] = pos2;
pmfg_edges[n_edges * 2] = (double)pos1;
pmfg_edges[n_edges * 2 + 1] = (double)pos2;
n_edges++;
}
}
n_edges_out[0] = (double)n_edges;
}
// Simple test kernel to verify OpenCL execution
__kernel void test_simple(
__global const double* input,
__global double* output,
const unsigned long sz)
{
size_t i = get_global_id(0);
if (i < sz) {
output[i] = input[i] * 2.0 + 42.0;
}
}