/* History:
   2025-11-21 Moved IML modules for version 14.3 and 15.1 into their own catalog.
   2025-11-25 Move update_logger module to v143_modules.
*/

/* IML modules for version 14.3 */
proc iml;

	**Updates logger from IML;
	start update_logger(error_code,
											section,
											message);
											
		**Update logger dataset;
		edit _log_record var {section message error_code};
			append;
		close _log_record;
		
		**Write to log;
		if error_code = 1 then do;
			
			error_class = "NOTE";
		end;
		else if error_code = 2 then do;
			
			error_class = "WARNING";
		end;
		else do;
			
			error_class = "ERROR";
		end;
			
		print_command = cat('%put ', error_class, ': ', message, ';');
		call execute(print_command);
	finish update_logger;

	/* Check the existence of a file before the call to "use" */
	start checkFile(f, section);
	
		ret = exist(f, "data");
		if ret = 0 then do;
			msg = cat("Data ", f, " does not exist");
			run update_logger(3, section, msg);
		end;
		return ret;
	
	finish checkFile;

	/* The unique function does not seem to always work correctly for floating point values.
        Below, x is a column vector. The code preserves the order of vector x.	*/
	start myunique(x, tol);
	
		ret = x[1, 1];
		n   = max(nrow(x), ncol(x));
		if (n > 1) then do;
			do i = 2 to n;
				tmp = abs(ret - x[i, 1]) >= tol;
				if all(tmp) then ret = ret // x[i];
			end;
		end;
	  
		return ret;
	
	finish myunique;
	
	/* Module to return the eigenvalues and eigenvectors from alternative v14.3 code */
	start eigenFromSVD(mat, tol, onlyEigen);
	
	    ** Compute the eigenvalues;
		ee = eigval(mat);

        ** Check for complex eigenvalues;
		if (ncol(ee) > 1) then do;
			run update_logger(3, "eigenFromSVD module", "Complex eigenvalues computed");
			return .;
		end;

        if onlyEigen then do;
			ret = [#"eigenvalues" = ee, #"eigenvectors" = 0];
			return ret;
		end;

		/* Check for eigenvalues of 0 */
		if any(ee = 0) then do;
			run update_logger(3, "eigenFromSVD module", "Some eigenvalues = 0");
			return .;
		end;

		** Get the unique eigenvalues (unique function not working ???);
		uee  = myunique(ee, tol);
		
		n    = nrow(mat);
		neig = max(nrow(uee), ncol(uee));
		emat = j(n, n, 0);  /* return matrix of eigenvectors */
		c2   = 0;

		** Loop over unique eigenvalues to get the eigenvectors;
		do i = 1 to neig;
			lambda = uee[i];  
			B      = mat - lambda # I(n);
			call svd(U, Q, V, B); 
			tmp 		  = Q < tol;			
			if sum(tmp) = 0 then do;
				run update_logger(3, "eigenFromSVD module", "No eigenvectors for eigenvalue");
				return .;
			end;
			evecs         = V[, loc(tmp)];  
			nc            = ncol(evecs);
			c1            = c2 + 1;
			c2            = c1 + nc - 1; 
			if (c2 > n) then do;
				run update_logger(3, "eigenFromSVD module", "Too many eigenvectors computed");
				return .;
			end;
			emat[, c1:c2] = evecs;
		end;
		if (c2 < n) then do;
			run update_logger(3, "eigenFromSVD module", "Too few eigenvectors computed");
			return .;
		end;
		
		** Test;
		/*
		do i = 1 to n;
		  lambda = ee[i];
		  evec   = emat[, i];
		  v1     = mat*evec;
		  v2     = lambda # evec;
		  tmp    = max(abs(v1 - v2));
		  print tmp;
		end;
		*/
		
	    ret = [#"eigenvalues" = ee, #"eigenvectors" = emat];
	
	    return ret;
	
	finish eigenFromSVD;
	
	start eigenDecomp(mat, alg, onlyEigen=);

	    ** alg = 1: call eigen function, alg = 2 use code below;
		
	    tol = 1e-8; /* tolerance */
	
	    if (alg = 1) then do;
			call eigen(eigvals, eigvecs, mat);
			ret = [#"eigenvalues" = eigvals, #"eigenvectors" = eigvecs];
		end;
		else do;	
			if isskipped(onlyEigen) then onlyEigen = 0;
			ret = eigenFromSVD(mat, tol, onlyEigen);
		end;
	
		return ret;
	
	finish eigenDecomp;
	
	/* Functions on matrices */
	start matIsSquare(A);
	
		ret = 0;
		if (nrow(A) = ncol(A)) then ret = 1;
	
		return ret;
	  
	finish matIsSquare;
	
	start matIsDiagonal(A);
	
		ret = 0;
		if matIsSquare(A) then do;
			diagA = diag(A);
			if all(A = diagA) then ret = 1;
		end;
		return ret;
	
	finish matIsDiagonal;
	
	start matIsSymmetric(A);
	
		ret = 0;
		if matIsSquare(A) & all(abs(A - A`) < 10**-8) then ret = 1;
		
		return ret;
		
	finish matIsSymmetric;
	
	start matIsColVec(A);
	
		if (ncol(A) = 1) then ret = 1;
		else ret = 0;
		
		return ret;
	
	finish matIsColVec;
	
	start matIsRowVec(A);
	
		if (nrow(A) = 1) then ret = 1;
		else ret = 0;
		
		return ret;
	
	finish matIsRowVec;

	/* Matrix multiplication functions */
	start matMultADAt(A, vec, alg);

		/* A is n by n matrix, vec is a vector of length n.
	    compute the product A * diag(vec) * A`. */
	
		if (alg = 1) then do;		
			ret = A * diag(vec) * A`;
		end;
		else do;
			n	= max(nrow(A), ncol(A));
			AD	= j(n, n, 0); /* To store A * D */
			ret = j(n, n, 0); /* return matrix */
			do i = 1 to n;
				AD[, i] = A[, i] # vec[i];
			end;
			do i = 1 to n;
				veci = AD[i, ];
				do j = i to n;
					vecj      = A[j, ];
					val       = sum(veci # vecj);
					ret[i, j] = val;
					ret[j, i] = val;
				end;
			end;
		end;
	
		return ret;
	
	finish matMultADAt;

	start matMultAB_general(A, B);
	
		nr  = nrow(A);
		nc  = ncol(B);
		ret = j(nr, nc, 0);
		Bt  = B`;
		do i = 1 to nr;
			veci = A[i, ];
			do j = 1 to nc;
				vecj      = Bt[j, ];
				ret[i, j] = sum(veci # vecj);
			end;
		end;
		
		return ret;
	
	finish matMultAB_general;
	
	/* A * B */
	start matMultAB(A, B, alg);

		if (alg = 1) then do;
			ret = A*B;
		end;
		else do;
			*Special case: matrix is zero;
			if all(A = 0) | all(B = 0) then do;
				ret = j(nrow(A), ncol(B), 0);			  
			end;
			else if matIsRowVec(A) & matIsColVec(B) then do;
				*Special case: dot product;
				Bt  = B`;
				ret = sum(A # Bt);
			end;
			else if matIsDiagonal(A) & matIsColVec(B) then do;
				*Special case: Diagonal matrix times a column vector;	
				ret = vecdiag(A) # B;
			end;
			else do;
				ret = matMultAB_general(A, B);
			end;
		end;
	
		return ret;
	
	finish matMultAB;
	
	/* A * A` */
	start matMultAAt(A, alg);

		if (alg = 1) then do;
			ret = A*A`;
		end;
		else do;
			* Special case when A is diagonal;
			if matIsDiagonal(A) then do;
				dvec = vecdiag(A);
				ret  = diag(dvec # dvec);
			end;
			else do;
				ret = matMultAB(A, A`, alg);
			end;
		end;
		
		return ret;
	
	finish matMultAAt;
	
	/* A * B * A` */
	start matMultABAt(A, B, alg);

		if (alg = 1) then do;
			ret = A*B*A`;
		end;
		else do;
			nr  = nrow(A);
			ret = j(nr, nr, 0);
			AB  = matMultAB(A, B, alg);
			do i = 1 to nr;
				veci = AB[i, ];
				do j = 1 to nr;
					vecj      = A[j, ];
					ret[i, j] = sum(veci # vecj);
				end;
			end;
		end;
	
		return ret;
	
	finish matMultABAt;
	
	/* A * B * C */
	start matMultABC(A, B, C, alg);
	
		if (alg = 1) then do;
			ret = A*B*C;
		end;
		else do;
			AB  = matMultAB(A, B, alg);
			ret = matMultAB(AB, C, alg);
		end;
		
		return ret;
	
	finish matMultABC;

	/* Function to compute random matrix from Wishart distrubution */
	start myRandWishart(df, S, alg);

		nr = nrow(S); 
		if df < nr then do;
			/* Stop with error */
			run update_logger(3, "myRandWishart module", "degrees of freedom < dimension of matrix");
			return .;
		end;

		if (alg = 1) then do;
			ret = randwishart(1, df, S);
			ret = shape(ret, nr, nr);
		end;
		else do;

			*Cholesky decomposition on the scale matrix S;
			L   = root(S); /* Upper triangle */
			
			*Generate df mvn vectors;
            m   = floor(df);
			Z   = j(m, nr, 0);
			mu  = j(1, nr, 0);
			sig = i(nr);
			do i = 1 to m;
				Z[i, ] = randnormal(1, mu, sig);
			end;
			free mu sig;
			tmp = matMultAB(Z, L, alg);
			ret = matMultAAt(tmp`, alg);	
		end;
		
		return ret;

	finish myRandWishart;


	reset storage=WORK.v143_modules;
	store module=_all_;
quit;
