qualia_core_db/specialized_libs/computational_economics/
welfare.rs

1//! Welfare economics primitives: social welfare functions, inequality and
2//! poverty metrics, cost-benefit analysis with distributional weights, and a
3//! needs/survival-floor allocation model.
4//!
5//! # Allocation class
6//!
7//! All kernels here are `AllocationClass::HotZeroHeap`. Sorting scratch uses
8//! fixed-capacity stack arrays (`[f64; MAX_POPULATION]`) or caller-provided
9//! output buffers — never `Vec`, `String`, or `Box`. No allocation occurs on
10//! any path, hot or cold, within this module.
11//!
12//! # Rights-affecting use
13//!
14//! Poverty metrics, distributional weights, and the survival-floor allocation
15//! model can inform decisions that affect entitlements, transfers, or access.
16//! Those uses MUST be paired with SHACL/deontic checks (see `deontic_logic.rs`
17//! opcodes `OP_OBLIGATE`/`OP_FORBID`/`OP_PERMIT`) and capacity-modalities
18//! review before any UI exposure or downstream action. The `repr(C)` report
19//! structs returned by rights-affecting functions carry diagnostics and
20//! assumptions, not just a scalar, so callers can audit the basis of a
21//! computed allocation. No function in this module performs or recommends an
22//! external action on its own.
23
24/// Maximum population size supported by stack-array scratch buffers.
25pub const MAX_POPULATION: usize = 256;
26
27/// Welfare kernel error vocabulary.
28///
29/// `repr(u8)` for ABI-stable reporting across WASM / edge / Webizen dispatch.
30#[derive(Debug, Clone, Copy, PartialEq, Eq)]
31#[repr(u8)]
32pub enum WelfareError {
33    /// Input failed validation (bad dimensions, negative where non-negative
34    /// is required, out-of-range parameter, mismatched lengths, etc.).
35    InvalidInput = 0,
36    /// Population or series is empty where a non-empty population is required.
37    InsufficientData = 1,
38    /// A non-finite value (NaN/inf) was encountered in input or during
39    /// computation.
40    NonFinite = 2,
41    /// A caller-supplied output buffer was too small for the request.
42    BufferTooSmall = 3,
43}
44
45impl WelfareError {
46    /// True when the failure is a caller-side buffer/dimension problem.
47    #[inline]
48    pub fn is_caller_error(self) -> bool {
49        matches!(
50            self,
51            WelfareError::InvalidInput
52                | WelfareError::BufferTooSmall
53                | WelfareError::InsufficientData
54        )
55    }
56}
57
58/// Report returned by rights-affecting welfare kernels.
59///
60/// `repr(C)` so it can cross the WASM / edge / GPU ABI as a fixed record and
61/// be audited by deontic / SHACL layers before any UI exposure. The scalar
62/// `value` is never the whole story: `assumptions` and `diagnostics` carry the
63/// basis of the computation.
64#[derive(Debug, Clone, Copy, PartialEq)]
65#[repr(C)]
66pub struct WelfareReport {
67    /// The primary scalar result (count, ratio, allocated amount, etc.).
68    /// Semantics are documented at each call site.
69    pub value: f64,
70    /// Secondary scalar (e.g. population size as f64, or a companion metric).
71    pub auxiliary: f64,
72    /// Bit-packed assumption flags (see `ASSUMPTION_*` constants).
73    pub assumptions: u32,
74    /// Diagnostic status code (0 = clean; non-zero mirrors `WelfareError` as
75    /// u8 plus kernel-specific high bits reserved for future use).
76    pub diagnostics: u32,
77}
78
79/// Assumption flag: at least one observation was clamped to the survival floor.
80pub const ASSUMPTION_FLOOR_CLAMPED: u32 = 1 << 0;
81/// Assumption flag: distributional weights were applied (CBA only).
82pub const ASSUMPTION_WEIGHTED: u32 = 1 << 1;
83/// Assumption flag: result is degenerate (e.g. zero-mean population for Gini).
84pub const ASSUMPTION_DEGENERATE: u32 = 1 << 2;
85/// Assumption flag: poverty line is above the maximum observed income.
86pub const ASSUMPTION_LINE_ABOVE_MAX: u32 = 1 << 3;
87
88impl WelfareReport {
89    /// Construct a clean report with no assumption flags.
90    pub const fn clean(value: f64, auxiliary: f64) -> Self {
91        Self {
92            value,
93            auxiliary,
94            assumptions: 0,
95            diagnostics: 0,
96        }
97    }
98
99    /// Construct a report with assumption flags set.
100    pub const fn with_assumptions(value: f64, auxiliary: f64, assumptions: u32) -> Self {
101        Self {
102            value,
103            auxiliary,
104            assumptions,
105            diagnostics: 0,
106        }
107    }
108}
109
110// ---------------------------------------------------------------------------
111// Validation helpers
112// ---------------------------------------------------------------------------
113
114#[inline]
115fn finite(x: f64) -> bool {
116    x.is_finite()
117}
118
119#[inline]
120fn finite_nonnegative(x: f64) -> bool {
121    x.is_finite() && x >= 0.0
122}
123
124#[inline]
125fn finite_strictly_positive(x: f64) -> bool {
126    x.is_finite() && x > 0.0
127}
128
129/// Validate a population slice: non-empty, within capacity, all finite and
130/// non-negative. Returns the population count on success.
131fn validate_population(incomes: &[f64]) -> Result<usize, WelfareError> {
132    if incomes.is_empty() {
133        return Err(WelfareError::InsufficientData);
134    }
135    if incomes.len() > MAX_POPULATION {
136        return Err(WelfareError::BufferTooSmall);
137    }
138    for &x in incomes {
139        if !finite_nonnegative(x) {
140            if !finite(x) {
141                return Err(WelfareError::NonFinite);
142            }
143            return Err(WelfareError::InvalidInput);
144        }
145    }
146    Ok(incomes.len())
147}
148
149/// Validate a population slice where strictly positive values are required
150/// (e.g. geometric-mean based Atkinson index).
151fn validate_positive_population(incomes: &[f64]) -> Result<usize, WelfareError> {
152    if incomes.is_empty() {
153        return Err(WelfareError::InsufficientData);
154    }
155    if incomes.len() > MAX_POPULATION {
156        return Err(WelfareError::BufferTooSmall);
157    }
158    for &x in incomes {
159        if !finite(x) {
160            return Err(WelfareError::NonFinite);
161        }
162        if x <= 0.0 {
163            return Err(WelfareError::InvalidInput);
164        }
165    }
166    Ok(incomes.len())
167}
168
169/// Copy `incomes` into a stack scratch array and sort ascending. Returns the
170/// count (== `incomes.len()`). Zero-heap: scratch is a fixed `[f64; MAX_POPULATION]`.
171fn copy_sort_ascending(incomes: &[f64]) -> Result<[f64; MAX_POPULATION], WelfareError> {
172    let n = validate_population(incomes)?;
173    let mut scratch = [0.0f64; MAX_POPULATION];
174    for (i, &x) in incomes.iter().enumerate() {
175        scratch[i] = x;
176    }
177    scratch[..n].sort_by(|a, b| a.total_cmp(b));
178    Ok(scratch)
179}
180
181// ---------------------------------------------------------------------------
182// Inequality metrics
183// ---------------------------------------------------------------------------
184
185/// Gini coefficient in `[0, 1]` via the mean-absolute-difference formula:
186///
187/// `G = sum_i sum_j |x_i - x_j| / (2 * n^2 * mean)`
188///
189/// Returns `WelfareError::InvalidInput` when the mean is zero (degenerate
190/// population of all zeros) since the Gini is undefined there. Inputs must be
191/// non-negative and finite; an empty population returns `InsufficientData`.
192///
193/// Allocation class: `HotZeroHeap`. The O(n^2) mean-absolute-difference sum
194/// is computed in place over the input slice — no sort is needed for this
195/// formula, so no scratch is used.
196pub fn gini_coefficient(incomes: &[f64]) -> Result<f64, WelfareError> {
197    let n = validate_population(incomes)?;
198    let mut sum: f64 = 0.0;
199    let mut total: f64 = 0.0;
200    for i in 0..n {
201        total += incomes[i];
202        for j in 0..n {
203            let d = incomes[i] - incomes[j];
204            sum += d.abs();
205        }
206    }
207    if total <= 0.0 {
208        // All-zero (or all-degenerate) population: Gini undefined.
209        return Err(WelfareError::InvalidInput);
210    }
211    let mean = total / n as f64;
212    let denom = 2.0 * (n as f64).powi(2) * mean;
213    Ok(sum / denom)
214}
215
216/// Write the Lorenz curve as interleaved `(cumulative_population_share,
217/// cumulative_income_share)` pairs into `out`: `[pop0, inc0, pop1, inc1, ...]`,
218/// sorted ascending by income. Returns the number of points written.
219///
220/// Each point corresponds to one observation (no binning), so `out` must hold
221/// at least `2 * n` `f64` slots. The first point is `(1/n, x_min/total)` and
222/// the last is `(1, 1)`.
223///
224/// Allocation class: `HotZeroHeap`. Sorting uses a fixed `[f64; MAX_POPULATION]`
225/// stack scratch array.
226pub fn lorenz_curve_into(incomes: &[f64], out: &mut [f64]) -> Result<usize, WelfareError> {
227    let n = validate_population(incomes)?;
228    if out.len() < 2 * n {
229        return Err(WelfareError::BufferTooSmall);
230    }
231    let sorted = copy_sort_ascending(incomes)?;
232    let mut total = 0.0f64;
233    for i in 0..n {
234        total += sorted[i];
235    }
236    if total <= 0.0 {
237        return Err(WelfareError::InvalidInput);
238    }
239    let mut cum = 0.0f64;
240    for i in 0..n {
241        cum += sorted[i];
242        out[2 * i] = (i + 1) as f64 / n as f64;
243        out[2 * i + 1] = cum / total;
244    }
245    Ok(n)
246}
247
248/// Atkinson inequality index with inequality-aversion parameter `epsilon > 0`.
249///
250/// - For `epsilon != 1`: `A = 1 - (geometric_mean / arithmetic_mean)^(1-epsilon)`.
251/// - For `epsilon == 1`: `A = 1 - geometric_mean / arithmetic_mean`.
252///
253/// All incomes must be strictly positive (the geometric mean is undefined for
254/// zero/negative values). `epsilon` must be finite and `> 0`. Returns a value
255/// in `[0, 1)`: `0` under perfect equality, increasing toward `1` as
256/// inequality rises and as `epsilon` (inequality aversion) rises.
257///
258/// Allocation class: `HotZeroHeap`.
259pub fn atkinson_inequality(incomes: &[f64], epsilon: f64) -> Result<f64, WelfareError> {
260    let n = validate_positive_population(incomes)?;
261    if !finite(epsilon) || epsilon <= 0.0 {
262        return Err(WelfareError::InvalidInput);
263    }
264
265    let mut sum = 0.0f64;
266    for i in 0..n {
267        if incomes[i] <= 0.0 {
268            return Err(WelfareError::InvalidInput);
269        }
270        sum += incomes[i];
271    }
272    let arithmetic_mean = sum / n as f64;
273    if arithmetic_mean <= 0.0 {
274        return Err(WelfareError::InvalidInput);
275    }
276
277    if (epsilon - 1.0).abs() < f64::EPSILON {
278        // Limit case: 1 - geo / arith
279        let mut log_sum = 0.0f64;
280        for i in 0..n {
281            log_sum += incomes[i].ln();
282        }
283        let geo = (log_sum / n as f64).exp();
284        let a = 1.0 - (geo / arithmetic_mean);
285        return Ok(a.clamp(0.0, 1.0));
286    }
287
288    // General case: 1 - ( power mean of order (1-eps) / arith )
289    let one_minus_eps = 1.0 - epsilon;
290    let mut sum_pow = 0.0f64;
291    for i in 0..n {
292        sum_pow += incomes[i].powf(one_minus_eps);
293    }
294    let mean_pow = sum_pow / n as f64;
295    let power_mean = if one_minus_eps.abs() > 1e-12 {
296        mean_pow.powf(1.0 / one_minus_eps)
297    } else {
298        arithmetic_mean
299    };
300    let a = 1.0 - (power_mean / arithmetic_mean);
301    if !a.is_finite() {
302        return Err(WelfareError::NonFinite);
303    }
304    Ok(a.clamp(0.0, 1.0))
305}
306
307/// Headcount poverty: returns `(count_poor, headcount_ratio)` where
308/// `count_poor` is the number of observations strictly below `poverty_line`
309/// and `headcount_ratio = count_poor / n`.
310///
311/// `poverty_line` must be finite and strictly positive. Incomes must be
312/// non-negative and finite.
313///
314/// Allocation class: `HotZeroHeap`.
315pub fn headcount_poverty(incomes: &[f64], poverty_line: f64) -> Result<(usize, f64), WelfareError> {
316    let n = validate_population(incomes)?;
317    if !finite_strictly_positive(poverty_line) {
318        return Err(WelfareError::InvalidInput);
319    }
320    let mut count = 0usize;
321    for &x in incomes {
322        if x < poverty_line {
323            count += 1;
324        }
325    }
326    let ratio = count as f64 / n as f64;
327    Ok((count, ratio))
328}
329
330/// Poverty gap ratio: `sum(max(0, line - income)) / (n * line)`.
331///
332/// Returns a value in `[0, 1]`: `0` when nobody is below the line, `1` when
333/// every observation is zero. `poverty_line` must be finite and strictly
334/// positive. Incomes must be non-negative and finite.
335///
336/// Allocation class: `HotZeroHeap`.
337pub fn poverty_gap_ratio(incomes: &[f64], poverty_line: f64) -> Result<f64, WelfareError> {
338    let n = validate_population(incomes)?;
339    if !finite_strictly_positive(poverty_line) {
340        return Err(WelfareError::InvalidInput);
341    }
342    let mut gap = 0.0f64;
343    for &x in incomes {
344        if x < poverty_line {
345            gap += poverty_line - x;
346        }
347    }
348    Ok(gap / (n as f64 * poverty_line))
349}
350
351// ---------------------------------------------------------------------------
352// Social welfare functions
353// ---------------------------------------------------------------------------
354
355/// Utilitarian (sum) social welfare: `sum_i u_i`.
356///
357/// Utilities must be finite. Negative utilities are permitted (welfare can be
358/// negative). An empty population returns `InsufficientData`.
359///
360/// Allocation class: `HotZeroHeap`.
361pub fn utilitarian_welfare(utilities: &[f64]) -> Result<f64, WelfareError> {
362    if utilities.is_empty() {
363        return Err(WelfareError::InsufficientData);
364    }
365    if utilities.len() > MAX_POPULATION {
366        return Err(WelfareError::BufferTooSmall);
367    }
368    let mut sum = 0.0f64;
369    for &u in utilities {
370        if !finite(u) {
371            return Err(WelfareError::NonFinite);
372        }
373        sum += u;
374    }
375    Ok(sum)
376}
377
378/// Rawlsian (minimax / maximin) social welfare: `min_i u_i`.
379///
380/// The welfare of a society is the welfare of its worst-off member. Utilities
381/// must be finite. An empty population returns `InsufficientData`.
382///
383/// Allocation class: `HotZeroHeap`.
384pub fn rawlsian_welfare(utilities: &[f64]) -> Result<f64, WelfareError> {
385    if utilities.is_empty() {
386        return Err(WelfareError::InsufficientData);
387    }
388    if utilities.len() > MAX_POPULATION {
389        return Err(WelfareError::BufferTooSmall);
390    }
391    let mut min = utilities[0];
392    if !finite(min) {
393        return Err(WelfareError::NonFinite);
394    }
395    for &u in &utilities[1..] {
396        if !finite(u) {
397            return Err(WelfareError::NonFinite);
398        }
399        if u < min {
400            min = u;
401        }
402    }
403    Ok(min)
404}
405
406/// Nash social welfare: the **product** of utilities, `prod_i u_i`.
407///
408/// This is the unnormalised Nash product (the sum of logs is the log of the
409/// product). It is *not* the geometric mean; divide by `n` externally if the
410/// geometric mean is required. All utilities must be strictly positive so the
411/// product is well-defined and non-degenerate; a single zero utility would
412/// collapse the product to zero and is rejected as `InvalidInput`. An empty
413/// population returns `InsufficientData`.
414///
415/// Allocation class: `HotZeroHeap`.
416pub fn nash_welfare(utilities: &[f64]) -> Result<f64, WelfareError> {
417    let n = validate_positive_population(utilities)?;
418    let mut product = 1.0f64;
419    for i in 0..n {
420        product *= utilities[i];
421    }
422    if !product.is_finite() {
423        return Err(WelfareError::NonFinite);
424    }
425    Ok(product)
426}
427
428// ---------------------------------------------------------------------------
429// Cost-benefit analysis
430// ---------------------------------------------------------------------------
431
432/// Net present value: `NPV = sum_{t=0}^{n_periods-1} (B_t - C_t) / (1+r)^t`.
433///
434/// `benefits` and `costs` must each have length `>= n_periods`; only the first
435/// `n_periods` entries are consumed. `discount_rate` must be finite and `> -1`
436/// (a rate of `-1` or below makes the discount factor non-positive). The
437/// period-0 cash flow is discounted by `(1+r)^0 = 1` (i.e. not discounted).
438///
439/// `n_periods` must be `> 0` and `<= MAX_POPULATION`.
440///
441/// Allocation class: `HotZeroHeap`.
442pub fn net_present_value(
443    benefits: &[f64],
444    costs: &[f64],
445    discount_rate: f64,
446    n_periods: usize,
447) -> Result<f64, WelfareError> {
448    if n_periods == 0 {
449        return Err(WelfareError::InsufficientData);
450    }
451    if n_periods > MAX_POPULATION {
452        return Err(WelfareError::BufferTooSmall);
453    }
454    if benefits.len() < n_periods || costs.len() < n_periods {
455        return Err(WelfareError::InvalidInput);
456    }
457    if !finite(discount_rate) || discount_rate <= -1.0 {
458        return Err(WelfareError::InvalidInput);
459    }
460    let one_plus_r = 1.0 + discount_rate;
461    let mut npv = 0.0f64;
462    // Documented convention: NPV = Σ_{t=0}^{n-1} (B_t − C_t)/(1+r)^t, so period 0
463    // is undiscounted ((1+r)^0 = 1) and each later period divides by another
464    // (1+r). (The prior code *multiplied* by (1+r) each period, computing a
465    // future value — the opposite of a present value.)
466    let mut discount = 1.0f64; // (1+r)^0
467    for t in 0..n_periods {
468        let b = benefits[t];
469        let c = costs[t];
470        if !finite(b) || !finite(c) {
471            return Err(WelfareError::NonFinite);
472        }
473        npv += (b - c) * discount;
474        discount /= one_plus_r;
475    }
476    if !npv.is_finite() {
477        return Err(WelfareError::NonFinite);
478    }
479    Ok(npv)
480}
481
482/// Distributional NPV: NPV with per-period distributional weights.
483///
484/// `NPV_w = sum_{t=0}^{n_periods-1} w_t * (B_t - C_t) / (1+r)^t`.
485///
486/// `weights` must have length `>= n_periods`. Weights must be finite and
487/// non-negative (a zero weight simply zeroes that period's contribution).
488/// `benefits`, `costs`, `discount_rate`, and `n_periods` follow the same
489/// rules as [`net_present_value`].
490///
491/// Returns a [`WelfareReport`] carrying the `ASSUMPTION_WEIGHTED` flag so
492/// downstream deontic / SHACL auditors can see that distributional weights
493/// were applied. `value` is the weighted NPV; `auxiliary` is the unweighted
494/// NPV for comparison.
495///
496/// Rights-affecting: distributional weights encode value judgements about
497/// whose benefits count how much. Pair with deontic / SHACL review before UI
498/// exposure.
499///
500/// Allocation class: `HotZeroHeap`.
501pub fn distributional_npv(
502    benefits: &[f64],
503    costs: &[f64],
504    weights: &[f64],
505    discount_rate: f64,
506    n_periods: usize,
507) -> Result<WelfareReport, WelfareError> {
508    if n_periods == 0 {
509        return Err(WelfareError::InsufficientData);
510    }
511    if n_periods > MAX_POPULATION {
512        return Err(WelfareError::BufferTooSmall);
513    }
514    if benefits.len() < n_periods || costs.len() < n_periods || weights.len() < n_periods {
515        return Err(WelfareError::InvalidInput);
516    }
517    if !finite(discount_rate) || discount_rate <= -1.0 {
518        return Err(WelfareError::InvalidInput);
519    }
520    let one_plus_r = 1.0 + discount_rate;
521    let mut weighted_npv = 0.0f64;
522    let mut unweighted_npv = 0.0f64;
523    let mut discount = 1.0f64;
524    for t in 0..n_periods {
525        let b = benefits[t];
526        let c = costs[t];
527        let w = weights[t];
528        if !finite(b) || !finite(c) || !finite(w) {
529            return Err(WelfareError::NonFinite);
530        }
531        if w < 0.0 {
532            return Err(WelfareError::InvalidInput);
533        }
534        let flow = (b - c) * discount;
535        weighted_npv += w * flow;
536        unweighted_npv += flow;
537        discount /= one_plus_r;
538    }
539    if !weighted_npv.is_finite() || !unweighted_npv.is_finite() {
540        return Err(WelfareError::NonFinite);
541    }
542    Ok(WelfareReport {
543        value: weighted_npv,
544        auxiliary: unweighted_npv,
545        assumptions: ASSUMPTION_WEIGHTED,
546        diagnostics: 0,
547    })
548}
549
550// ---------------------------------------------------------------------------
551// Needs / survival-floor allocation model
552// ---------------------------------------------------------------------------
553
554/// Allocate a fixed budget `budget` across `needs` so that each recipient
555/// first receives their survival floor `floors[i]`, then the residual is
556/// distributed proportionally to surplus need `needs[i] - floors[i]`.
557///
558/// Composes with deontic and capacity modalities: the returned
559/// [`WelfareReport`] carries assumption flags so a downstream deontic / SHACL
560/// layer can verify that floors were honoured (`ASSUMPTION_FLOOR_CLAMPED` is
561/// set when any allocation was clamped to the floor) and that the budget was
562/// sufficient (`ASSUMPTION_DEGENERATE` is set when the budget could not cover
563/// all floors, in which case floors are scaled proportionally).
564///
565/// `value` is the total actually allocated; `auxiliary` is the residual after
566/// floors (the amount distributed proportionally). Results are written into
567/// the caller-owned `out` buffer (`out.len() >= n`).
568///
569/// Rights-affecting: this model can determine access to subsistence resources.
570/// Pair with `OP_OBLIGATE`/`OP_FORBID` deontic checks and capacity-modalities
571/// review before any UI exposure or downstream transfer.
572///
573/// Allocation class: `HotZeroHeap`.
574pub fn survival_floor_allocation_into(
575    needs: &[f64],
576    floors: &[f64],
577    budget: f64,
578    out: &mut [f64],
579) -> Result<WelfareReport, WelfareError> {
580    if needs.is_empty() {
581        return Err(WelfareError::InsufficientData);
582    }
583    let n = needs.len();
584    if n > MAX_POPULATION {
585        return Err(WelfareError::BufferTooSmall);
586    }
587    if floors.len() != n || out.len() < n {
588        return Err(WelfareError::InvalidInput);
589    }
590    if !finite(budget) || budget < 0.0 {
591        return Err(WelfareError::InvalidInput);
592    }
593
594    let mut total_floor = 0.0f64;
595    let mut total_surplus_need = 0.0f64;
596    for i in 0..n {
597        if !finite(needs[i]) || !finite(floors[i]) {
598            return Err(WelfareError::NonFinite);
599        }
600        if needs[i] < 0.0 || floors[i] < 0.0 {
601            return Err(WelfareError::InvalidInput);
602        }
603        if floors[i] > needs[i] {
604            // Floor cannot exceed total need.
605            return Err(WelfareError::InvalidInput);
606        }
607        total_floor += floors[i];
608        total_surplus_need += needs[i] - floors[i];
609    }
610
611    let mut assumptions: u32 = 0;
612    let mut total_allocated = 0.0f64;
613    let mut residual_distributed = 0.0f64;
614
615    if budget >= total_floor {
616        // Floors fully covered; distribute residual proportionally to surplus need.
617        let residual = budget - total_floor;
618        for i in 0..n {
619            let base = floors[i];
620            let extra = if total_surplus_need > 0.0 {
621                residual * (needs[i] - floors[i]) / total_surplus_need
622            } else {
623                0.0
624            };
625            let alloc = base + extra;
626            out[i] = alloc;
627            total_allocated += alloc;
628            residual_distributed += extra;
629        }
630    } else {
631        // Insufficient budget: scale floors proportionally so no recipient is
632        // arbitrarily zeroed out. This is a degenerate (under-budget) case.
633        assumptions |= ASSUMPTION_DEGENERATE;
634        let scale = if total_floor > 0.0 {
635            budget / total_floor
636        } else {
637            0.0
638        };
639        for i in 0..n {
640            let alloc = floors[i] * scale;
641            out[i] = alloc;
642            total_allocated += alloc;
643        }
644        assumptions |= ASSUMPTION_FLOOR_CLAMPED;
645    }
646
647    if !total_allocated.is_finite() || !residual_distributed.is_finite() {
648        return Err(WelfareError::NonFinite);
649    }
650
651    Ok(WelfareReport {
652        value: total_allocated,
653        auxiliary: residual_distributed,
654        assumptions,
655        diagnostics: 0,
656    })
657}
658
659// ---------------------------------------------------------------------------
660// Tests
661// ---------------------------------------------------------------------------
662
663#[cfg(test)]
664mod tests {
665    use super::*;
666
667    fn approx(a: f64, b: f64) -> bool {
668        (a - b).abs() < 1e-9
669    }
670
671    // --- Gini ---------------------------------------------------------------
672
673    #[test]
674    fn gini_perfect_equality_is_zero() {
675        let incomes = [10.0, 10.0, 10.0];
676        let g = gini_coefficient(&incomes).unwrap();
677        assert!(approx(g, 0.0));
678    }
679
680    #[test]
681    fn gini_extreme_inequality_approaches_one_minus_one_over_n() {
682        // [0, 0, 100]: G = 1 - 1/n = 2/3 for n = 3.
683        let incomes = [0.0, 0.0, 100.0];
684        let g = gini_coefficient(&incomes).unwrap();
685        let expected = 1.0 - 1.0 / 3.0;
686        assert!(approx(g, expected));
687    }
688
689    #[test]
690    fn gini_two_person_split_is_one_half() {
691        // [0, 1]: G = 1 - 1/2 = 0.5.
692        let incomes = [0.0, 1.0];
693        let g = gini_coefficient(&incomes).unwrap();
694        assert!(approx(g, 0.5));
695    }
696
697    #[test]
698    fn gini_empty_is_insufficient_data() {
699        assert_eq!(gini_coefficient(&[]), Err(WelfareError::InsufficientData));
700    }
701
702    #[test]
703    fn gini_all_zero_is_invalid() {
704        assert_eq!(
705            gini_coefficient(&[0.0, 0.0]),
706            Err(WelfareError::InvalidInput)
707        );
708    }
709
710    #[test]
711    fn gini_nan_is_non_finite() {
712        assert_eq!(
713            gini_coefficient(&[1.0, f64::NAN]),
714            Err(WelfareError::NonFinite)
715        );
716    }
717
718    #[test]
719    fn gini_negative_is_invalid() {
720        assert_eq!(
721            gini_coefficient(&[1.0, -1.0]),
722            Err(WelfareError::InvalidInput)
723        );
724    }
725
726    // --- Lorenz -------------------------------------------------------------
727
728    #[test]
729    fn lorenz_curve_cumulative_shares() {
730        let incomes = [1.0, 2.0, 3.0];
731        let mut out = [0.0f64; 6];
732        let n = lorenz_curve_into(&incomes, &mut out).unwrap();
733        assert_eq!(n, 3);
734        // Sorted: [1, 2, 3], total = 6.
735        // Points: (1/3, 1/6), (2/3, 3/6=1/2), (3/3=1, 6/6=1).
736        assert!(approx(out[0], 1.0 / 3.0));
737        assert!(approx(out[1], 1.0 / 6.0));
738        assert!(approx(out[2], 2.0 / 3.0));
739        assert!(approx(out[3], 0.5));
740        assert!(approx(out[4], 1.0));
741        assert!(approx(out[5], 1.0));
742    }
743
744    #[test]
745    fn lorenz_buffer_too_small() {
746        let incomes = [1.0, 2.0, 3.0];
747        let mut out = [0.0f64; 5]; // need 6
748        assert_eq!(
749            lorenz_curve_into(&incomes, &mut out),
750            Err(WelfareError::BufferTooSmall)
751        );
752    }
753
754    #[test]
755    fn lorenz_empty_is_insufficient_data() {
756        let mut out = [0.0f64; 2];
757        assert_eq!(
758            lorenz_curve_into(&[], &mut out),
759            Err(WelfareError::InsufficientData)
760        );
761    }
762
763    // --- Atkinson -----------------------------------------------------------
764
765    #[test]
766    fn atkinson_equality_is_zero() {
767        let incomes = [10.0, 10.0, 10.0];
768        let a = atkinson_inequality(&incomes, 0.5).unwrap();
769        assert!(approx(a, 0.0));
770    }
771
772    #[test]
773    fn atkinson_higher_epsilon_is_more_inequality_averse() {
774        let incomes = [1.0, 2.0, 10.0];
775        let a_low = atkinson_inequality(&incomes, 0.5).unwrap();
776        let a_high = atkinson_inequality(&incomes, 2.0).unwrap();
777        assert!(a_high > a_low, "higher epsilon must yield higher A");
778    }
779
780    #[test]
781    fn atkinson_epsilon_one_uses_log_form() {
782        let incomes = [1.0, 2.0, 4.0];
783        let a = atkinson_inequality(&incomes, 1.0).unwrap();
784        // geometric mean = (1*2*4)^(1/3) = 8^(1/3) = 2
785        // arithmetic mean = 7/3
786        // A = 1 - 2 / (7/3) = 1 - 6/7 = 1/7
787        assert!(approx(a, 1.0 / 7.0));
788    }
789
790    #[test]
791    fn atkinson_zero_income_is_invalid() {
792        assert_eq!(
793            atkinson_inequality(&[0.0, 1.0], 0.5),
794            Err(WelfareError::InvalidInput)
795        );
796    }
797
798    #[test]
799    fn atkinson_nonpositive_epsilon_is_invalid() {
800        assert_eq!(
801            atkinson_inequality(&[1.0, 2.0], 0.0),
802            Err(WelfareError::InvalidInput)
803        );
804        assert_eq!(
805            atkinson_inequality(&[1.0, 2.0], -1.0),
806            Err(WelfareError::InvalidInput)
807        );
808    }
809
810    // --- Headcount ----------------------------------------------------------
811
812    #[test]
813    fn headcount_counts_poor_and_ratio() {
814        let incomes = [5.0, 15.0, 25.0];
815        let (count, ratio) = headcount_poverty(&incomes, 10.0).unwrap();
816        assert_eq!(count, 1);
817        assert!(approx(ratio, 1.0 / 3.0));
818    }
819
820    #[test]
821    fn headcount_none_poor() {
822        let incomes = [20.0, 30.0];
823        let (count, ratio) = headcount_poverty(&incomes, 10.0).unwrap();
824        assert_eq!(count, 0);
825        assert!(approx(ratio, 0.0));
826    }
827
828    #[test]
829    fn headcount_line_at_income_excludes_boundary() {
830        // Strictly below the line: an income equal to the line is not poor.
831        let incomes = [10.0, 20.0];
832        let (count, _) = headcount_poverty(&incomes, 10.0).unwrap();
833        assert_eq!(count, 0);
834    }
835
836    #[test]
837    fn headcount_zero_line_is_invalid() {
838        assert_eq!(
839            headcount_poverty(&[1.0, 2.0], 0.0),
840            Err(WelfareError::InvalidInput)
841        );
842    }
843
844    // --- Poverty gap --------------------------------------------------------
845
846    #[test]
847    fn poverty_gap_ratio_basic() {
848        let incomes = [5.0, 15.0, 25.0];
849        let g = poverty_gap_ratio(&incomes, 10.0).unwrap();
850        // gap = 5; n*line = 30; 5/30 = 1/6.
851        assert!(approx(g, 1.0 / 6.0));
852    }
853
854    #[test]
855    fn poverty_gap_none_poor_is_zero() {
856        let incomes = [20.0, 30.0];
857        let g = poverty_gap_ratio(&incomes, 10.0).unwrap();
858        assert!(approx(g, 0.0));
859    }
860
861    #[test]
862    fn poverty_gap_all_zero_is_one() {
863        let incomes = [0.0, 0.0];
864        let g = poverty_gap_ratio(&incomes, 10.0).unwrap();
865        assert!(approx(g, 1.0));
866    }
867
868    // --- Utilitarian / Rawlsian / Nash --------------------------------------
869
870    #[test]
871    fn utilitarian_is_sum() {
872        let u = [1.0, 2.0, 3.0];
873        assert!(approx(utilitarian_welfare(&u).unwrap(), 6.0));
874    }
875
876    #[test]
877    fn utilitarian_allows_negative() {
878        let u = [-1.0, 2.0, 3.0];
879        assert!(approx(utilitarian_welfare(&u).unwrap(), 4.0));
880    }
881
882    #[test]
883    fn rawlsian_is_min() {
884        let u = [1.0, 2.0, 3.0];
885        assert!(approx(rawlsian_welfare(&u).unwrap(), 1.0));
886    }
887
888    #[test]
889    fn rawlsian_negative_min() {
890        let u = [-5.0, 2.0, 3.0];
891        assert!(approx(rawlsian_welfare(&u).unwrap(), -5.0));
892    }
893
894    #[test]
895    fn nash_is_product() {
896        let u = [1.0, 2.0, 3.0];
897        assert!(approx(nash_welfare(&u).unwrap(), 6.0));
898    }
899
900    #[test]
901    fn nash_zero_utility_is_invalid() {
902        assert_eq!(nash_welfare(&[0.0, 2.0]), Err(WelfareError::InvalidInput));
903    }
904
905    #[test]
906    fn nash_negative_utility_is_invalid() {
907        assert_eq!(nash_welfare(&[-1.0, 2.0]), Err(WelfareError::InvalidInput));
908    }
909
910    #[test]
911    fn welfare_empty_is_insufficient_data() {
912        assert_eq!(
913            utilitarian_welfare(&[]),
914            Err(WelfareError::InsufficientData)
915        );
916        assert_eq!(rawlsian_welfare(&[]), Err(WelfareError::InsufficientData));
917        assert_eq!(nash_welfare(&[]), Err(WelfareError::InsufficientData));
918    }
919
920    #[test]
921    fn welfare_nan_is_non_finite() {
922        assert_eq!(
923            utilitarian_welfare(&[1.0, f64::NAN]),
924            Err(WelfareError::NonFinite)
925        );
926        assert_eq!(
927            rawlsian_welfare(&[1.0, f64::NAN]),
928            Err(WelfareError::NonFinite)
929        );
930    }
931
932    // --- NPV ----------------------------------------------------------------
933
934    #[test]
935    fn npv_basic_discounting() {
936        let benefits = [10.0, 10.0];
937        let costs = [5.0, 5.0];
938        let r = 0.1;
939        let npv = net_present_value(&benefits, &costs, r, 2).unwrap();
940        // Documented convention (period 0 undiscounted): 5/(1.1)^0 + 5/(1.1)^1.
941        // (The prior expected value 5/1.1 + 5/1.21 discounted period 0 too, which
942        // contradicts both the doc comment and `npv_period_zero_not_discounted`.)
943        let expected = 5.0 + 5.0 / 1.1;
944        assert!(approx(npv, expected));
945    }
946
947    #[test]
948    fn npv_period_zero_not_discounted() {
949        let benefits = [100.0, 0.0];
950        let costs = [0.0, 0.0];
951        let npv = net_present_value(&benefits, &costs, 0.5, 2).unwrap();
952        assert!(approx(npv, 100.0));
953    }
954
955    #[test]
956    fn npv_zero_rate_is_sum() {
957        let benefits = [10.0, 10.0, 10.0];
958        let costs = [1.0, 2.0, 3.0];
959        let npv = net_present_value(&benefits, &costs, 0.0, 3).unwrap();
960        assert!(approx(npv, 24.0));
961    }
962
963    #[test]
964    fn npv_length_mismatch_is_invalid() {
965        let benefits = [10.0];
966        let costs = [5.0, 5.0];
967        assert_eq!(
968            net_present_value(&benefits, &costs, 0.1, 2),
969            Err(WelfareError::InvalidInput)
970        );
971    }
972
973    #[test]
974    fn npv_zero_periods_is_insufficient_data() {
975        assert_eq!(
976            net_present_value(&[1.0], &[1.0], 0.1, 0),
977            Err(WelfareError::InsufficientData)
978        );
979    }
980
981    #[test]
982    fn npv_rate_at_minus_one_is_invalid() {
983        assert_eq!(
984            net_present_value(&[1.0], &[1.0], -1.0, 1),
985            Err(WelfareError::InvalidInput)
986        );
987    }
988
989    // --- Distributional NPV -------------------------------------------------
990
991    #[test]
992    fn distributional_npv_applies_weights() {
993        let benefits = [10.0, 10.0];
994        let costs = [5.0, 5.0];
995        let weights = [1.0, 2.0];
996        let r = 0.1;
997        let report = distributional_npv(&benefits, &costs, &weights, r, 2).unwrap();
998        // t=0: *1.0 , t=1: /1.1
999        let expected_weighted = 1.0 * 5.0 + 2.0 * 5.0 / 1.1;
1000        let expected_unweighted = 5.0 + 5.0 / 1.1;
1001        assert!(approx(report.value, expected_weighted));
1002        assert!(approx(report.auxiliary, expected_unweighted));
1003        assert!(report.assumptions & ASSUMPTION_WEIGHTED != 0);
1004    }
1005
1006    #[test]
1007    fn distributional_npv_negative_weight_is_invalid() {
1008        assert_eq!(
1009            distributional_npv(&[10.0], &[5.0], &[-1.0], 0.1, 1),
1010            Err(WelfareError::InvalidInput)
1011        );
1012    }
1013
1014    #[test]
1015    fn distributional_npv_weights_length_mismatch() {
1016        assert_eq!(
1017            distributional_npv(&[10.0, 10.0], &[5.0, 5.0], &[1.0], 0.1, 2),
1018            Err(WelfareError::InvalidInput)
1019        );
1020    }
1021
1022    // --- Survival-floor allocation ------------------------------------------
1023
1024    #[test]
1025    fn allocation_covers_floors_then_distributes_residual() {
1026        let needs = [10.0, 20.0, 30.0];
1027        let floors = [4.0, 5.0, 6.0];
1028        let budget = 25.0; // total_floor = 15, residual = 10
1029        let mut out = [0.0f64; 3];
1030        let report = survival_floor_allocation_into(&needs, &floors, budget, &mut out).unwrap();
1031        // surplus needs: 6, 15, 24 -> total 45
1032        // extras: 10*6/45, 10*15/45, 10*24/45
1033        let expected = [
1034            4.0 + 10.0 * 6.0 / 45.0,
1035            5.0 + 10.0 * 15.0 / 45.0,
1036            6.0 + 10.0 * 24.0 / 45.0,
1037        ];
1038        for i in 0..3 {
1039            assert!(approx(out[i], expected[i]));
1040        }
1041        assert!(approx(report.value, budget));
1042        assert!(approx(report.auxiliary, 10.0));
1043        assert_eq!(report.assumptions, 0);
1044    }
1045
1046    #[test]
1047    fn allocation_under_budget_scales_floors() {
1048        let needs = [10.0, 20.0];
1049        let floors = [4.0, 6.0];
1050        let budget = 5.0; // total_floor = 10 > budget -> degenerate
1051        let mut out = [0.0f64; 2];
1052        let report = survival_floor_allocation_into(&needs, &floors, budget, &mut out).unwrap();
1053        // scale = 5/10 = 0.5 -> [2, 3]
1054        assert!(approx(out[0], 2.0));
1055        assert!(approx(out[1], 3.0));
1056        assert!(approx(report.value, 5.0));
1057        assert!(report.assumptions & ASSUMPTION_DEGENERATE != 0);
1058        assert!(report.assumptions & ASSUMPTION_FLOOR_CLAMPED != 0);
1059    }
1060
1061    #[test]
1062    fn allocation_floor_exceeds_need_is_invalid() {
1063        assert_eq!(
1064            survival_floor_allocation_into(&[5.0], &[10.0], 100.0, &mut [0.0]),
1065            Err(WelfareError::InvalidInput)
1066        );
1067    }
1068
1069    #[test]
1070    fn allocation_buffer_too_small() {
1071        let needs = [10.0, 20.0];
1072        let floors = [4.0, 6.0];
1073        let mut out = [0.0f64; 1];
1074        assert_eq!(
1075            survival_floor_allocation_into(&needs, &floors, 100.0, &mut out),
1076            Err(WelfareError::InvalidInput)
1077        );
1078    }
1079
1080    #[test]
1081    fn allocation_empty_is_insufficient_data() {
1082        assert_eq!(
1083            survival_floor_allocation_into(&[], &[], 100.0, &mut []),
1084            Err(WelfareError::InsufficientData)
1085        );
1086    }
1087
1088    // --- Error classification -----------------------------------------------
1089
1090    #[test]
1091    fn error_is_caller_error_classification() {
1092        assert!(WelfareError::InvalidInput.is_caller_error());
1093        assert!(WelfareError::InsufficientData.is_caller_error());
1094        assert!(WelfareError::BufferTooSmall.is_caller_error());
1095        assert!(!WelfareError::NonFinite.is_caller_error());
1096    }
1097
1098    #[test]
1099    fn welfare_report_clean_constructor() {
1100        let r = WelfareReport::clean(1.0, 2.0);
1101        assert_eq!(r.value, 1.0);
1102        assert_eq!(r.auxiliary, 2.0);
1103        assert_eq!(r.assumptions, 0);
1104        assert_eq!(r.diagnostics, 0);
1105    }
1106
1107    #[test]
1108    fn welfare_report_with_assumptions_constructor() {
1109        let r = WelfareReport::with_assumptions(1.0, 2.0, ASSUMPTION_WEIGHTED);
1110        assert!(r.assumptions & ASSUMPTION_WEIGHTED != 0);
1111    }
1112
1113    #[test]
1114    fn over_capacity_population_is_buffer_too_small() {
1115        let big = [1.0f64; MAX_POPULATION + 1];
1116        // Use a small slice that exceeds capacity.
1117        let slice = &big[..MAX_POPULATION + 1];
1118        assert_eq!(gini_coefficient(slice), Err(WelfareError::BufferTooSmall));
1119    }
1120}
qualia_core_db/specialized_libs/computational_economics/welfare.rs

qualia_core_db/specialized_libs/computational_economics/
welfare.rs
