*! version 1.0.0  20aug2025  I I Bolotov
program define clsp, eclass byable(recall)
	version 16.0
	/*
		Convex Least Squares Programming (CLSP) is a modular two-step estimator 
		for solving underdetermined, ill-posed, or structurally constrained     
		least-squares problems. It combines pseudoinverse-based estimation with 
		convex-programming correction methods inspired by Lasso, Ridge, and     
		Elastic Net to ensure numerical stability, constraint enforcement, and  
		interpretability. The package also provides numerical stability analysis
		and CLSP-specific diagnostics, including partial R^2, normalized RMSE   
		(NRMSE), Monte Carlo t-tests for mean NRMSE, and condition-number-based 
		confidence bands.                                                       

		Author: Ilya Bolotov, MBA, Ph.D.                                        
		Date: 20 August 2025                                                    
	*/
	tempname vars CLSP warnings error
	// replay last result                                                       
	if replay() {
		if _by() {
			error 190
		}
		cap conf mat e(z)
		if _rc {
			di as err "results of clsp not found"
			exit 301
		}
		_summary									   // print summary
		exit 0
	}
	// syntax                                                                   
	syntax																	///
	[anything] [if] [in], [													///
		b(string) Constraints(string) Slackvars(string) MODel(string)		///
		m(int 0) p(int 0) i(int 1) j(int 1) ZERODiagonal r(int 1) Z(name)	///
		RCOND(real  0) TOLerance(real -1) ITERationlimit(int 0) noFINAL		///
		alpha(numlist >=-1 <=1 sort) CVXopt(string asis) MISSingokay *		///
	]
	// adjust and preprocess options for Python's CLSP module                   
	if           ustrregexm( `"`anything'"',       "^estat",     1)           {
		_estat `=ustrregexrf(`"`anything'"',       "^estat", "", 1)', `options'
		exit 0
	}
	if ("`b'"                                                == "")           {
		di as err "option b() required"
		exit 198
	}
	if ("`b'"                                                != "")           {
			   cap unab varlist :     `b',         max(1)
		loc rc =   _rc
		if _rc cap conf mat           `b'
		else       loc  b             `varlist'
		if _rc     &  `rc' != 103 {
			di as err                 "b: `b'"
			di as err                 "invalid varname or matrix"
			exit cond(`rc' == 103, 103, 198)
		}
		else       if `rc' == 103 error 103
	}
	if ("`constraints'"                                      != "")           {
			   cap unab varlist :     `constraints'
		if _rc cap conf mat           `constraints'
		else       loc  constraints   `varlist'
		if _rc                                                                {
			di as err                 "constraints: `constraints'"
			di as err                 "invalid varlist or matrix"
			exit 198
		}
	}
	if ("`slackvars'"                                        != "")           {
			   cap unab varlist :     `slackvars'
		if _rc cap conf mat           `slackvars'
		else       loc  slackvars     `varlist'
		if _rc                                                                {
			di as err                 "slackvars: `slackvars'"
			di as err                 "invalid varlist or matrix"
			exit 198
		}
	}
	if ("`model'"                                            != "")           {
			   cap unab varlist :     `model'
		if _rc cap conf mat           `model'
		else       loc  model         `varlist'
		if _rc                                                                {
			di as err                 "model: `model'"
			di as err                 "invalid varlist or matrix"
			exit 198
		}
	}
	if (`m' <  0 | `p' <  0 | `i' <  0 | `j' <  0)                            {
		di as err "m(), p(), i(), and j() must be non-negative"
		exit 198
	}
	if (`r' <  0 | `iterationlimit'          <  0)                            {
		di as err "r() and iterationlimit() must be at least 1"
		exit 198
	}
	if ("`z'"                                                != "")         ///
	conf mat  `z'
	if ! (`rcond'     == -1 | (`rcond'     >= 0 & `rcond'     < 1))           {
		di as err "rcond() must be -1 or lie in [0,1)"
		exit 198
	}
	if ! (`tolerance' == -1 | (`tolerance' >= 0 & `tolerance' < 1))           {
		di as err "tolerance() must be -1 or lie in [0,1)"
		exit 198
	}
	if ("`alpha'"                                            != "")           {
		numlist " `alpha'"
		loc        alpha       `r(numlist)'
		foreach    a in        `alpha'                            {
			if ! (`a' == -1 | (`a' >= 0 & `a' <= 1)) {
				di as err "each value in alpha() must be -1 or lie in [0,1]"
				exit 198
			}
			if   (`a' == -1                        ) {
				loc alpha "None"
				continue, break
			}
		}
		if `: word count of `alpha''  > 1									///
			loc alpha =      "[" + ustrregexra("`alpha'", "\s+", ", ") + "]"
	}
	else    loc alpha     "None"
	loc   problem "'general'"
	mata: if (ustrregexm(`"`anything'"',                       "ap|tm", 1)) ///
		  st_local("problem",                                 "'ap'"      );;
	mata: if (ustrregexm(`"`anything'"',                     "cmls|rp", 1)) ///
		  st_local("problem",                               "'cmls'"      );;
	loc   _a       ", selectvar='`touse'', missingval=nan), dtype=float64"
	loc   C_array  "None"
	if  "`constraints'"  != ""                                                {
		mata: if (! hasmissing( _st_varindex(tokens("`constraints'"    )))) ///
			  st_local("C_array", "_(array(Data.get('`constraints''`_a'))");;
		mata: if (J(0,0,.)        !=      st_matrix("`constraints'"      )) ///
			  st_local("C_array", "array(Matrix.get('`constraints''    ))");;
	}
	loc   S_array  "None"
	if  "`slackvars'"    != ""                                                {
		mata: if (! hasmissing( _st_varindex(tokens("`slackvars'"      )))) ///
			  st_local("S_array", "_(array(Data.get('`slackvars''  `_a'))");;
		mata: if (J(0,0,.)        !=      st_matrix("`slackvars'"        )) ///
			  st_local("S_array", "array(Matrix.get('`slackvars''      ))");;
	}
	loc   M_array  "None"
	if  "`model'"        != ""                                                {
		mata: if (! hasmissing( _st_varindex(tokens("`model'"          )))) ///
			  st_local("M_array", "_(array(Data.get('`model''      `_a'))");;
		mata: if (J(0,0,.)        !=      st_matrix("`model'"            )) ///
			  st_local("M_array", "array(Matrix.get('`model''          ))");;
	}
	loc   b_array  "None"
	if  "`b'"            != ""                                                {
		mata: if (! hasmissing( _st_varindex(tokens("`b'"              )))) ///
			  st_local("b_array", "_(array(Data.get('`b''          `_a'))");;
		mata: if (J(0,0,.)        !=      st_matrix("`b'"                )) ///
			  st_local("b_array", "array(Matrix.get('`b''              ))");;
	}
	loc   m              =cond( `m',                "`m'",            "None")
	loc   p              =cond( `p',                "`p'",            "None")
	loc   i              =cond( `i',                 `i',                  1)
	loc   j              =cond( `j',                 `j',                  1)
	loc   zerodiagonal   =cond("`zerodiagonal'"!= "",        "True", "False")
	loc   r              =cond( `r',                 `r',                  1)
	loc   Z_array  "None"
	if  "`z'"            != ""                                                {
		mata: if (J(0,0,.)          !=    st_matrix("`z'"                )) ///
			  st_local("Z_array", "array(Matrix.get('`z''              ))");;
	}
	loc   rcond          =cond( `rcond'        >   0,     "`rcond'",		///
						  cond( `rcond'        <   0,        "True", "False"))
	loc   tolerance      =cond( `tolerance'    >=  0, "`tolerance'",  "None")
	loc   iterationlimit =cond( `iterationlimit',"`iterationlimit'",  "None")
	loc   final          =cond("`final'"       == "",        "True", "False")
	if `"`cvxopt'"'      != ""												///
	loc   cvxopt         =", " + ustrregexrf(`"`cvxopt'"',  "^\s*,",      "")
	// estimate with the help of Python's CLSP module                           
	sca                 `vars'    = ""
	foreach name in     `constraints' `slackvars' `model' `b' `z'             {
		cap conf var    `name'
		if ! _rc sca    `vars'    =   `vars' +  " `name'"
	}
	tempvar              touse
	mark                `touse'  `if' `in'
	if ("`missingokay'"  == "")												///
	markout             `touse'  `=   `vars''
	if ("`b_array'"  == "None")                                               {
		di as err "RHS error: b() is misspecified"
		error 7102
	}
	if ("`C_array'"  == "None"  &        "`M_array'"  == "None"  &			///
		"`m'"        == "None"  &        "`p'"        == "None")              {
		di as err "LHS error: all of constraints(), model(), m(), and p() " ///
							 "are misspecified"
		error 7102
	}
	python: `CLSP', `warnings', `error' = _(None, False,                    ///
	                                         fun=CLSP().solve,              ///
	                              **dict(problem=`problem',                 ///
	                                           C=`C_array',                 ///
	                                           S=`S_array',                 ///
	                                           M=`M_array',                 ///
	                                           b=`b_array',                 ///
	                                           m=`m',                       ///
	                                           p=`p',                       ///
	                                           i=`i',                       ///
	                                           j=`j',                       ///
	                               zero_diagonal=`zerodiagonal',            ///
	                                           r=`r',                       ///
	                                           Z=`Z_array',                 ///
	                                       rcond=`rcond',                   ///
	                                   tolerance=`tolerance',               ///
	                             iteration_limit=`iterationlimit',          ///
	                                       final=`final',                   ///
	                                       alpha=`alpha'    `cvxopt'      ))
	python:         `warnings' = [SFIToolkit.displayln('{txt}note: '+k    ) ///
	                              for k in dict.fromkeys(str(w.message    ) ///
	                              for w in `warnings'                     )]
	python:                       Macro.setLocal("`error'", `error'       )
	if "``error''"  != ""                                                     {
		di as err   "``error''"
		exit 7102
	}
	mata: st_eclear()                                  // store results
	eret post,         esample(`touse')
	qui  count   if    e(sample)
	eret sca N      =  r(N)
	foreach  s   in    iteration_limit r tolerance alpha kappaC kappaB		///
					   kappaA r2_partial nrmse nrmse_partial                  {
		loc  code             "`CLSP'.`s'"
		_store  `s',   py(`code'       ) r(e) t(scalar)
	}
	eret loc cmdline   clsp  `0'
	eret loc cmd       clsp
	eret loc estat_cmd clsp estat
	eret loc title     CLSP estimation
	foreach  s   in    final seed distribution                                {
		loc  code   =  cond(! inlist("`s'", "distribution"      ),			///
							  "`CLSP'.`s'",									///
							  "`CLSP'.`s'.__code__.co_names[-1]")
		_store  `s',   py(`code', True ) r(e) t(macro )
	}
	foreach  s   in    A C_idx b Z zhat z x y z_lower z_upper x_lower		///
											  x_upper y_lower y_upper         {
		loc  code             "`CLSP'.`s'"
		_store  `s',   py(`code'       ) r(e) t(matrix)
	}
	python:  del   `warnings', `error'
	mata: st_global("s(clsp)","`CLSP'" )
	_summary, noread								   // store summary
end

program define _estat, eclass
	tempname tmp warnings error bar
	// assert last result                                                       
	cap conf mat e(z)
	if _rc {
		di as err "results of clsp not found"
		exit 301
	}
	// syntax                                                                   
	syntax																	///
	[anything], [															///
		RESET THRESHold(real 0) BAR HBAR MATrix(name) NCOLs(int 1)			///
		SAMPLEsize(int 50) SEED(int 0) DISTribution(string) PARTial			///
		SIMulate Level(cilevel)												///
	]
	// adjust and preprocess options for Python's CLSP module                   
	if ("`anything'"                                         == "")           {
		di as err "subcommand required"
		exit 321
	}
	if (`threshold'                          <  0)                            {
		di as err "threshold() must be non-negative"
		exit 198
	}
	if ("`bar'"               != "" & "`hbar'"               != "")           {
		di as err "must specify either bar or hbar option"
		exit 198
	}
	loc   reset          =cond("`reset'"       != "",        "True", "False")
	loc   matrix         =cond("`matrix'"      != "",    "`matrix'","rmsa_i")
	loc   seed           =cond( `seed',                    "`seed'",  "None")
	loc   distribution   =cond("`distribution'"!= "",                       ///
												  "'`distribution''", "None")
	loc   partial        =cond("`partial'"     != "",        "True", "False")
	loc   simulate       =cond("`simulate'"    != "",        "True", "False")
	// compute the structural correlogram of the CLSP constraint system         
	if  ustrregexm( `"`anything'"',                "^corr",      1)           {
		_read
		python:  `tmp', `warnings', `error' = _(None, False,                ///
		                                         fun=`s(clsp)'.corr,        ///
		                                **dict(reset=`reset',               ///
		                                   threshold=`threshold'          ))
		python:         `warnings' = [SFIToolkit.displayln('{txt}note: '+k) ///
		                              for k in dict.fromkeys(str(w.message) ///
		                              for w in `warnings'                 )]
		python:                       Macro.setLocal("`error'", `error'   )
		if "``error''"  != ""												  {
			di as err      "``error''"
			exit 7102
		}
		loc CLSP            `s(clsp)'				   // store results
		foreach  s   in    rmsa                                               {
			loc  code      "  `CLSP'.`s'"
			_store  `s',   py(`code'       ) r(e) t(scalar)
		}
		foreach  s   in    rmsa_i rmsa_dkappaC rmsa_dkappaB rmsa_dkappaA	///
						   rmsa_dnrmse rmsa_dzhat rmsa_dz rmsa_dx             {
			loc  code      "  `CLSP'.`s'"
			_store  `s',   py(`code'       ) r(e) t(matrix)
		}
		python:  del   `warnings', `error'
		mata: st_global("s(clsp)","`CLSP'" )
		_summary, noread							   // store summary
		if "`bar'"      != ""  |											///
		   "`hbar'"     != ""                                                 {
			loc b        = cond("`bar'" != "",    "bar", "hbar")
			tempname f
			mkf     `f'
			frame   `f':   qui svmat    e(`matrix'), n(`matrix')
			frame   `f':   qui g constraint    = _n
			frame   `f':   if `: colsof e(`matrix')' < 2 ren `matrix'* `matrix'
			loc j    = 1							   // plot result
			frame   `f':   foreach var of varl               `matrix'*        {
				graph `b' `var', over(constraint) title("`var'")			///
								 name(_clsp`j', replace) nodraw
				loc  g "`g' _clsp`j'"
				loc  ++j
			}
			graph combine   `g', cols(`ncols'   )							///
								 name(`matrix', replace) nodraw
			graph drop      `g'
			graph display   `matrix'
		}
		exit 0
	}
	// perform bootstrap or Monte Carlo t-tests on the NRMSE statistic from the 
	// CLSP estimator                                                           
	if  ustrregexm( `"`anything'"',                "^ttest",     1)           {
		_read
		python:  `tmp', `warnings', `error' = _(None, False,                ///
		                                         fun=`s(clsp)'.ttest,       ///
		                                **dict(reset=`reset',               ///
		                                 sample_size=`samplesize',          ///
		                                        seed=`seed',                ///
		                                distribution=`distribution',        ///
		                                     partial=`partial',             ///
		                                    simulate=`simulate'           ))
		python:         `warnings' = [SFIToolkit.displayln('{txt}note: '+k) ///
		                              for k in dict.fromkeys(str(w.message) ///
		                              for w in `warnings'                 )]
		python:                       Macro.setLocal("`error'", `error'   )
		if "``error''"  != ""												  {
			di as err      "``error''"
			exit 7102
		}
		loc CLSP            `s(clsp)'				   // store results
		if "`seed'"                  != "None"								///
		mata: st_global("e(seed)",         "`seed'"                           )
		if "`distribution'"          != "None"								///
		mata: st_global("e(distribution)", subinstr("`distribution'","'","",.))
		foreach  s   in    nrmse_ttest                                        {
			loc  code      "  `CLSP'.`s'"
			_store  `s',   py(`code'       ) r(e) t(matrix)
		}
		tempname f
		mkf     `f'
		frame   `f':   qui svmat    e(nrmse_ttest), n(nrmse_ttest)
		frame   `f':   ren            nrmse_ttest*    nrmse_ttest
		frame   `f':   qui sum
		loc t   `r(N)' `r(mean)' `r(sd)' `e(nrmse)'
		python:  del   `warnings', `error'
		mata: st_global("s(clsp)","`CLSP'" )
		_summary, noread							   // store summary
		ttesti  `t', l(`level')						   // print ttest
		exit 0
	}
	di as err "{bf:estat `anything'} not valid"
	exit 321
end

program define _read, sclass
	version 16.0
	tempname tmp
	// scalars, macros, and matrices in e()                                     
	cap python:     `tmp'      = CLSP()
	foreach  s   in    iteration_limit r tolerance alpha kappaC kappaB		///
					   kappaA rmsa r2_partial nrmse nrmse_partial             {
		cap python: `tmp'.`s'  =      Scalar.getValue('e(`s')')
	}
	foreach  s   in    final seed distribution                                {
			 if      "`s'"     == "final"									///
		cap python: `tmp'.`s'  = bool(Macro.getGlobal('e(`s')'))
		else if      "`s'"     == "seed"									///
		cap python: `tmp'.`s'  =  int(Macro.getGlobal('e(`s')'))
		else if      "`s'"     == "distribution"                              {
			python: `tmp'.rng  = random.default_rng(`tmp'.seed)
			python: `tmp'.`s'  = lambda n:`tmp'.rng.`e(`s')'(size=(n, 1))
		}
		else																///
		cap python: `tmp'.`s'  =      Macro.getGlobal('e(`s')')
	}
	foreach  s   in    A C_idx b Z zhat z x y rmsa_i rmsa_dkappaC			///
					   rmsa_dkappaB rmsa_dkappaA rmsa_dnrmse rmsa_dzhat		///
					   rmsa_dz rmsa_dx nrmse_ttest z_lower z_upper x_lower	///
					   x_upper y_lower y_upper                                {
			 if      "`s'"     == "C_idx"									///
		cap python: `tmp'.`s'  =    [int(i) for i in Matrix.get('e(C_idx)')[0]]
		else if      "`s'"     == "b"										///
		cap python: `tmp'.`s'  =     array(Matrix.get('e(rhs)'), dtype=float64)
		else																///
		cap python: `tmp'.`s'  =     array(Matrix.get('e(`s')'), dtype=float64)
	}
	cap sret loc clsp   `tmp'						   // store object
end

program define _reorder, eclass
	version 16.0
	loc rows =  e(C_idx)[1, 1]
	loc cols =  e(C_idx)[1, 2]
	// scalars, macros, and matrices in e()                                     
	foreach  s   in    N iteration_limit r tolerance alpha kappaC kappaB	///
					   kappaA rmsa r2_partial nrmse nrmse_partial             {
		mata:          st_numscalar("e(`s')",      st_numscalar("e(`s')"),  ///
									anyof(("iteration_limit"    ),"`s'" )   ///
									?      "hidden" :"visible"          )
	}
	foreach  s   in    distribution seed final title estat_cmd cmd cmdline    {
		mata:             st_global("e(`s')",         st_global("e(`s')"),  ///
									anyof(("N/A"                ),"`s'" )   ///
									?      "hidden" :"visible"          )
	}
	foreach  s   in    y_upper y_lower x_upper x_lower z_upper z_lower		///
					   nrmse_ttest rmsa_dx rmsa_dz rmsa_dzhat rmsa_dnrmse	///
					   rmsa_dkappaA rmsa_dkappaB rmsa_dkappaC rmsa_i y x z	///
					   zhat Z rhs C_idx A                                     {
		mata:             st_matrix("e(`s')",         st_matrix("e(`s')"),  ///
									anyof(("Z","rhs","C_idx","A"),"`s'" ) | ///
									ustrregexm("`s'","(rmsa|ttest)",   1)   ///
									?      "hidden" :"visible"          )
		cap conf mat e(`s')
		if      _rc                   continue
		if      inlist(  "`s'",            "x","x_lower","x_upper"      )	///
		mata:   st_local( "x" ,     cols(st_matrix("e(`s')")) ==  1         ///
			                               ?      `","`s'""'      :  "" )
		if               "`x'" == ""  continue   /* skip if non-vector       */
		mata:   tmp =           st_matrixrowstripe("e(`s')"             )
			 if inlist(  "`s'", "A","rhs"                               )	///
		mata:   tmp[.,          1] =             J( rows(tmp),    1, "M");  ///
			    tmp[1..`rows',  1] =             J(`rows',        1, "C");
		else if inlist(  "`s'", "C_idx"                                 )	///
		mata:   tmp =                             ("",               "C");
		else if inlist(  "`s'", "Z","zhat","z","z_lower","z_upper"   `x')	///
		mata:   tmp[.,          1] =             J( rows(tmp),    1, "S");  ///
			    tmp[1..`cols',  1] =             J(`cols',        1, "C");
		else if inlist(  "`s'", "y","y_lower","y_upper"                 )	///
		mata:   tmp[.,          1] =             J( rows(tmp),    1, "S");
		mata:                   st_matrixrowstripe("e(`s')",                ///
			                    /* row    stripes     */ tmp            )
		mata:   tmp =           st_matrixcolstripe("e(`s')"             )
			 if inlist(  "`s'", "A","Z"                                 )	///
		mata:   tmp[.,          1] =             J( rows(tmp),    1, "S");  ///
			    tmp[1..`cols',  1] =             J(`cols',        1, "C");
		else if inlist(  "`s'", "C_idx"                                 )	///
		mata:   tmp =                             ("","rows")\("","cols");
		else if inlist(  "`s'", "rmsa_dzhat","rmsa_dz","rmsa_dx"        )	///
		mata:   tmp[.,          1] =             J( rows(tmp),    1, "S");  ///
			    tmp[1..`cols',  1] =             J(`cols',        1, "C");
		else																///
		mata:   tmp =                             ("",             "`s'");
		mata:                   st_matrixcolstripe("e(`s')",                ///
			                    /* column stripes     */ tmp            )
	}
	// macros and matrices in r()                                               
	foreach  s   in    inverse                                                {
		mata:             st_global("r(`s')",         st_global("r(`s')"))
	}
	foreach  s   in    y_upper y_lower x_upper x_lower z_upper z_lower		///
					   rmsa_dx rmsa_dz rmsa_dzhat rmsa_dnrmse rmsa_dkappaA	///
					   rmsa_dkappaB rmsa_dkappaC rmsa_i                       {
		mata:             st_matrix("r(`s')",         st_matrix("r(`s')"))
		cap conf mat r(`s')
		if _rc   continue
		mata:                   st_matrixrowstripe("r(`s')",  ("","`s'" ))
		mata:                   st_matrixcolstripe("r(`s')",                ///
			                                      ("","min" )\("","max" )\  ///
			                                      ("","mean")\("","sd"  ))
	}
end

program define _store, eclass
	version 16.0
	tempname tmp
	// syntax
	syntax name, PYthon(string asis) Result(string) Type(string)
	loc result =             substr(strlower("`result'"),             1, 1)
	loc type   =          strproper(strlower("`type'"                    ))
	// scalars, macros, and matrices in e() and r()                             
	mata:     st_local("_type",    ("scalar",      "local",   "matrix"   )[ ///
		               selectindex(("Scalar",      "Macro",   "Matrix"   )  ///
		                           :==                          "`type'")])
	mata:     st_local("_attr",    ("setValue",    "setLocal", "store"   )[ ///
		               selectindex(("Scalar",      "Macro",   "Matrix"   )  ///
		                           :== "`type'"                         )])
	mata:     st_local("_fun",     ("numscalar",   "global",  "matrix"   )[ ///
		               selectindex(("Scalar",      "Macro",   "Matrix"   )  ///
		                           :==                          "`type'")])
	cap python:      `type'.`_attr'('`tmp'',                  _( `python'))
	if _rc																	///
	mata:    st_`_fun'(   "`result'(`namelist')",  .                      )
	if       "`namelist'" == "b"                   loc namelist  "rhs"
	if       "`namelist'" == "C_idx"  |										///
			 "`result'"   == "r"											///
	mata:    st_matrix(   "`tmp'",                    st_matrix("`tmp'" )')
	mata:    st_`_fun'(   "`result'(`namelist')",                           ///
		     st_`=cond(   "`_fun'"!="global","`_fun'","local")'("`tmp'" ) )
end

program define _summary, sclass
	version 16.0
	tempname tmp warnings error
	if "`1'" == ""    _read
	python:  `tmp', `warnings', `error' = _(None, False,                    ///
	                                         fun=`s(clsp)'.summary,         ///
	                              **dict(display=False                    ))
	python:         `warnings' = [SFIToolkit.displayln('{txt}note: '+k)     ///
		                          for k in dict.fromkeys(str(w.message)     ///
		                          for w in `warnings'                     )]
	python:                       Macro.setLocal("`error'", `error'       )
	if "``error''"  != ""                                                     {
		di as err      "``error''"
		exit 7102
	}
	mata:    st_rclear()                               // store results
	foreach  s   in    inverse                                                {
		loc  code      "`tmp'['`s'']"
		cap _store `s', py(`code', True ) r(r) t(macro )
	}
	foreach  s   in    rmsa_i rmsa_dkappaC rmsa_dkappaB rmsa_dkappaA		///
					   rmsa_dnrmse rmsa_dzhat rmsa_dz rmsa_dx				///
					   z_lower z_upper x_lower x_upper y_lower y_upper        {
		loc  code      "`tmp'['`s'']"
		cap _store `s', py(`code'       ) r(r) t(matrix)
	}
	foreach  v   in           `s(clsp)' `tmp' `warnings' `error'              {
		cap  python:    del   `v'
	}
	cap sret loc clsp  ""
	_reorder
	loc p  "partial"
	di as  res  _n "Estimator Configuration:"		   // print results
	mata:  display(_sprint("Generalized inverse",  st_global("r(inverse)"  )))
	mata:  display(_sprint("Iterations (r)",    st_numscalar("e(r)"        )))
	mata:  display(_sprint("Tolerance",         st_numscalar("e(tolerance)")))
	mata:  display(_sprint("Final correction",     st_global("e(final)"    )))
	mata:  display(                                usubinstr(               ///
				  _sprint("Regularization (a)", st_numscalar("e(alpha)")),  ///
				  "(a)",  "(α)",  .                                         ))
	di as  res  _n "Numerical Stability:"
	mata:  display(_sprint("kappaC",            st_numscalar("e(kappaC)"   )))
	mata:  display(_sprint("kappaB",            st_numscalar("e(kappaB)"   )))
	mata:  display(_sprint("kappaA",            st_numscalar("e(kappaA)"   )))
	cap conf sca e(rmsa)
	if ! _rc																///
	mata:  display(_sprint("rmsa",              st_numscalar("e(rmsa)"     )))
	cap conf mat r(rmsa_dz)
	if ! _rc																///
	foreach s in   rmsa_i rmsa_dkappaC rmsa_dkappaB rmsa_dkappaA			///
				   rmsa_dnrmse rmsa_dzhat rmsa_dz rmsa_dx                     {
		mata: display(_sprint("`s'",               st_matrix("r(`s')"      )))
	}
	di as  res  _n "Goodness of Fit:"
	mata:  display(_sprint("r2_`p'",            st_numscalar("e(r2_`p')"   )))
	mata:  display(_sprint("nrmse",             st_numscalar("e(nrmse)"    )))
	mata:  display(_sprint("nrmse_`p'",         st_numscalar("e(nrmse_`p')")))
	cap conf mat r(z_lower)
	if ! _rc																///
	foreach s in   z_lower z_upper x_lower x_upper y_lower y_upper            {
		mata: display(_sprint("`s'",               st_matrix("r(`s')"      )))
	}
end

* Mata code ***************                                                     
version 16.0

loc RS           real scalar
loc SS           string scalar
loc TR           transmorphic rowvector

mata:
mata set matastrict on

`SS' _sprint(`SS' l,    `TR' v)
{
    if (orgtype(v) == "scalar")
    return("  {txt}" + sprintf("%-20s", l + ":")                              +
             "{res}" + sprintf( "%16s", _format(v,    10)))
    else
    return("  {txt}" + sprintf("%-16s", l + ":") +  "min="                    +
             "{res}" + sprintf( "%10s", _format(v[1], 10))                    +
           "  {txt}"                             +  "max="                    +
             "{res}" + sprintf( "%10s", _format(v[2], 10))                    +
           "  {txt}"                             + "mean="                    +
             "{res}" + sprintf( "%10s", _format(v[3], 10))                    +
           "  {txt}"                             +   "sd="                    +
             "{res}" + sprintf( "%10s", _format(v[4], 10)))
}
`SS' _format(`TR' x, |  `RS' w)
{
    `RS'  dp, ep
    `SS'  fixed
    w     = w != . ? w : 10
    dp    = max((3, trunc(w / 2) - 1))
    ep    = max((2, dp           - 1))

    if (isstring(x)) return(sprintf("%"+strofreal(w)+"s",                   x))
    fixed =                 sprintf("%"+strofreal(w)+"."+strofreal(dp)+"f", x)
    if ( missing(x)) return(sprintf("%"+strofreal(w)+"."+strofreal(dp)+"f", .))
    if ( !  (abs(x)  >= 10^(-(dp + 1))      & abs(x)    <  10^(w - dp)      &
          strlen(fixed)    <=      w))
                     return(sprintf("%"+strofreal(w)+"."+strofreal(ep)+"e", x))
    else             return(fixed)
}
end

* Python 3 code ***********                                                     
version 16.0
python:
# Stata Function Interface
from sfi import  Data, Macro, Matrix, Scalar, SFIToolkit
# Python Modules
from   warnings  import catch_warnings, simplefilter
try:
    from   numpy import nan, array, ndarray, float64, random
    from   clsp  import CLSP
except:
    SFIToolkit.errprint('pyclsp not found; install from PyPI')
    SFIToolkit.error(7102)

def _(x, flag=False, fun=None, **kwargs):              # preprocess output
    if  fun is None:
        return  (str(x)           if flag                         else
                 float('nan')     if x is None  or (isinstance(x, ndarray)
                                                and x.size  == 0) else
                 x.reshape(-1, 1) if               (isinstance(x, ndarray)
                                                and x.ndim  == 1) else
                 array([a.ravel() for  a in x],             dtype=float64)
                                  if                isinstance(x, list   )
                                  and           all(isinstance(a, ndarray)
                                  for  a in x)                    else
                 list(x.values()) if                isinstance(x, dict   )
                                                                  else  x)
    else:
        result, warnings, error = None, [], ''
        with catch_warnings(record=True) as warnings:
            simplefilter('always')
            try:
                result = fun(**kwargs)
            except Exception as e:
                error  = str(   e)
        return  result, warnings, error
end