added pairwise grad function

Author: hexqnt

src/models/arctangent_step.rs

@@ -17,28 +17,45 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     amplitude * (slope * (x - x0)).atan() + offset
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 4);
+
+    let amplitude = param[0];
+    let slope = param[1];
+    let x0 = param[2];
+    let offset = param[3];
+    let z = slope * (x - x0);
+    let atan_z = z.atan();
+    let inv_den = 1.0 / (1.0 + z * z);
+
+    grad[0] = atan_z;
+    grad[1] = amplitude * (x - x0) * inv_den;
+    grad[2] = -amplitude * slope * inv_den;
+    grad[3] = 1.0;
+
+    amplitude * atan_z + offset
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let amplitude = param[0];
-    let slope = param[1];
-    let x0 = param[2];
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 4];
     let mut index = 0;
     while index < x_values.len() {
-        let x = x_values[index];
-        let z = slope * (x - x0);
-        let atan_z = z.atan();
-        let inv_den = 1.0 / (1.0 + z * z);
-        let residual = value_first[index];
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
 
-        gradient[0] += residual * atan_z;
-        gradient[1] += residual * (amplitude * (x - x0) * inv_den);
-        gradient[2] += residual * (-amplitude * slope * inv_den);
-        gradient[3] += residual;
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
+        gradient[3] += upstream * point_grad[3];
         index += 1;
     }
 }

src/models/arrhenius.rs

@@ -16,22 +16,37 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     prefactor * (temp_coeff / x).exp()
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 2);
+
+    let x = positive_x(x);
+    let prefactor = param[0];
+    let temp_coeff = param[1];
+    let exp_term = (temp_coeff / x).exp();
+
+    grad[0] = exp_term;
+    grad[1] = prefactor * exp_term / x;
+
+    prefactor * exp_term
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let prefactor = param[0];
-    let temp_coeff = param[1];
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 2];
     let mut index = 0;
     while index < x_values.len() {
-        let x = positive_x(x_values[index]);
-        let exp_term = (temp_coeff / x).exp();
-        let residual = value_first[index];
-        gradient[0] += residual * exp_term;
-        gradient[1] += residual * (prefactor * exp_term / x);
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
         index += 1;
     }
 }

src/models/bi_exponential.rs

@@ -17,29 +17,47 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     a1 * (-k1 * x).exp() + a2 * (-k2 * x).exp() + offset
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 5);
+
+    let a1 = param[0];
+    let k1 = param[1];
+    let a2 = param[2];
+    let k2 = param[3];
+    let offset = param[4];
+    let exp1 = (-k1 * x).exp();
+    let exp2 = (-k2 * x).exp();
+
+    grad[0] = exp1;
+    grad[1] = -a1 * x * exp1;
+    grad[2] = exp2;
+    grad[3] = -a2 * x * exp2;
+    grad[4] = 1.0;
+
+    a1 * exp1 + a2 * exp2 + offset
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let a1 = param[0];
-    let k1 = param[1];
-    let a2 = param[2];
-    let k2 = param[3];
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 5];
     let mut index = 0;
     while index < x_values.len() {
-        let x = x_values[index];
-        let exp1 = (-k1 * x).exp();
-        let exp2 = (-k2 * x).exp();
-        let residual = value_first[index];
-
-        gradient[0] += residual * exp1;
-        gradient[1] += residual * (-a1 * x * exp1);
-        gradient[2] += residual * exp2;
-        gradient[3] += residual * (-a2 * x * exp2);
-        gradient[4] += residual;
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
+        gradient[3] += upstream * point_grad[3];
+        gradient[4] += upstream * point_grad[4];
         index += 1;
     }
 }

src/models/damped_sinusoid.rs

@@ -18,31 +18,49 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     amplitude * (-damping * x).exp() * (omega * x + phi).sin() + offset
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 5);
+
+    let amplitude = param[0];
+    let damping = param[1];
+    let omega = param[2];
+    let phi = param[3];
+    let offset = param[4];
+    let exp_part = (-damping * x).exp();
+    let angle = omega * x + phi;
+    let sin_part = angle.sin();
+    let cos_part = angle.cos();
+
+    grad[0] = exp_part * sin_part;
+    grad[1] = -amplitude * x * exp_part * sin_part;
+    grad[2] = amplitude * exp_part * cos_part * x;
+    grad[3] = amplitude * exp_part * cos_part;
+    grad[4] = 1.0;
+
+    amplitude * exp_part * sin_part + offset
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let amplitude = param[0];
-    let damping = param[1];
-    let omega = param[2];
-    let phi = param[3];
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 5];
     let mut index = 0;
     while index < x_values.len() {
-        let x = x_values[index];
-        let exp_part = (-damping * x).exp();
-        let angle = omega * x + phi;
-        let sin_part = angle.sin();
-        let cos_part = angle.cos();
-        let residual = value_first[index];
-
-        gradient[0] += residual * exp_part * sin_part;
-        gradient[1] += residual * (-amplitude * x * exp_part * sin_part);
-        gradient[2] += residual * (amplitude * exp_part * cos_part * x);
-        gradient[3] += residual * (amplitude * exp_part * cos_part);
-        gradient[4] += residual;
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
+        gradient[3] += upstream * point_grad[3];
+        gradient[4] += upstream * point_grad[4];
         index += 1;
     }
 }

src/models/dispatch.rs

@@ -92,6 +92,7 @@ pub(crate) fn add_model_grad_unscaled(
     gradient: &mut [f64],
 ) {
     debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
     if family.is_polynomial() {
         polynomial::add_value_grad(x_values, param, value_first, gradient);

src/models/emg.rs

@@ -37,12 +37,35 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     }
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 5);
+    grad.fill(0.0);
+    value_at(param, x)
+}
+
 pub(super) fn add_value_grad(
-    _x_values: &[f64],
-    _param: &[f64],
-    _value_first: &[f64],
-    _gradient: &mut [f64],
+    x_values: &[f64],
+    param: &[f64],
+    value_first: &[f64],
+    gradient: &mut [f64],
 ) {
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
+
+    let mut point_grad = [0.0; 5];
+    let mut index = 0;
+    while index < x_values.len() {
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
+        gradient[3] += upstream * point_grad[3];
+        gradient[4] += upstream * point_grad[4];
+        index += 1;
+    }
 }
 
 pub(super) fn add_value_grad_raw_hessian(

src/models/exponential_basic.rs

@@ -15,23 +15,40 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     offset + amplitude * (-decay_rate * x).exp()
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 3);
+
+    let offset = param[0];
+    let amplitude = param[1];
+    let decay_rate = param[2];
+    let exp_part = (-decay_rate * x).exp();
+
+    grad[0] = 1.0;
+    grad[1] = exp_part;
+    grad[2] = -amplitude * x * exp_part;
+
+    offset + amplitude * exp_part
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let amplitude = param[1];
-    let decay_rate = param[2];
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 3];
     let mut index = 0;
     while index < x_values.len() {
-        let x = x_values[index];
-        let exp_part = (-decay_rate * x).exp();
-        let residual = value_first[index];
-        gradient[0] += residual;
-        gradient[1] += residual * exp_part;
-        gradient[2] += residual * (-amplitude * x * exp_part);
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
         index += 1;
     }
 }

src/models/exponential_half_life.rs

@@ -19,27 +19,43 @@ pub(super) fn value_at(param: &[f64], x: f64) -> f64 {
     offset + amplitude * exponent.exp()
 }
 
+#[inline]
+pub(super) fn value_grad_at(param: &[f64], x: f64, grad: &mut [f64]) -> f64 {
+    debug_assert_eq!(grad.len(), 3);
+
+    let offset = param[0];
+    let amplitude = param[1];
+    let half_life_raw = param[2];
+    let (half_life, d_c_raw) = positive_param_with_derivative(half_life_raw);
+    let exponent = -LN_2 * x / half_life;
+    let pow = exponent.exp();
+    let d_model_d_c = amplitude * pow * LN_2 * x / (half_life * half_life);
+
+    grad[0] = 1.0;
+    grad[1] = pow;
+    grad[2] = d_model_d_c * d_c_raw;
+
+    offset + amplitude * pow
+}
+
 pub(super) fn add_value_grad(
     x_values: &[f64],
     param: &[f64],
     value_first: &[f64],
     gradient: &mut [f64],
 ) {
-    let amplitude = param[1];
-    let half_life_raw = param[2];
-    let (half_life, d_c_raw) = positive_param_with_derivative(half_life_raw);
+    debug_assert_eq!(x_values.len(), value_first.len());
+    debug_assert_eq!(gradient.len(), param.len());
 
+    let mut point_grad = [0.0; 3];
     let mut index = 0;
     while index < x_values.len() {
-        let x = x_values[index];
-        let exponent = -LN_2 * x / half_life;
-        let pow = exponent.exp();
-        let residual = value_first[index];
-        let d_model_d_c = amplitude * pow * LN_2 * x / (half_life * half_life);
-
-        gradient[0] += residual;
-        gradient[1] += residual * pow;
-        gradient[2] += residual * d_model_d_c * d_c_raw;
+        let upstream = value_first[index];
+        value_grad_at(param, x_values[index], &mut point_grad);
+
+        gradient[0] += upstream * point_grad[0];
+        gradient[1] += upstream * point_grad[1];
+        gradient[2] += upstream * point_grad[2];
         index += 1;
     }
 }