/* MSC-P-032 segmentation stability. MIT License. Secondary implementation.
   The declared linear congruential recursion is reproduced exactly, so this
   program draws the same resamples as the Python, R and SPSS references. The
   clustering itself is delegated to PROC FASTCLUS, whose assignment order,
   tie-breaking and convergence rule are its own, so the recovery values may
   differ in the last decimals. This file is an independent check of the
   stability conclusion and of the flags, not a digit-for-digit reproduction. */
%let candidates = 2 3 4;
%let proposed = 4;
%let resamples = 40;
%let recovery_min = 0.75;
%let margin_min = 0.10;
%let lcg_seed = 20260912;
%let lcg_multiplier = 1103515245;
%let lcg_increment = 12345;
%let lcg_modulus = 2147483648;
/* SAS numbers are doubles, exact for integers only up to 2^53. Written directly,
   &lcg_multiplier * state reaches 2.4e18 and would be rounded, giving a different
   recursion from the Python reference. The multiplier is therefore split; the
   split is exact, since 16838 * 65536 + 20077 is the multiplier itself. */
%let lcg_high = 16838;
%let lcg_low = 20077;

filename p032csv "public/datasets/msc-p032-stability-panel.csv";
data _null_;
  infile p032csv obs=1 lrecl=32767 truncover;
  input;
  if strip(_infile_) ne "customer_id,recency_days,frequency_12m,avg_basket_eur" then do;
    put "ERROR: exact schema required"; abort cancel;
  end;
run;
data p032;
  infile p032csv dsd firstobs=2 truncover lrecl=32767 end=eof;
  input @;
  if countc(_infile_, ',') ne 3 then do; put "ERROR: exactly four cells required per row"; abort cancel; end;
  length customer_id $8;
  input customer_id $ recency_days frequency_12m avg_basket_eur;
  if customer_id ne cats('T', put(_n_, z4.)) then do; put "ERROR: customers must be ordered T0001, T0002, ... without gaps"; abort cancel; end;
  array positive[3] recency_days frequency_12m avg_basket_eur;
  do i = 1 to 3;
    if missing(positive[i]) or positive[i] <= 0 then do; put "ERROR: recency, frequency and basket must be finite and strictly positive"; abort cancel; end;
  end;
  drop i;
  log_recency = log(recency_days);
  frequency = frequency_12m;
  log_basket = log(avg_basket_eur);
  order = _n_;
  if eof then call symputx("customers", _n_);
run;
%if &customers < 100 %then %do; %put ERROR: too few customers for the declared stability analysis; %abort cancel; %end;
/* Declared standardization: mean zero, unit standard deviation over the whole file. */
proc stdize data=p032 out=scaled method=std;
  var log_recency frequency log_basket;
run;
/* Declared structureless reference: each feature permuted independently by the
   same recursion, so a recovery level can be read against resampling noise. */
data reference;
  set scaled end=last;
  array feature[3] log_recency frequency log_basket;
  array store[3, 100000] _temporary_;
  retain state &lcg_seed;
  do j = 1 to 3; store[j, _n_] = feature[j]; end;
  if last then do;
    do j = 1 to 3;
      do i = &customers to 2 by -1;
        high = mod(&lcg_high * state, &lcg_modulus);
        state = mod(high * 65536 + &lcg_low * state + &lcg_increment, &lcg_modulus);
        swap = mod(state, i) + 1;
        hold = store[j, i]; store[j, i] = store[j, swap]; store[j, swap] = hold;
      end;
    end;
    do i = 1 to &customers;
      order = i;
      log_recency = store[1, i]; frequency = store[2, i]; log_basket = store[3, i];
      output;
    end;
  end;
  keep order log_recency frequency log_basket;
run;

%macro recovery(source=, k=, out=);
  proc fastclus data=&source maxclusters=&k maxiter=100 converge=0 out=base(keep=order cluster rename=(cluster=segment)) noprint;
    var log_recency frequency log_basket;
  run;
  data _null_; call symputx("state", &lcg_seed); run;
  proc datasets nolist nowarn; delete recoveries; quit;
  proc sql noprint; create table recoveries (segment num, jaccard num); quit;
  %do b = 1 %to &resamples;
    data draw;
      state = symgetn("state");
      do i = 1 to &customers;
        high = mod(&lcg_high * state, &lcg_modulus);
        state = mod(high * 65536 + &lcg_low * state + &lcg_increment, &lcg_modulus);
        pick = mod(state, &customers) + 1;
        output;
      end;
      call symputx("state", state);
      keep pick;
    run;
    proc sql noprint;
      create table sample as select s.* from draw d inner join &source s on d.pick = s.order;
    quit;
    proc fastclus data=sample maxclusters=&k maxiter=100 converge=0 mean=centroids noprint;
      var log_recency frequency log_basket;
    run;
    proc fastclus data=&source seed=centroids maxclusters=&k maxiter=0 out=relabelled(keep=order cluster rename=(cluster=rebuilt)) noprint;
      var log_recency frequency log_basket;
    run;
    proc sql noprint;
      create table pairs as
        select b.segment, r.rebuilt, count(*) as inter
        from base b inner join relabelled r on b.order = r.order
        group by b.segment, r.rebuilt;
      create table sizes as select segment, count(*) as n_segment from base group by segment;
      create table rsizes as select rebuilt, count(*) as n_rebuilt from relabelled group by rebuilt;
      create table best as
        select p.segment, max(p.inter / (s.n_segment + t.n_rebuilt - p.inter)) as jaccard
        from pairs p inner join sizes s on p.segment = s.segment
                     inner join rsizes t on p.rebuilt = t.rebuilt
        group by p.segment;
      insert into recoveries select segment, jaccard from best;
    quit;
  %end;
  proc means data=recoveries noprint nway;
    class segment; var jaccard;
    output out=&out(drop=_type_ _freq_) mean=mean_recovery;
  run;
%mend recovery;

%macro run_all;
  data summary; length flag $10; stop; run;
  %local index k;
  %do index = 1 %to %sysfunc(countw(&candidates));
    %let k = %scan(&candidates, &index);
    %recovery(source=scaled, k=&k, out=observed&k);
    %recovery(source=reference, k=&k, out=null&k);
    proc sql noprint;
      select min(mean_recovery) into :weakest from observed&k;
      select min(mean_recovery) into :nullweak from null&k;
    quit;
    data one;
      length flag $10;
      segments = &k;
      weakest = &weakest;
      null_weakest = &nullweak;
      margin = weakest - null_weakest;
      flag = ifc(weakest >= &recovery_min and margin >= &margin_min, "STABLE", "NOT_STABLE");
    run;
    proc append base=summary data=one force; run;
  %end;
  proc sql noprint;
    select flag into :proposed_flag from summary where segments = &proposed;
  quit;
  data verdict;
    length verdict $32;
    verdict = ifc("&proposed_flag" = "STABLE", "PROPOSED_SEGMENTATION_STABLE", "PROPOSED_SEGMENTATION_NOT_STABLE");
  run;
  proc print data=summary noobs; run;
  proc print data=verdict noobs; run;
%mend run_all;
%run_all;