/* MSC-P-031 customer segmentation. MIT License. Secondary implementation.
   PROC FASTCLUS is seeded with the same four declared observations, but its
   assignment order, its tie-breaking and its convergence rule are its own, so its
   last decimals may differ from the Python, R and SPSS references. This file is an
   independent check of the segmentation, not a digit-for-digit reproduction. */
filename p031csv "public/datasets/msc-p031-segmentation-panel.csv";
data _null_;
  infile p031csv obs=1 lrecl=32767 truncover;
  input;
  if strip(_infile_) ne "customer_id,recency_days,frequency_12m,avg_basket_eur,category_breadth,digital_share,next_quarter_revenue_eur" then do;
    put "ERROR: exact schema required"; abort cancel;
  end;
run;
data p031;
  infile p031csv dsd firstobs=2 truncover lrecl=32767 end=eof;
  input @;
  if countc(_infile_, ',') ne 6 then do; put "ERROR: exactly seven cells required per row"; abort cancel; end;
  length customer_id $8;
  input customer_id $ recency_days frequency_12m avg_basket_eur category_breadth digital_share next_quarter_revenue_eur;
  if customer_id ne cats('S', put(_n_, z4.)) then do; put "ERROR: customers must be ordered S0001, S0002, ... without gaps"; abort cancel; end;
  array positive[4] recency_days frequency_12m avg_basket_eur category_breadth;
  do i = 1 to 4;
    if missing(positive[i]) or positive[i] <= 0 then do; put "ERROR: recency, frequency, basket and breadth must be finite and strictly positive"; abort cancel; end;
  end;
  drop i;
  if missing(digital_share) or digital_share <= 0 or digital_share >= 1 then do; put "ERROR: the digital share must lie strictly between zero and one"; abort cancel; end;
  if missing(next_quarter_revenue_eur) or next_quarter_revenue_eur < 0 then do; put "ERROR: the outcome must be finite and non-negative"; abort cancel; end;
  log_recency = log(recency_days);
  frequency = frequency_12m;
  log_basket = log(avg_basket_eur);
  breadth = category_breadth;
  digital = digital_share;
  order = _n_;
  if eof then call symputx("customers", _n_);
run;
%if &customers < 100 %then %do; %put ERROR: too few customers for the declared segmentation; %abort cancel; %end;
/* Declared standardization: mean zero, unit standard deviation over the whole file. */
proc stdize data=p031 out=scaled method=std;
  var log_recency frequency log_basket breadth digital;
run;
/* Declared start: the observations at four fixed ranks of the composite score. */
data scored;
  set scaled;
  composite = log_recency + frequency + log_basket + breadth + digital;
run;
proc sort data=scored out=ranked; by composite order; run;
data seeds;
  set ranked;
  position = _n_;
  if position in (%sysfunc(round(0.125 * (&customers - 1))) + 1,
                  %sysfunc(round(0.375 * (&customers - 1))) + 1,
                  %sysfunc(round(0.625 * (&customers - 1))) + 1,
                  %sysfunc(round(0.875 * (&customers - 1))) + 1);
  keep log_recency frequency log_basket breadth digital;
run;
proc fastclus data=scaled seed=seeds maxclusters=4 maxiter=100 converge=0 out=clustered mean=centroids noprint;
  var log_recency frequency log_basket breadth digital;
run;
proc means data=clustered noprint nway;
  class cluster; var next_quarter_revenue_eur;
  output out=segment_outcome(drop=_type_) n=size mean=outcome_mean;
run;
proc sql noprint;
  select min(size) / &customers into :smallest_share from segment_outcome;
  select max(outcome_mean) / min(outcome_mean) into :ratio from segment_outcome;
  select mean(next_quarter_revenue_eur) into :grand from p031;
  select sum(size * (outcome_mean - &grand) ** 2) into :between from segment_outcome;
  select sum((next_quarter_revenue_eur - &grand) ** 2) into :total from p031;
quit;
/* Separation and stability are reported by the reference implementations; this
   program checks balance and usefulness, and prints the segment profile. */
data summary;
  customers = &customers;
  smallest_share = &smallest_share;
  balance_flag = ifc(smallest_share >= 0.05, "PASS", "FAIL");
  outcome_ratio_high_low = &ratio;
  outcome_variance_explained = &between / &total;
  usefulness_flag = ifc(outcome_ratio_high_low >= 1.50, "PASS", "FAIL");
  verdict = ifc(balance_flag = "PASS" and usefulness_flag = "PASS",
                "SEGMENTATION_READABLE_FOR_ACTION", "DIAGNOSTIC_BLOCKS_SEGMENTATION_READING");
run;
proc print data=summary noobs; run;
proc print data=centroids noobs; run;