#Developed 12/21/2025
"print_psych.cfa" <-
function(x,digits=2,all=FALSE,cut=NULL,sort=FALSE,suppress.warnings=TRUE,...)  {
if(!is.matrix(x) && !is.null(x$fa) && is.list(x$fa)) x <-x$fa   #handles the output from fa.poly
if(!is.matrix(x))  {if (!is.null(x$fn) ) {if(x$fn == "principal") {cat("Principal Components Analysis") } else {
 cat("Factor Analysis using method = ",x$fm )}  }} else {load <- x
                                                         class(x)<- NULL
                                                         print(round(x,digits))
                                                         
                                                         return()} 
   cat("\nCall: ")
   print(x$Call)
     
   	rload <- load <- x$loadings  
   	 if(is.null(cut)) cut <- 0   #caving into recommendations to print all loadings
 
   	nitems <- dim(load)[1]
 	nfactors <- dim(load)[2]
   	ncol <- dim(load)[2]
   	   cat('Confirmatory factor (Structure)  loadings\n')
    
     
     rloads <- round(load,digits)
	    	fx <- format(rloads,digits=digits)     	
	print(rloads)
	if(nfactors >1 ) {
	     fx <- format(rloads,digits=digits)
	    	nc <- nchar(rloads[1,3], type = "c")
	    	 fx.1 <- fx[,1,drop=FALSE]    #drop = FALSE  preserves the rownames for single factors
	    	 #fx.2 <- fx[,3:(2+ncol),drop=FALSE]
	    	 load.2 <- as.matrix(load)
	    	 
	    	 
	    	         	#fx[abs(load[,]) < cut,] <- paste(rep(" ", nc), collapse = "")
	} else {load.2<- as.matrix(load)}
	
	
	if(nfactors > 1) {if(is.null(x$Phi)) {h2 <- rowSums(load.2^2)} else {h2 <- diag(load.2 %*% x$Phi %*% t(load.2)) }} else {h2 <-load.2^2}
         #	if(!is.null(x$uniquenesses)) {u2 <- x$uniquenesses[u2.order]}  else 
         	{u2 <- (1 - h2)}
vtotal <- sum(h2 + u2)	
if(nfactors > 1)  {vx <- colSums(load^2) } else {vx <- sum(load^2)}
 names(vx) <- colnames(load)
          varex <- rbind("SS loadings" =   vx)
          varex <- rbind(varex, "Proportion Var" =  vx/vtotal)
         if (nfactors > 1) {
                      varex <- rbind(varex, "Cumulative Var"=  cumsum(vx/vtotal))
                       varex <- rbind(varex, "Proportion Explained"=  vx/sum(vx))
                      varex <- rbind(varex, "Cumulative Proportion"=  cumsum(vx/sum(vx))) 
                      }     	                                            	    	
	    	print(round(varex, digits))
	    	
	    if(!is.null(x$Phi))	{
	    	 cat('With correlations of \n')
	    	lowerMat(x$Phi)
	    	}
	    	
	    	#now print the results fron stats
	    	xx <- x #to save the original call
	    	x <- x$stats  #allowing calls using the calls from print.psych.fa
	if(!is.null(x$dof))  {cat("\nThe degrees of freedom of the model are ",x$dof)}
 if(!is.null(x$rms)) {cat("\nThe root mean square of the residuals (RMSR) is ", round(x$rms,digits),"\n") }
    if(!is.null(x$crms)) {cat("The df corrected root mean square of the residuals is ", round(x$crms,digits),"\n",...) }
    
     if((!is.null(x$nh)) && (!is.na(x$nh))) {cat("\nThe harmonic n.obs is " ,round(x$nh)) }
     if((!is.null(x$chi)) && (!is.na(x$chi))) {cat(" with the empirical chi square ", round(x$chi,digits), " with prob < ", signif(x$EPVAL,digits),"\n" ,...)  }
   	  
   	 if(!is.na(x$n.obs)) {cat("The total n.obs was ",x$n.obs, " with Likelihood Chi Square = ",round(x$STATISTIC,digits), " with prob < ", signif(x$PVAL,digits),"\n",...)}
  
     
   	if(!is.null(x$TLI)) {cat("\nTucker Lewis Index of factoring reliability = ",round(x$TLI,digits+1))}
   	if(!is.null(x$RMSEA)) {cat("\nRMSEA index = ",round(x$RMSEA[1],digits+1), " and the", (x$RMSEA[4])*100,"% confidence intervals are ",round(x$RMSEA[2:3],digits+1),...)  }
   	if(!is.null(x$BIC)) {cat("\nBIC = ",round(x$BIC,digits))}	    	

#get information from factor stats
if(!is.null(x$R2)) { stats.df <- t(data.frame(sqrt(x$R2),x$R2,2*x$R2 -1))

 rownames(stats.df) <- c("Correlation of (regression) scores with factors  ","Multiple R square of scores with factors ","Minimum correlation of possible factor scores ")
         colnames(stats.df) <- colnames(xx$loadings)
 } else {stats.df <- NULL}

 if(!is.null(stats.df)) { cat("\nMeasures of factor score adequacy             \n")
      print(round(stats.df,digits))}
}