# Grid Refinement Algorithm Comparison: Fortran vs Rust
## Overview
Comparing the TWOPNT grid refinement algorithm from the original Fortran implementation with the Rust rewrite attempt.
## Key Differences Found
### 1. **Algorithm Structure**
**Fortran Original (`refine` subroutine)**:
- **Location**: Lines ~3800-4200 in `twopnt_core.f90`
- **Structure**: Single comprehensive subroutine with multiple phases
- **Error Handling**: Extensive validation with detailed error messages
**Rust Rewrite**:
- **Location**: `adaptive_grid_twopoint.rs`
- **Structure**: Object-oriented with `Grid` struct and methods
- **Error Handling**: Basic assertions and logging
### 2. **Mathematical Criteria**
**Fortran**:
```fortran
! Component significance test
if (.not. abs(range)>toler0*max(one,maxmag)) cycle active_components
! Gradient threshold (toler1 = 0.2)
if (toler1 * range<differ) vary1(k) = vary1(k) + 1
! Curvature threshold (toler2 = 0.2)
if (toler2 * range < differ) vary2(k) = vary2(k) + 1
```
**Rust**:
```rust
// Component significance test
if !(range_y_j.abs() > abs_tol * f64::max(1.0, maxmag_j)) {
continue;
}
// Gradient threshold (d = 0.2 equivalent)
if dy_i.abs() > delta {
vary_y[i] += 1;
}
// Curvature threshold (g = 0.2 equivalent)
if differ_dy > gamma {
vary_dy_dx[i] += 1;
}
```
### 3. **Critical Implementation Differences**
#### **Parameter Mapping**:
- **Fortran**: `toler0`, `toler1`, `toler2` (0.2 default)
- **Rust**: `abs_tol`, `d`, `g` (user-provided)
#### **Gradient Calculation**:
**Fortran**:
```fortran
pure real(RK) function grad(u,_x,comp,point)
grad = (u(comp,point+1)-u(comp,point)) / (_x(point+1)-_x(point))
end function grad
```
**Rust**:
```rust
fn grad(&self, y_j: DVector<f64>, i: usize) -> f64 {
let dy = y_j[i + 1] - y_j[i];
let dx = x_vector[i + 1] - x_vector[i];
dy / dx
}
```
#### **Curvature Analysis**:
**Fortran** (more sophisticated):
```fortran
! Calculate gradient at first point and store in right
right = temp
for i in 1..n-1 {
left = right;
right = list_dy_i[i];
differ_dy = (right - left).abs();
// Check against gamma threshold
}
```
**Rust** (simplified):
```rust
let temp = list_dy_i[0];
let right = temp;
for i in 1..n-1 {
let left = right;
let right = list_dy_i[i];
let differ_dy = (right - left).abs();
// Check against gamma threshold
}
```
### 4. **Weight Calculation**
**Fortran**:
```fortran
! Weight intervals based on violations
do k = 1, vars%points - 1
weight(k) = vary1(k)
if (1<k) weight(k) = weight(k) + vary2(k)
if (k<vars%points - 1) weight(k) = weight(k) + vary2(k + 1)
if (0<weight(k)) most = most + 1
end do
```
**Rust**:
```rust
for i in 0..self.n - 1 {
weights[i] = vary_y[i];
if i > 0 {
weights[i] += vary_dy_dx[i]
}
if i < self.n - 1 {
weights[i] += vary_dy_dx[i + 1]
}
if weights[i] > 0 {
most += 1
}
}
```
### 5. **Point Insertion Strategy**
**Fortran** (backward iteration):
```fortran
new_points: do old = vars%points, 2, -1
! Copy right boundary
vars%_x(new) = vars%_x(old)
u(:,new) = u(:,old)
new = new - 1
! Interpolate if marked
if (mark(old-1)) then
vars%_x(new) = half*(vars%_x(old)+vars%_x(old-1))
u(:,new) = half*(u(:,old)+u(:,old-1))
new = new - 1
end if
end do new_points
```
**Rust** (similar backward iteration):
```rust
for i_old in (1..n_former).rev() {
new_grid[n_new - 1] = x[i_old];
new_initial_guess.column_mut(n_new - 1).copy_from(&y.column(i_old));
n_new -= 1;
if mark[i_old - 1] {
new_grid[n_new - 1] = 0.5 * (x[i_old] + x[i_old - 1]);
let dy_i = y_min_1 + 0.5 * (y_i - y_min_1);
new_initial_guess.column_mut(n_new - 1).copy_from(&dy_i);
n_new -= 1;
}
}
```
## Major Issues in Rust Implementation
### 1. **Missing Active Component Check**
- **Fortran**: Only processes `vars%active(j)` components
- **Rust**: Processes all components without filtering
### 2. **Incomplete Curvature Analysis**
- **Fortran**: Proper left/right gradient comparison
- **Rust**: Simplified implementation missing gradient range calculation
### 3. **Parameter Inconsistency**
- **Fortran**: Uses relative thresholds (fraction of range)
- **Rust**: Uses absolute thresholds (delta, gamma)
### 4. **Missing Degeneracy Checks**
- **Fortran**: Extensive validation for degenerate intervals
- **Rust**: Basic monotonicity check only
### 5. **Error Handling**
- **Fortran**: Comprehensive error reporting with specific failure modes
- **Rust**: Basic logging without detailed diagnostics
## Recommendations for Rust Implementation
1. **Add Active Component Filtering**:
```rust
if !self.active_components[j] { continue; }
```
2. **Fix Curvature Calculation**:
```rust
let derivative_range = grad_max - grad_min;
let gamma = derivative_range * self.g;
```
3. **Implement Relative Thresholds**:
```rust
let delta = self.d * range_y_j; ```
4. **Add Comprehensive Validation**:
```rust
```
5. **Improve Error Reporting**:
```rust
```
## Conclusion
The Rust implementation captures the basic structure but misses several critical details from the Fortran original:
- Active component filtering
- Proper relative threshold calculation
- Complete curvature analysis
- Comprehensive error handling
The mathematical core is similar, but the implementation details significantly affect the algorithm's robustness and accuracy.
