1use std::fmt::Debug;
21use std::mem::size_of_val;
22use std::sync::Arc;
23
24use arrow::array::{
25 Array, AsArray, BooleanArray, Float64Array, NullBufferBuilder, UInt64Array,
26 downcast_array,
27};
28use arrow::compute::{and, filter, is_not_null};
29use arrow::datatypes::{FieldRef, Float64Type, UInt64Type};
30use arrow::{
31 array::ArrayRef,
32 datatypes::{DataType, Field},
33};
34use datafusion_expr::{EmitTo, GroupsAccumulator};
35use datafusion_functions_aggregate_common::aggregate::groups_accumulator::accumulate::accumulate_multiple;
36use log::debug;
37
38use crate::covariance::CovarianceAccumulator;
39use crate::stddev::StddevAccumulator;
40use datafusion_common::{Result, ScalarValue};
41use datafusion_expr::{
42 Accumulator, AggregateUDFImpl, Documentation, Signature, Volatility,
43 function::{AccumulatorArgs, StateFieldsArgs},
44 utils::format_state_name,
45};
46use datafusion_functions_aggregate_common::stats::StatsType;
47use datafusion_macros::user_doc;
48
49make_udaf_expr_and_func!(
50 Correlation,
51 corr,
52 y x,
53 "Correlation between two numeric values.",
54 corr_udaf
55);
56
57#[user_doc(
58 doc_section(label = "Statistical Functions"),
59 description = "Returns the coefficient of correlation between two numeric values.",
60 syntax_example = "corr(expression1, expression2)",
61 sql_example = r#"```sql
62> SELECT corr(column1, column2) FROM table_name;
63+--------------------------------+
64| corr(column1, column2) |
65+--------------------------------+
66| 0.85 |
67+--------------------------------+
68```"#,
69 standard_argument(name = "expression1", prefix = "First"),
70 standard_argument(name = "expression2", prefix = "Second")
71)]
72#[derive(Debug, PartialEq, Eq, Hash)]
73pub struct Correlation {
74 signature: Signature,
75}
76
77impl Default for Correlation {
78 fn default() -> Self {
79 Self::new()
80 }
81}
82
83impl Correlation {
84 pub fn new() -> Self {
86 Self {
87 signature: Signature::exact(
88 vec![DataType::Float64, DataType::Float64],
89 Volatility::Immutable,
90 )
91 .with_parameter_names(vec!["y".to_string(), "x".to_string()])
92 .expect("valid parameter names for corr"),
93 }
94 }
95}
96
97impl AggregateUDFImpl for Correlation {
98 fn name(&self) -> &str {
99 "corr"
100 }
101
102 fn signature(&self) -> &Signature {
103 &self.signature
104 }
105
106 fn return_type(&self, _arg_types: &[DataType]) -> Result<DataType> {
107 Ok(DataType::Float64)
108 }
109
110 fn accumulator(&self, _acc_args: AccumulatorArgs) -> Result<Box<dyn Accumulator>> {
111 Ok(Box::new(CorrelationAccumulator::try_new()?))
112 }
113
114 fn state_fields(&self, args: StateFieldsArgs) -> Result<Vec<FieldRef>> {
115 let name = args.name;
116 Ok(vec![
117 Field::new(format_state_name(name, "count"), DataType::UInt64, true),
118 Field::new(format_state_name(name, "mean1"), DataType::Float64, true),
119 Field::new(format_state_name(name, "m2_1"), DataType::Float64, true),
120 Field::new(format_state_name(name, "mean2"), DataType::Float64, true),
121 Field::new(format_state_name(name, "m2_2"), DataType::Float64, true),
122 Field::new(
123 format_state_name(name, "algo_const"),
124 DataType::Float64,
125 true,
126 ),
127 ]
128 .into_iter()
129 .map(Arc::new)
130 .collect())
131 }
132
133 fn documentation(&self) -> Option<&Documentation> {
134 self.doc()
135 }
136
137 fn groups_accumulator_supported(&self, _args: AccumulatorArgs) -> bool {
138 true
139 }
140
141 fn create_groups_accumulator(
142 &self,
143 _args: AccumulatorArgs,
144 ) -> Result<Box<dyn GroupsAccumulator>> {
145 debug!("GroupsAccumulator is created for aggregate function `corr(c1, c2)`");
146 Ok(Box::new(CorrelationGroupsAccumulator::new()))
147 }
148}
149
150#[derive(Debug)]
152pub struct CorrelationAccumulator {
153 covar: CovarianceAccumulator,
154 stddev1: StddevAccumulator,
155 stddev2: StddevAccumulator,
156}
157
158impl CorrelationAccumulator {
159 pub fn try_new() -> Result<Self> {
161 Ok(Self {
162 covar: CovarianceAccumulator::try_new(StatsType::Population)?,
163 stddev1: StddevAccumulator::try_new(StatsType::Population)?,
164 stddev2: StddevAccumulator::try_new(StatsType::Population)?,
165 })
166 }
167}
168
169impl Accumulator for CorrelationAccumulator {
170 fn update_batch(&mut self, values: &[ArrayRef]) -> Result<()> {
171 let values = if values[0].null_count() != 0 || values[1].null_count() != 0 {
177 let mask = and(&is_not_null(&values[0])?, &is_not_null(&values[1])?)?;
178 let values1 = filter(&values[0], &mask)?;
179 let values2 = filter(&values[1], &mask)?;
180
181 vec![values1, values2]
182 } else {
183 values.to_vec()
184 };
185
186 self.covar.update_batch(&values)?;
187 self.stddev1.update_batch(&values[0..1])?;
188 self.stddev2.update_batch(&values[1..2])?;
189 Ok(())
190 }
191
192 fn evaluate(&mut self) -> Result<ScalarValue> {
193 let covar = self.covar.evaluate()?;
194 let stddev1 = self.stddev1.evaluate()?;
195 let stddev2 = self.stddev2.evaluate()?;
196
197 let mean1 = self.covar.get_mean1();
200 let mean2 = self.covar.get_mean2();
201
202 if mean1.is_nan() && mean2.is_nan() {
204 return Ok(ScalarValue::Float64(Some(f64::NAN)));
205 }
206 let n = self.covar.get_count();
207 if mean1.is_nan() || mean2.is_nan() || n < 2 {
208 return Ok(ScalarValue::Float64(None));
209 }
210
211 if let ScalarValue::Float64(Some(c)) = covar
212 && let ScalarValue::Float64(Some(s1)) = stddev1
213 && let ScalarValue::Float64(Some(s2)) = stddev2
214 {
215 if s1 == 0_f64 || s2 == 0_f64 {
216 return Ok(ScalarValue::Float64(None));
217 } else {
218 return Ok(ScalarValue::Float64(Some(c / s1 / s2)));
219 }
220 }
221
222 Ok(ScalarValue::Float64(None))
223 }
224
225 fn size(&self) -> usize {
226 size_of_val(self) - size_of_val(&self.covar) + self.covar.size()
227 - size_of_val(&self.stddev1)
228 + self.stddev1.size()
229 - size_of_val(&self.stddev2)
230 + self.stddev2.size()
231 }
232
233 fn state(&mut self) -> Result<Vec<ScalarValue>> {
234 Ok(vec![
235 ScalarValue::from(self.covar.get_count()),
236 ScalarValue::from(self.covar.get_mean1()),
237 ScalarValue::from(self.stddev1.get_m2()),
238 ScalarValue::from(self.covar.get_mean2()),
239 ScalarValue::from(self.stddev2.get_m2()),
240 ScalarValue::from(self.covar.get_algo_const()),
241 ])
242 }
243
244 fn merge_batch(&mut self, states: &[ArrayRef]) -> Result<()> {
245 let states_c = [
246 Arc::clone(&states[0]),
247 Arc::clone(&states[1]),
248 Arc::clone(&states[3]),
249 Arc::clone(&states[5]),
250 ];
251 let states_s1 = [
252 Arc::clone(&states[0]),
253 Arc::clone(&states[1]),
254 Arc::clone(&states[2]),
255 ];
256 let states_s2 = [
257 Arc::clone(&states[0]),
258 Arc::clone(&states[3]),
259 Arc::clone(&states[4]),
260 ];
261
262 self.covar.merge_batch(&states_c)?;
263 self.stddev1.merge_batch(&states_s1)?;
264 self.stddev2.merge_batch(&states_s2)?;
265 Ok(())
266 }
267
268 fn retract_batch(&mut self, values: &[ArrayRef]) -> Result<()> {
269 let values = if values[0].null_count() != 0 || values[1].null_count() != 0 {
270 let mask = and(&is_not_null(&values[0])?, &is_not_null(&values[1])?)?;
271 let values1 = filter(&values[0], &mask)?;
272 let values2 = filter(&values[1], &mask)?;
273
274 vec![values1, values2]
275 } else {
276 values.to_vec()
277 };
278
279 self.covar.retract_batch(&values)?;
280 self.stddev1.retract_batch(&values[0..1])?;
281 self.stddev2.retract_batch(&values[1..2])?;
282 Ok(())
283 }
284
285 fn supports_retract_batch(&self) -> bool {
286 true
287 }
288}
289
290#[derive(Default)]
291pub struct CorrelationGroupsAccumulator {
292 count: Vec<u64>,
296 sum_x: Vec<f64>,
298 sum_y: Vec<f64>,
300 sum_xy: Vec<f64>,
302 sum_xx: Vec<f64>,
304 sum_yy: Vec<f64>,
306}
307
308impl CorrelationGroupsAccumulator {
309 pub fn new() -> Self {
310 Default::default()
311 }
312}
313
314fn accumulate_correlation_states(
320 group_indices: &[usize],
321 state_arrays: (
322 &UInt64Array, &Float64Array, &Float64Array, &Float64Array, &Float64Array, &Float64Array, ),
329 mut value_fn: impl FnMut(usize, u64, &[f64]),
330) {
331 let (counts, sum_x, sum_y, sum_xy, sum_xx, sum_yy) = state_arrays;
332
333 assert_eq!(counts.null_count(), 0);
334 assert_eq!(sum_x.null_count(), 0);
335 assert_eq!(sum_y.null_count(), 0);
336 assert_eq!(sum_xy.null_count(), 0);
337 assert_eq!(sum_xx.null_count(), 0);
338 assert_eq!(sum_yy.null_count(), 0);
339
340 let counts_values = counts.values().as_ref();
341 let sum_x_values = sum_x.values().as_ref();
342 let sum_y_values = sum_y.values().as_ref();
343 let sum_xy_values = sum_xy.values().as_ref();
344 let sum_xx_values = sum_xx.values().as_ref();
345 let sum_yy_values = sum_yy.values().as_ref();
346
347 for (idx, &group_idx) in group_indices.iter().enumerate() {
348 let row = [
349 sum_x_values[idx],
350 sum_y_values[idx],
351 sum_xy_values[idx],
352 sum_xx_values[idx],
353 sum_yy_values[idx],
354 ];
355 value_fn(group_idx, counts_values[idx], &row);
356 }
357}
358
359impl GroupsAccumulator for CorrelationGroupsAccumulator {
375 fn update_batch(
376 &mut self,
377 values: &[ArrayRef],
378 group_indices: &[usize],
379 opt_filter: Option<&BooleanArray>,
380 total_num_groups: usize,
381 ) -> Result<()> {
382 self.count.resize(total_num_groups, 0);
383 self.sum_x.resize(total_num_groups, 0.0);
384 self.sum_y.resize(total_num_groups, 0.0);
385 self.sum_xy.resize(total_num_groups, 0.0);
386 self.sum_xx.resize(total_num_groups, 0.0);
387 self.sum_yy.resize(total_num_groups, 0.0);
388
389 let array_x = downcast_array::<Float64Array>(&values[0]);
390 let array_y = downcast_array::<Float64Array>(&values[1]);
391
392 accumulate_multiple(
393 group_indices,
394 &[&array_x, &array_y],
395 opt_filter,
396 |group_index, batch_index, columns| {
397 let x = columns[0].value(batch_index);
398 let y = columns[1].value(batch_index);
399 self.count[group_index] += 1;
400 self.sum_x[group_index] += x;
401 self.sum_y[group_index] += y;
402 self.sum_xy[group_index] += x * y;
403 self.sum_xx[group_index] += x * x;
404 self.sum_yy[group_index] += y * y;
405 },
406 );
407
408 Ok(())
409 }
410
411 fn evaluate(&mut self, emit_to: EmitTo) -> Result<ArrayRef> {
412 let counts = emit_to.take_needed(&mut self.count);
414 let sum_xs = emit_to.take_needed(&mut self.sum_x);
415 let sum_ys = emit_to.take_needed(&mut self.sum_y);
416 let sum_xys = emit_to.take_needed(&mut self.sum_xy);
417 let sum_xxs = emit_to.take_needed(&mut self.sum_xx);
418 let sum_yys = emit_to.take_needed(&mut self.sum_yy);
419
420 let n = counts.len();
421 let mut values = Vec::with_capacity(n);
422 let mut nulls = NullBufferBuilder::new(n);
423
424 for i in 0..n {
433 let count = counts[i];
434 let sum_x = sum_xs[i];
435 let sum_y = sum_ys[i];
436 let sum_xy = sum_xys[i];
437 let sum_xx = sum_xxs[i];
438 let sum_yy = sum_yys[i];
439
440 if sum_x.is_nan() && sum_y.is_nan() {
443 values.push(f64::NAN);
445 nulls.append_non_null();
446 continue;
447 } else if count < 2 || sum_x.is_nan() || sum_y.is_nan() {
448 values.push(0.0);
450 nulls.append_null();
451 continue;
452 }
453
454 let mean_x = sum_x / count as f64;
455 let mean_y = sum_y / count as f64;
456
457 let numerator = sum_xy - sum_x * mean_y;
458 let denominator =
459 ((sum_xx - sum_x * mean_x) * (sum_yy - sum_y * mean_y)).sqrt();
460
461 if denominator == 0.0 {
462 values.push(0.0);
463 nulls.append_null();
464 } else {
465 values.push(numerator / denominator);
466 nulls.append_non_null();
467 }
468 }
469
470 Ok(Arc::new(Float64Array::new(values.into(), nulls.finish())))
471 }
472
473 fn state(&mut self, emit_to: EmitTo) -> Result<Vec<ArrayRef>> {
474 let count = emit_to.take_needed(&mut self.count);
476 let sum_x = emit_to.take_needed(&mut self.sum_x);
477 let sum_y = emit_to.take_needed(&mut self.sum_y);
478 let sum_xy = emit_to.take_needed(&mut self.sum_xy);
479 let sum_xx = emit_to.take_needed(&mut self.sum_xx);
480 let sum_yy = emit_to.take_needed(&mut self.sum_yy);
481
482 Ok(vec![
483 Arc::new(UInt64Array::from(count)),
484 Arc::new(Float64Array::from(sum_x)),
485 Arc::new(Float64Array::from(sum_y)),
486 Arc::new(Float64Array::from(sum_xy)),
487 Arc::new(Float64Array::from(sum_xx)),
488 Arc::new(Float64Array::from(sum_yy)),
489 ])
490 }
491
492 fn convert_to_state(
493 &self,
494 values: &[ArrayRef],
495 opt_filter: Option<&BooleanArray>,
496 ) -> Result<Vec<ArrayRef>> {
497 assert_eq!(values.len(), 2, "two arguments to convert_to_state");
498 let array_x = downcast_array::<Float64Array>(&values[0]);
499 let array_y = downcast_array::<Float64Array>(&values[1]);
500
501 let len = array_x.len();
502 let mut counts = Vec::with_capacity(len);
503 let mut sum_x = Vec::with_capacity(len);
504 let mut sum_y = Vec::with_capacity(len);
505 let mut sum_xy = Vec::with_capacity(len);
506 let mut sum_xx = Vec::with_capacity(len);
507 let mut sum_yy = Vec::with_capacity(len);
508
509 for row in 0..len {
510 let included = array_x.is_valid(row)
511 && array_y.is_valid(row)
512 && opt_filter
513 .is_none_or(|filter| filter.is_valid(row) && filter.value(row));
514 if included {
515 let x = array_x.value(row);
516 let y = array_y.value(row);
517 counts.push(1);
518 sum_x.push(x);
519 sum_y.push(y);
520 sum_xy.push(x * y);
521 sum_xx.push(x * x);
522 sum_yy.push(y * y);
523 } else {
524 counts.push(0);
525 sum_x.push(0.0);
526 sum_y.push(0.0);
527 sum_xy.push(0.0);
528 sum_xx.push(0.0);
529 sum_yy.push(0.0);
530 }
531 }
532
533 Ok(vec![
534 Arc::new(UInt64Array::from(counts)),
535 Arc::new(Float64Array::from(sum_x)),
536 Arc::new(Float64Array::from(sum_y)),
537 Arc::new(Float64Array::from(sum_xy)),
538 Arc::new(Float64Array::from(sum_xx)),
539 Arc::new(Float64Array::from(sum_yy)),
540 ])
541 }
542 fn merge_batch(
543 &mut self,
544 values: &[ArrayRef],
545 group_indices: &[usize],
546 total_num_groups: usize,
547 ) -> Result<()> {
548 self.count.resize(total_num_groups, 0);
550 self.sum_x.resize(total_num_groups, 0.0);
551 self.sum_y.resize(total_num_groups, 0.0);
552 self.sum_xy.resize(total_num_groups, 0.0);
553 self.sum_xx.resize(total_num_groups, 0.0);
554 self.sum_yy.resize(total_num_groups, 0.0);
555
556 let partial_counts = values[0].as_primitive::<UInt64Type>();
558 let partial_sum_x = values[1].as_primitive::<Float64Type>();
559 let partial_sum_y = values[2].as_primitive::<Float64Type>();
560 let partial_sum_xy = values[3].as_primitive::<Float64Type>();
561 let partial_sum_xx = values[4].as_primitive::<Float64Type>();
562 let partial_sum_yy = values[5].as_primitive::<Float64Type>();
563
564 accumulate_correlation_states(
565 group_indices,
566 (
567 partial_counts,
568 partial_sum_x,
569 partial_sum_y,
570 partial_sum_xy,
571 partial_sum_xx,
572 partial_sum_yy,
573 ),
574 |group_index, count, values| {
575 self.count[group_index] += count;
576 self.sum_x[group_index] += values[0];
577 self.sum_y[group_index] += values[1];
578 self.sum_xy[group_index] += values[2];
579 self.sum_xx[group_index] += values[3];
580 self.sum_yy[group_index] += values[4];
581 },
582 );
583
584 Ok(())
585 }
586
587 fn size(&self) -> usize {
588 self.count.capacity() * size_of::<u64>()
589 + self.sum_x.capacity() * size_of::<f64>()
590 + self.sum_y.capacity() * size_of::<f64>()
591 + self.sum_xy.capacity() * size_of::<f64>()
592 + self.sum_xx.capacity() * size_of::<f64>()
593 + self.sum_yy.capacity() * size_of::<f64>()
594 }
595}
596
597#[cfg(test)]
598mod tests {
599 use super::*;
600
601 #[test]
602 fn test_accumulate_correlation_states() {
603 let group_indices = vec![0, 1, 0, 1];
605 let counts = UInt64Array::from(vec![1, 2, 3, 4]);
606 let sum_x = Float64Array::from(vec![10.0, 20.0, 30.0, 40.0]);
607 let sum_y = Float64Array::from(vec![1.0, 2.0, 3.0, 4.0]);
608 let sum_xy = Float64Array::from(vec![10.0, 40.0, 90.0, 160.0]);
609 let sum_xx = Float64Array::from(vec![100.0, 400.0, 900.0, 1600.0]);
610 let sum_yy = Float64Array::from(vec![1.0, 4.0, 9.0, 16.0]);
611
612 let mut accumulated = vec![];
613 accumulate_correlation_states(
614 &group_indices,
615 (&counts, &sum_x, &sum_y, &sum_xy, &sum_xx, &sum_yy),
616 |group_idx, count, values| {
617 accumulated.push((group_idx, count, values.to_vec()));
618 },
619 );
620
621 let expected = vec![
622 (0, 1, vec![10.0, 1.0, 10.0, 100.0, 1.0]),
623 (1, 2, vec![20.0, 2.0, 40.0, 400.0, 4.0]),
624 (0, 3, vec![30.0, 3.0, 90.0, 900.0, 9.0]),
625 (1, 4, vec![40.0, 4.0, 160.0, 1600.0, 16.0]),
626 ];
627 assert_eq!(accumulated, expected);
628
629 let counts = UInt64Array::from(vec![Some(1), None, Some(3), Some(4)]);
631 let sum_x = Float64Array::from(vec![10.0, 20.0, 30.0, 40.0]);
632 let sum_y = Float64Array::from(vec![1.0, 2.0, 3.0, 4.0]);
633 let sum_xy = Float64Array::from(vec![10.0, 40.0, 90.0, 160.0]);
634 let sum_xx = Float64Array::from(vec![100.0, 400.0, 900.0, 1600.0]);
635 let sum_yy = Float64Array::from(vec![1.0, 4.0, 9.0, 16.0]);
636
637 let result = std::panic::catch_unwind(|| {
638 accumulate_correlation_states(
639 &group_indices,
640 (&counts, &sum_x, &sum_y, &sum_xy, &sum_xx, &sum_yy),
641 |_, _, _| {},
642 )
643 });
644 assert!(result.is_err());
645 }
646
647 #[test]
648 fn convert_to_state_roundtrips_through_merge() -> Result<()> {
649 let x = Arc::new(Float64Array::from(vec![
650 Some(1.0),
651 Some(2.0),
652 None,
653 Some(4.0),
654 Some(8.0),
655 Some(16.0),
656 Some(32.0),
657 ])) as ArrayRef;
658 let y = Arc::new(Float64Array::from(vec![
659 Some(2.0),
660 Some(4.0),
661 Some(6.0),
662 None,
663 Some(16.0),
664 Some(32.0),
665 Some(64.0),
666 ])) as ArrayRef;
667 let filter = BooleanArray::from(vec![
668 Some(true),
669 Some(false),
670 Some(true),
671 Some(true),
672 None,
673 Some(true),
674 Some(true),
675 ]);
676 let values = vec![x, y];
677 let group_indices = vec![0, 1, 0, 1, 0, 0, 0];
678
679 let mut direct = CorrelationGroupsAccumulator::new();
680 direct.update_batch(&values, &group_indices, Some(&filter), 2)?;
681 let direct = direct.evaluate(EmitTo::All)?;
682
683 let converter = CorrelationGroupsAccumulator::new();
684 let state = converter.convert_to_state(&values, Some(&filter))?;
685 let mut merged = CorrelationGroupsAccumulator::new();
686 merged.merge_batch(&state, &group_indices, 2)?;
687 let merged = merged.evaluate(EmitTo::All)?;
688
689 assert_eq!(
690 direct.as_any().downcast_ref::<Float64Array>().unwrap(),
691 merged.as_any().downcast_ref::<Float64Array>().unwrap()
692 );
693 Ok(())
694 }
695
696 #[test]
697 fn convert_to_state_preserves_empty_and_filtered_rows() -> Result<()> {
698 let converter = CorrelationGroupsAccumulator::new();
699 let empty_values = vec![
700 Arc::new(Float64Array::from(Vec::<Option<f64>>::new())) as ArrayRef,
701 Arc::new(Float64Array::from(Vec::<Option<f64>>::new())) as ArrayRef,
702 ];
703 let state = converter.convert_to_state(&empty_values, None)?;
704 for state_array in &state {
705 assert_eq!(state_array.len(), 0);
706 assert_eq!(state_array.null_count(), 0);
707 }
708
709 let values = vec![
710 Arc::new(Float64Array::from(vec![Some(1.0), Some(2.0), None])) as ArrayRef,
711 Arc::new(Float64Array::from(vec![Some(2.0), None, Some(4.0)])) as ArrayRef,
712 ];
713 let filter = BooleanArray::from(vec![Some(false), None, Some(false)]);
714 let group_indices = vec![0, 1, 0];
715 let state = converter.convert_to_state(&values, Some(&filter))?;
716 for state_array in &state {
717 assert_eq!(state_array.len(), values[0].len());
718 assert_eq!(state_array.null_count(), 0);
719 }
720
721 let counts = state[0].as_any().downcast_ref::<UInt64Array>().unwrap();
722 assert_eq!(counts, &UInt64Array::from(vec![0, 0, 0]));
723
724 let mut merged = CorrelationGroupsAccumulator::new();
725 merged.merge_batch(&state, &group_indices, 2)?;
726 let result = merged.evaluate(EmitTo::All)?;
727 let result = result.as_any().downcast_ref::<Float64Array>().unwrap();
728 assert_eq!(result.len(), 2);
729 assert_eq!(result.null_count(), 2);
730 Ok(())
731 }
732}