#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
# 
#   http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

#  
# THIS SCRIPT CHOOSES A GLM REGRESSION MODEL IN A STEPWISE ALGIRITHM USING AIC
# EACH GLM REGRESSION IS SOLVED USING NEWTON/FISHER SCORING WITH TRUST REGIONS
#
# INPUT PARAMETERS:
# ---------------------------------------------------------------------------------------------
# NAME  TYPE     DEFAULT   MEANING
# ---------------------------------------------------------------------------------------------
# X     String    ---      Location to read the matrix X of feature vectors
# Y     String    ---      Location to read response matrix Y with 1 column
# B     String    ---      Location to store estimated regression parameters (the betas)
# S     String    ---      Location to write the selected features ordered as computed by the algorithm
# O     String    " "      Location to write the printed statistics; by default is standard output
# link  Int       2        Link function code: 1 = log, 2 = Logit, 3 = Probit, 4 = Cloglog
# yneg  Double    0.0      Response value for Bernoulli "No" label, usually 0.0 or -1.0
# icpt  Int       0        Intercept presence, X columns shifting and rescaling:
#                          0 = no intercept, no shifting, no rescaling;
#                          1 = add intercept, but neither shift nor rescale X;
#                          2 = add intercept, shift & rescale X columns to mean = 0, variance = 1
# tol   Double    0.000001 Tolerance (epsilon)
# disp  Double    0.0      (Over-)dispersion value, or 0.0 to estimate it from data
# moi   Int       200      Maximum number of outer (Newton / Fisher Scoring) iterations
# mii   Int       0        Maximum number of inner (Conjugate Gradient) iterations, 0 = no maximum
# thr   Double    0.01     Threshold to stop the algorithm: if the decrease in the value of AIC falls below thr
#                          no further features are being checked and the algorithm stops 
# fmt   String   "text"    The betas matrix output format, such as "text" or "csv"
# ---------------------------------------------------------------------------------------------
# OUTPUT: Matrix beta, whose size depends on icpt:
#     icpt=0: ncol(X) x 1;  icpt=1: (ncol(X) + 1) x 1;  icpt=2: (ncol(X) + 1) x 2
#
# In addition, in the last run of GLM some statistics are provided in CSV format, one comma-separated name-value
# pair per each line, as follows:
#
# NAME                  MEANING
# -------------------------------------------------------------------------------------------
# TERMINATION_CODE      A positive integer indicating success/failure as follows:
#                       1 = Converged successfully; 2 = Maximum number of iterations reached; 
#                       3 = Input (X, Y) out of range; 4 = Distribution/link is not supported
# BETA_MIN              Smallest beta value (regression coefficient), excluding the intercept
# BETA_MIN_INDEX        Column index for the smallest beta value
# BETA_MAX              Largest beta value (regression coefficient), excluding the intercept
# BETA_MAX_INDEX        Column index for the largest beta value
# INTERCEPT             Intercept value, or NaN if there is no intercept (if icpt=0)
# DISPERSION            Dispersion used to scale deviance, provided as "disp" input parameter
#                       or estimated (same as DISPERSION_EST) if the "disp" parameter is <= 0
# DISPERSION_EST        Dispersion estimated from the dataset
# DEVIANCE_UNSCALED     Deviance from the saturated model, assuming dispersion == 1.0
# DEVIANCE_SCALED       Deviance from the saturated model, scaled by the DISPERSION value
# -------------------------------------------------------------------------------------------
#
# HOW TO INVOKE THIS SCRIPT - EXAMPLE:
# hadoop jar SystemML.jar -f StepGLM.dml -nvargs X=INPUT_DIR/X Y=INPUT_DIR/Y B=OUTPUT_DIR/betas
# 			 S=OUTPUT_DIR_S/selected O=OUTPUT_DIR/stats link=2 yneg=-1.0 icpt=2 tol=0.00000001 
#			 disp=1.0 moi=100 mii=10 thr=0.01 fmt=csv  
#
# THE StepGLM SCRIPT CURRENTLY SUPPORTS BERNOULLI DISTRIBUTION FAMILY AND THE FOLLOWING LINK FUNCTIONS ONLY!
# 	- LOG
#	- LOGIT
#	- PROBIT
#	- CLOGLOG
	
fileX = $X;
fileY = $Y;
fileB = $B;
intercept_status = ifdef ($icpt, 0);   
thr = ifdef ($thr, 0.01); 
bernoulli_No_label = ifdef ($yneg, 0.0);    # $yneg = 0.0;
distribution_type = 2;

bernoulli_No_label = as.double (bernoulli_No_label);

# currently only the forward selection strategy in supported: start from one feature and iteratively add 
# features until AIC improves
dir = "forward";

print("BEGIN STEPWISE GLM SCRIPT");
print ("Reading X and Y...");
X_orig = read (fileX);
Y = read (fileY);

if (distribution_type == 2 & ncol(Y) == 1) {
	is_Y_negative = (Y == bernoulli_No_label);
	Y = append (1 - is_Y_negative, is_Y_negative);
	count_Y_negative = sum (is_Y_negative);
	if (count_Y_negative == 0) {
		stop ("StepGLM Input Error: all Y-values encode Bernoulli YES-label, none encode NO-label");
	}
	if (count_Y_negative == nrow(Y)) {
		stop ("StepGLM Input Error: all Y-values encode Bernoulli NO-label, none encode YES-label");
	}
}

num_records = nrow (X_orig);
num_features = ncol (X_orig);

# BEGIN STEPWISE GENERALIZED LINEAR MODELS 

if (dir == "forward") {  
	
	continue = TRUE;
	columns_fixed = matrix (0, rows = 1, cols = num_features);
	columns_fixed_ordered = matrix (0, rows = 1, cols = 1);
  
	# X_global stores the best model found at each step 
	X_global = matrix (0, rows = num_records, cols = 1);
 
	if (intercept_status == 0) {
		# Compute AIC of an empty model with no features and no intercept (all Ys are zero)
		[AIC_best] = glm (X_global, Y, 0, num_features, columns_fixed_ordered, " ");
	} else {
		# compute AIC of an empty model with only intercept (all Ys are constant)
		all_ones = matrix (1, rows = num_records, cols = 1);
		[AIC_best] = glm (all_ones, Y, 0, num_features, columns_fixed_ordered, " ");
	}
	print ("Best AIC without any features: " + AIC_best);
  
	# First pass to examine single features
	AICs = matrix (AIC_best, rows = 1, cols = num_features);
	parfor (i in 1:num_features) { 	
		[AIC_1] = glm (X_orig[,i], Y, intercept_status, num_features, columns_fixed_ordered, " ");
		AICs[1,i] = AIC_1;
	}
  
	# Determine the best AIC 
	column_best = 0;	
	for (k in 1:num_features) {
		AIC_cur = as.scalar (AICs[1,k]);
		if ( (AIC_cur < AIC_best) & ((AIC_best - AIC_cur) > abs (thr * AIC_best)) ) {
			column_best = k;
			AIC_best = as.scalar(AICs[1,k]);
		}
	}
  
	if (column_best == 0) {
		print ("AIC of an empty model is " + AIC_best + " and adding no feature achieves more than " + (thr * 100) + "% decrease in AIC!");
		if (intercept_status == 0) {
			# Compute AIC of an empty model with no features and no intercept (all Ys are zero)
			[AIC_best] = glm (X_global, Y, 0, num_features, columns_fixed_ordered, fileB);
		} else {
			# compute AIC of an empty model with only intercept (all Ys are constant)
			###all_ones = matrix (1, rows = num_records, cols = 1);
			[AIC_best] = glm (all_ones, Y, 0, num_features, columns_fixed_ordered, fileB);
		}
	};
  
	print ("Best AIC " + AIC_best + " achieved with feature: " + column_best);	
	columns_fixed[1,column_best] = 1;
	columns_fixed_ordered[1,1] = column_best;
	X_global = X_orig[,column_best];		
  
	while (continue) {
		# Subsequent passes over the features
		parfor (i in 1:num_features) { 
			if (as.scalar(columns_fixed[1,i]) == 0) {	
        
				# Construct the feature matrix
				X = append (X_global, X_orig[,i]);
        
				[AIC_2] = glm (X, Y, intercept_status, num_features, columns_fixed_ordered, " ");
				AICs[1,i] = AIC_2;
			}		
		}
    
		# Determine the best AIC
		for (k in 1:num_features) {
			AIC_cur = as.scalar (AICs[1,k]);
			if ( (AIC_cur < AIC_best) & ((AIC_best - AIC_cur) > abs (thr * AIC_best)) & (as.scalar(columns_fixed[1,k]) == 0) ) {
				column_best = k;
				AIC_best = as.scalar(AICs[1,k]);
			}
		}
    
		# Append best found features (i.e., columns) to X_global
		if (as.scalar(columns_fixed[1,column_best]) == 0) { # new best feature found
			print ("Best AIC " + AIC_best + " achieved with feature: " + column_best);
			columns_fixed[1,column_best] = 1;
			columns_fixed_ordered = append (columns_fixed_ordered, as.matrix(column_best));
			if (ncol(columns_fixed_ordered) == num_features) { # all features examined
				X_global = append (X_global, X_orig[,column_best]);
				continue = FALSE;
			} else {
				X_global = append (X_global, X_orig[,column_best]);
			}
		} else {
		continue = FALSE;
		}
	}
  
	# run GLM with selected set of features
	print ("Running GLM with selected features...");
	[AIC] = glm (X_global, Y, intercept_status, num_features, columns_fixed_ordered, fileB);
  
} else {
	stop ("Currently only forward selection strategy is supported!");
}


################### UDFS USED IN THIS SCRIPT ##################

glm = function (Matrix[Double] X, Matrix[Double] Y, Int intercept_status, Double num_features_orig, Matrix[Double] Selected, String fileB) return (Double AIC) {
		
	# distribution family code: 1 = Power, 2 = Bernoulli/Binomial; currently only Bernouli distribution family is supported!		
	distribution_type = 2;                		# $dfam = 2;
	variance_as_power_of_the_mean = 0.0;  		# $vpow = 0.0;
	# link function code: 0 = canonical (depends on distribution), 1 = Power, 2 = Logit, 3 = Probit, 4 = Cloglog, 5 = Cauchit;
	# currently only log (link = 1), logit (link = 2), probit (link = 3), and cloglog (link = 4) are supported!
	link_type = ifdef ($link, 2);         		# $link = 2;
	link_as_power_of_the_mean = 0.0;      		# $lpow = 0.0;

	dispersion = ifdef ($disp, 0.0);            # $disp = 0.0;
	eps = ifdef ($tol, 0.000001);               # $tol  = 0.000001;
	max_iteration_IRLS = ifdef ($moi, 200);     # $moi  = 200;
	max_iteration_CG = ifdef ($mii, 0);         # $mii  = 0;

	variance_as_power_of_the_mean = as.double (variance_as_power_of_the_mean);
	link_as_power_of_the_mean = as.double (link_as_power_of_the_mean);

	dispersion = as.double (dispersion);
	eps = as.double (eps);              

	# Default values for output statistics:
	regularization = 0.0;
	termination_code     = 0.0;
	min_beta             = 0.0 / 0.0;
	i_min_beta           = 0.0 / 0.0;
	max_beta             = 0.0 / 0.0;
	i_max_beta           = 0.0 / 0.0;
	intercept_value      = 0.0 / 0.0;
	dispersion           = 0.0 / 0.0;
	estimated_dispersion = 0.0 / 0.0;
	deviance_nodisp      = 0.0 / 0.0;
	deviance             = 0.0 / 0.0;				  
                  
	#####   INITIALIZE THE PARAMETERS   #####
                  
    num_records  = nrow (X);
    num_features = ncol (X);
    zeros_r = matrix (0, rows = num_records, cols = 1);
    ones_r = 1 + zeros_r;
                  
    # Introduce the intercept, shift and rescale the columns of X if needed
                  
    if (intercept_status == 1 | intercept_status == 2) { # add the intercept column
		X = append (X, ones_r);
        num_features = ncol (X);
	}
                  
    scale_lambda = matrix (1, rows = num_features, cols = 1);
    if (intercept_status == 1 | intercept_status == 2) {
		scale_lambda [num_features, 1] = 0;
    }
                  
    if (intercept_status == 2) {  # scale-&-shift X columns to mean 0, variance 1
		# Important assumption: X [, num_features] = ones_r
        avg_X_cols = t(colSums(X)) / num_records;
        var_X_cols = (t(colSums (X ^ 2)) - num_records * (avg_X_cols ^ 2)) / (num_records - 1);
        is_unsafe = (var_X_cols <= 0);
        scale_X = 1.0 / sqrt (var_X_cols * (1 - is_unsafe) + is_unsafe);
        scale_X [num_features, 1] = 1;
        shift_X = - avg_X_cols * scale_X;
        shift_X [num_features, 1] = 0;
        rowSums_X_sq = (X ^ 2) %*% (scale_X ^ 2) + X %*% (2 * scale_X * shift_X) + sum (shift_X ^ 2);
	} else {
        scale_X = matrix (1, rows = num_features, cols = 1);
        shift_X = matrix (0, rows = num_features, cols = 1);
        rowSums_X_sq = rowSums (X ^ 2);
    }	
                  
    # Henceforth we replace "X" with "X %*% (SHIFT/SCALE TRANSFORM)" and rowSums(X ^ 2)
    # with "rowSums_X_sq" in order to preserve the sparsity of X under shift and scale.
    # The transform is then associatively applied to the other side of the expression,
    # and is rewritten via "scale_X" and "shift_X" as follows:
    #
    # ssX_A  = (SHIFT/SCALE TRANSFORM) %*% A    --- is rewritten as:
    # ssX_A  = diag (scale_X) %*% A;
    # ssX_A [num_features, ] = ssX_A [num_features, ] + t(shift_X) %*% A;
    #
    # tssX_A = t(SHIFT/SCALE TRANSFORM) %*% A   --- is rewritten as:
    # tssX_A = diag (scale_X) %*% A + shift_X %*% A [num_features, ];
                  
    # Initialize other input-dependent parameters
                  
    lambda = scale_lambda * regularization;
    if (max_iteration_CG == 0) {
		max_iteration_CG = num_features;
    }
                  
    # Set up the canonical link, if requested [Then we have: Var(mu) * (d link / d mu) = const]
                  
    if (link_type == 0) {
		if (distribution_type == 1) {
			link_type = 1;
            link_as_power_of_the_mean = 1.0 - variance_as_power_of_the_mean;
		} else { 
			if (distribution_type == 2) {
				link_type = 2;
			}   
		}   
	}
                                
    # For power distributions and/or links, we use two constants,
    # "variance as power of the mean" and "link_as_power_of_the_mean",
    # to specify the variance and the link as arbitrary powers of the
    # mean.  However, the variance-powers of 1.0 (Poisson family) and
    # 2.0 (Gamma family) have to be treated as special cases, because
    # these values integrate into logarithms.  The link-power of 0.0
    # is also special as it represents the logarithm link.
                  
    num_response_columns = ncol (Y);
    is_supported = 0;              
	if (num_response_columns == 2 & distribution_type == 2 & link_type >= 1 & link_type <= 4) { # BERNOULLI DISTRIBUTION
		is_supported = 1;			  
	}
	if (num_response_columns == 1 & distribution_type == 2) {
		print ("Error: Bernoulli response matrix has not been converted into two-column format.");
    }

	if (is_supported == 1) {
                    
		#####   INITIALIZE THE BETAS   #####
                    
        [beta, saturated_log_l, isNaN] = 
			glm_initialize (X, Y, distribution_type, variance_as_power_of_the_mean, link_type, link_as_power_of_the_mean, intercept_status, max_iteration_CG);
					  
		# print(" --- saturated logLik " + saturated_log_l);
					  
        if (isNaN == 0) {
                      
			#####  START OF THE MAIN PART  #####
                      
            sum_X_sq = sum (rowSums_X_sq);
            trust_delta = 0.5 * sqrt (num_features) / max (sqrt (rowSums_X_sq));
            ###  max_trust_delta = trust_delta * 10000.0;
            log_l = 0.0;
            deviance_nodisp = 0.0;
            new_deviance_nodisp = 0.0;
            isNaN_log_l = 2;
            newbeta = beta;
            g = matrix (0.0, rows = num_features, cols = 1);
            g_norm = sqrt (sum ((g + lambda * beta) ^ 2));
            accept_new_beta = 1;
            reached_trust_boundary = 0;
            neg_log_l_change_predicted = 0.0;
            i_IRLS = 0;
                      
            # print ("BEGIN IRLS ITERATIONS...");
                      
            ssX_newbeta = diag (scale_X) %*% newbeta;
            ssX_newbeta [num_features, ] = ssX_newbeta [num_features, ] + t(shift_X) %*% newbeta;
            all_linear_terms = X %*% ssX_newbeta;
                      
            [new_log_l, isNaN_new_log_l] = glm_log_likelihood_part
				(all_linear_terms, Y, distribution_type, variance_as_power_of_the_mean, link_type, link_as_power_of_the_mean);
					  
            if (isNaN_new_log_l == 0) {
				new_deviance_nodisp = 2.0 * (saturated_log_l - new_log_l);
                new_log_l = new_log_l - 0.5 * sum (lambda * newbeta ^ 2);
            }
                      
            while (termination_code == 0) {
				accept_new_beta = 1;
                        
                if (i_IRLS > 0) {
					if (isNaN_log_l == 0) {
						accept_new_beta = 0;
                    }
                          
                    # Decide whether to accept a new iteration point and update the trust region
                    # See Alg. 4.1 on p. 69 of "Numerical Optimization" 2nd ed. by Nocedal and Wright
                          
                    rho = (- new_log_l + log_l) / neg_log_l_change_predicted;
					if (rho < 0.25 | isNaN_new_log_l == 1) {
						trust_delta = 0.25 * trust_delta;
					}
					if (rho > 0.75 & isNaN_new_log_l == 0 & reached_trust_boundary == 1) {
						trust_delta = 2 * trust_delta;
                            
						### if (trust_delta > max_trust_delta) {
						###     trust_delta = max_trust_delta;
						### }
					}
					if (rho > 0.1 & isNaN_new_log_l == 0) {
						accept_new_beta = 1;
					}
				}
                        
                if (accept_new_beta == 1) {
					beta = newbeta;  log_l = new_log_l;  deviance_nodisp = new_deviance_nodisp;  isNaN_log_l = isNaN_new_log_l;
                          
					[g_Y, w] = glm_dist (all_linear_terms, Y, distribution_type, variance_as_power_of_the_mean, link_type, link_as_power_of_the_mean);
                          
					# We introduced these variables to avoid roundoff errors:
					#     g_Y = y_residual / (y_var * link_grad);
					#     w   = 1.0 / (y_var * link_grad * link_grad);
                          
					gXY = - t(X) %*% g_Y;
					g = diag (scale_X) %*% gXY + shift_X %*% gXY [num_features, ];
					g_norm = sqrt (sum ((g + lambda * beta) ^ 2));
				}
                        
                [z, neg_log_l_change_predicted, num_CG_iters, reached_trust_boundary] = 
					get_CG_Steihaug_point (X, scale_X, shift_X, w, g, beta, lambda, trust_delta, max_iteration_CG);
                        
				newbeta = beta + z;
                        
				ssX_newbeta = diag (scale_X) %*% newbeta;
				ssX_newbeta [num_features, ] = ssX_newbeta [num_features, ] + t(shift_X) %*% newbeta;
				all_linear_terms = X %*% ssX_newbeta;
                        
				[new_log_l, isNaN_new_log_l] = glm_log_likelihood_part
					(all_linear_terms, Y, distribution_type, variance_as_power_of_the_mean, link_type, link_as_power_of_the_mean);
                        
				if (isNaN_new_log_l == 0) {
					new_deviance_nodisp = 2.0 * (saturated_log_l - new_log_l);
					new_log_l = new_log_l - 0.5 * sum (lambda * newbeta ^ 2);
				}
                        
				log_l_change = new_log_l - log_l;               # R's criterion for termination: |dev - devold|/(|dev| + 0.1) < eps
                       
				if (reached_trust_boundary == 0 & isNaN_new_log_l == 0 & 
				   (2.0 * abs (log_l_change) < eps * (deviance_nodisp + 0.1) | abs (log_l_change) < (abs (log_l) + abs (new_log_l)) * 0.00000000000001) ) {
					termination_code = 1;
				}
				rho = - log_l_change / neg_log_l_change_predicted;
				z_norm = sqrt (sum (z * z));
                                              
				i_IRLS = i_IRLS + 1;                        
                        
				if (i_IRLS == max_iteration_IRLS) {
					termination_code = 2;
				}
			}
                      
            beta = newbeta;
            log_l = new_log_l;
            deviance_nodisp = new_deviance_nodisp;
                     
            #---------------------------- last part

			if (termination_code != 1) {
				print ("One of the runs of GLM did not converged in " + i_IRLS + " steps!");
			}
                      
            ##### COMPUTE AIC ##### 
            
			if (distribution_type == 2 & link_type >= 1 & link_type <= 4) {			
				AIC = -2 * log_l;
				if (sum (X) != 0) {
					AIC = AIC + 2 * num_features;	
				}
			} else {
				stop ("Currently only the Bernoulli distribution family the following link functions are supported: log, logit, probit, and cloglog!");
			}
					           
            if (fileB != " ") {
				fileO = ifdef ($O, " ");
				fileS = $S;
				fmt  = ifdef ($fmt, "text");	
			
				# Output which features give the best AIC and are being used for linear regression 
				write (Selected, fileS, format=fmt);
		
				ssX_beta = diag (scale_X) %*% beta;
                ssX_beta [num_features, ] = ssX_beta [num_features, ] + t(shift_X) %*% beta;
				if (intercept_status == 2) {
					beta_out = append (ssX_beta, beta);
                } else {
					beta_out = ssX_beta;
                }
                
                if (intercept_status == 0 & num_features == 1) {
					p = sum (X == 1);
					if (p == num_records) {
						beta_out = beta_out[1,];
					}					
                } 

								
                if (intercept_status == 1 | intercept_status == 2) {
					intercept_value = as.scalar (beta_out [num_features, 1]);
                    beta_noicept = beta_out [1 : (num_features - 1), 1];
                } else {
					beta_noicept = beta_out [1 : num_features, 1];
                }
                min_beta = min (beta_noicept);
                max_beta = max (beta_noicept);
                tmp_i_min_beta = rowIndexMin (t(beta_noicept))
                i_min_beta = as.scalar (tmp_i_min_beta [1, 1]);
                tmp_i_max_beta = rowIndexMax (t(beta_noicept))
                i_max_beta = as.scalar (tmp_i_max_beta [1, 1]);
                        
                #####  OVER-DISPERSION PART  #####
                      
                all_linear_terms = X %*% ssX_beta;
				[g_Y, w] = glm_dist (all_linear_terms, Y, distribution_type, variance_as_power_of_the_mean, link_type, link_as_power_of_the_mean);
                        
                pearson_residual_sq = g_Y ^ 2 / w;
                pearson_residual_sq = replace (target = pearson_residual_sq, pattern = 0.0/0.0, replacement = 0);
                # pearson_residual_sq = (y_residual ^ 2) / y_var;
                        
                if (num_records > num_features) {
					estimated_dispersion = sum (pearson_residual_sq) / (num_records - num_features);
                }
                if (dispersion <= 0) {
					dispersion = estimated_dispersion;
                }
                deviance = deviance_nodisp / dispersion;
                        
                #####  END OF THE MAIN PART  #####
                        
                str = "BETA_MIN," + min_beta;
                str = append (str, "BETA_MIN_INDEX," + i_min_beta);
                str = append (str, "BETA_MAX," + max_beta);
                str = append (str, "BETA_MAX_INDEX," + i_max_beta);
                str = append (str, "INTERCEPT," + intercept_value);
                str = append (str, "DISPERSION," + dispersion);
                str = append (str, "DISPERSION_EST," + estimated_dispersion);
                str = append (str, "DEVIANCE_UNSCALED," + deviance_nodisp);
                str = append (str, "DEVIANCE_SCALED," + deviance);
                        
                if (fileO != " ") {
					write (str, fileO);
                } 
				else {
					print (str);
                }
			
				# Prepare the output matrix
				print ("Writing the output matrix...");
                if (intercept_status == 0 & num_features == 1) { 
					if (p == num_records) {
						beta_out_tmp = matrix (0, rows = num_features_orig + 1, cols = 1); 
						beta_out_tmp[num_features_orig + 1,] = beta_out;
						beta_out = beta_out_tmp;
						write (beta_out, fileB, format=fmt);
						stop ("");
					} else if (sum (X) == 0){
						beta_out = matrix (0, rows = num_features_orig, cols = 1);
						write (beta_out, fileB, format=fmt);
						stop ("");
					}
				}

				no_selected = ncol (Selected);
				max_selected = max (Selected);
				last = max_selected + 1;	
		
				if (intercept_status != 0) {
		
					Selected_ext = append (Selected, as.matrix (last));			
					P1 = table (seq (1, ncol (Selected_ext)), t(Selected_ext)); 

					if (intercept_status == 2) {
			
						P1_ssX_beta = P1 * ssX_beta;
						P2_ssX_beta = colSums (P1_ssX_beta);
						P1_beta = P1 * beta;
						P2_beta = colSums (P1_beta);
				
						if (max_selected < num_features_orig) {
						
							P2_ssX_beta = append (P2_ssX_beta, matrix (0, rows=1, cols=(num_features_orig - max_selected)));
							P2_beta = append (P2_beta, matrix (0, rows=1, cols=(num_features_orig - max_selected)));
							
							P2_ssX_beta[1, num_features_orig+1] = P2_ssX_beta[1, max_selected + 1]; 
							P2_ssX_beta[1, max_selected + 1] = 0;
							
							P2_beta[1, num_features_orig+1] = P2_beta[1, max_selected + 1]; 
							P2_beta[1, max_selected + 1] = 0;

						}
						beta_out = append (t(P2_ssX_beta), t(P2_beta));
				
					} else {
			
						P1_beta = P1 * beta;
						P2_beta = colSums (P1_beta);
				
						if (max_selected < num_features_orig) {
							P2_beta = append (P2_beta, matrix (0, rows=1, cols=(num_features_orig - max_selected)));
							P2_beta[1, num_features_orig+1] = P2_beta[1, max_selected + 1] ; 
							P2_beta[1, max_selected + 1] = 0;
						}
						beta_out = t(P2_beta);
				
					}
				} else {
		
					P1 = table (seq (1, no_selected), t(Selected)); 
					P1_beta = P1 * beta;
					P2_beta = colSums (P1_beta);	

					if (max_selected < num_features_orig) {
						P2_beta = append (P2_beta, matrix (0, rows=1, cols=(num_features_orig - max_selected)));
					}		

					beta_out = t(P2_beta);	
				}
	
				write ( beta_out, fileB, format=fmt );
			
			}
                      
		} else { 
			stop ("Input matrices X and/or Y are out of range!"); 
        }
	} else { 
		stop ("Response matrix with " + num_response_columns + " columns, distribution family (" + distribution_type + ", " + variance_as_power_of_the_mean
               + ") and link family (" + link_type + ", " + link_as_power_of_the_mean + ") are NOT supported together.");
    }
}

glm_initialize = function (Matrix[double] X, Matrix[double] Y, int dist_type, double var_power, int link_type, double link_power, int icept_status, int max_iter_CG)
  return (Matrix[double] beta, double saturated_log_l, int isNaN)
{
    saturated_log_l = 0.0;
    isNaN = 0;
    y_corr = Y [, 1];
    if (dist_type == 2) {
      n_corr = rowSums (Y);
      is_n_zero = (n_corr == 0);
      y_corr = Y [, 1] / (n_corr + is_n_zero) + (0.5 - Y [, 1]) * is_n_zero;    
    }
    linear_terms = y_corr;
    if (dist_type == 1 & link_type == 1) { # POWER DISTRIBUTION
      if          (link_power ==  0) {
        if (sum (y_corr < 0) == 0) {
          is_zero_y_corr = (y_corr == 0);
          linear_terms = log (y_corr + is_zero_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        } else { isNaN = 1; }
      } else { if (link_power ==  1.0) {
        linear_terms = y_corr;
      } else { if (link_power == -1.0) {
        linear_terms = 1.0 / y_corr;
      } else { if (link_power ==  0.5) {
        if (sum (y_corr < 0) == 0) {
          linear_terms = sqrt (y_corr);
        } else { isNaN = 1; }
      } else { if (link_power >   0) {
        if (sum (y_corr < 0) == 0) {
          is_zero_y_corr = (y_corr == 0);
          linear_terms = (y_corr + is_zero_y_corr) ^ link_power - is_zero_y_corr;
        } else { isNaN = 1; }
      } else {
        if (sum (y_corr <= 0) == 0) {
          linear_terms = y_corr ^ link_power;
        } else { isNaN = 1; }
      }}}}}
    }
    if (dist_type == 2 & link_type >= 1 & link_type <= 5)
    { # BINOMIAL/BERNOULLI DISTRIBUTION
      if          (link_type == 1 & link_power == 0)  { # Binomial.log
        if (sum (y_corr < 0) == 0) {
          is_zero_y_corr = (y_corr == 0);
          linear_terms = log (y_corr + is_zero_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        } else { isNaN = 1; }
      } else { if (link_type == 1 & link_power >  0)  { # Binomial.power_nonlog pos
        if (sum (y_corr < 0) == 0) {
          is_zero_y_corr = (y_corr == 0);
          linear_terms = (y_corr + is_zero_y_corr) ^ link_power - is_zero_y_corr;
        } else { isNaN = 1; }
      } else { if (link_type == 1)                      { # Binomial.power_nonlog neg
        if (sum (y_corr <= 0) == 0) {
          linear_terms = y_corr ^ link_power;
        } else { isNaN = 1; }
      } else { 
        is_zero_y_corr = (y_corr <= 0);
        is_one_y_corr  = (y_corr >= 1.0);
        y_corr = y_corr * (1.0 - is_zero_y_corr) * (1.0 - is_one_y_corr) + 0.5 * (is_zero_y_corr + is_one_y_corr);
        if (link_type == 2)                           { # Binomial.logit
          linear_terms = log (y_corr / (1.0 - y_corr)) 
          + is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        } else { if (link_type == 3)                  { # Binomial.probit
          y_below_half = y_corr + (1.0 - 2.0 * y_corr) * (y_corr > 0.5);
          t = sqrt (- 2.0 * log (y_below_half));
          approx_inv_Gauss_CDF = - t + (2.515517 + t * (0.802853 + t * 0.010328)) / (1.0 + t * (1.432788 + t * (0.189269 + t * 0.001308)));
          linear_terms = approx_inv_Gauss_CDF * (1.0 - 2.0 * (y_corr > 0.5))
          + is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        } else { if (link_type == 4)                  { # Binomial.cloglog
          linear_terms = log (- log (1.0 - y_corr))
          - log (- log (0.5)) * (is_zero_y_corr + is_one_y_corr)
          + is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        } else { if (link_type == 5)                  { # Binomial.cauchit
          linear_terms = tan ((y_corr - 0.5) * 3.1415926535897932384626433832795)
          + is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
        }}  }}}}}
    }
    
    if (isNaN == 0) {
      [saturated_log_l, isNaN] = 
        glm_log_likelihood_part (linear_terms, Y, dist_type, var_power, link_type, link_power);
    }
    
    if ((dist_type == 1 & link_type == 1 & link_power == 0) |
          (dist_type == 2 & link_type >= 2))
    {    
      desired_eta = 0.0;
    } else { if (link_type == 1 & link_power == 0) {
      desired_eta = log (0.5);
    } else { if (link_type == 1) {
      desired_eta = 0.5 ^ link_power;
    } else {
      desired_eta = 0.5;
    }}}
    
    beta = matrix (0.0, rows = ncol(X), cols = 1);
    
    if (desired_eta != 0) {
      if (icept_status == 1 | icept_status == 2) {
        beta [nrow(beta), 1] = desired_eta;
      } else {
        # We want: avg (X %*% ssX_transform %*% beta) = desired_eta
        # Note that "ssX_transform" is trivial here, hence ignored
        
        beta = straightenX (X, 0.000001, max_iter_CG);  
        beta = beta * desired_eta;
      }   }   }


glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
                     int dist_type, double var_power, int link_type, double link_power)
  return (Matrix[double] g_Y, Matrix[double] w)
# ORIGINALLY we returned more meaningful vectors, namely:
# Matrix[double] y_residual    : y - y_mean, i.e. y observed - y predicted
# Matrix[double] link_gradient : derivative of the link function
# Matrix[double] var_function  : variance without dispersion, i.e. the V(mu) function
# BUT, this caused roundoff errors, so we had to compute "directly useful" vectors
# and skip over the "meaningful intermediaries".  Now we output these two variables:
#     g_Y = y_residual / (var_function * link_gradient);
#     w   = 1.0 / (var_function * link_gradient ^ 2);
{
    num_records = nrow (linear_terms);
    zeros_r = matrix (0.0, rows = num_records, cols = 1);
    ones_r = 1 + zeros_r;
    g_Y  = zeros_r;
    w  = zeros_r;
    
    # Some constants
    
    one_over_sqrt_two_pi = 0.39894228040143267793994605993438;
    ones_2 = matrix (1.0, rows = 1, cols = 2);
    p_one_m_one = ones_2;
    p_one_m_one [1, 2] = -1.0;
    m_one_p_one = ones_2;
    m_one_p_one [1, 1] = -1.0;
    zero_one = ones_2;
    zero_one [1, 1] = 0.0;
    one_zero = ones_2;
    one_zero [1, 2] = 0.0;
    flip_pos = matrix (0, rows = 2, cols = 2);
    flip_neg = flip_pos;
    flip_pos [1, 2] = 1;
    flip_pos [2, 1] = 1;
    flip_neg [1, 2] = -1;
    flip_neg [2, 1] = 1;
    
    if (dist_type == 1 & link_type == 1) { # POWER DISTRIBUTION
      y_mean = zeros_r;
      if          (link_power ==  0) {
        y_mean = exp (linear_terms);
        y_mean_pow = y_mean ^ (1 - var_power);
        w   = y_mean_pow * y_mean;
        g_Y = y_mean_pow * (Y - y_mean);
      } else { if (link_power ==  1.0) {
        y_mean = linear_terms;
        w   = y_mean ^ (- var_power);
        g_Y = w * (Y - y_mean);
      } else {
        y_mean = linear_terms ^ (1.0 / link_power);
        c1  = (1 - var_power) / link_power - 1;
        c2  = (2 - var_power) / link_power - 2;
        g_Y = (linear_terms ^ c1) * (Y - y_mean) / link_power;
        w   = (linear_terms ^ c2) / (link_power ^ 2);
      }   }}
    if (dist_type == 2 & link_type >= 1 & link_type <= 5)
    { # BINOMIAL/BERNOULLI DISTRIBUTION
      if (link_type == 1) { # BINOMIAL.POWER LINKS
        if (link_power == 0)  { # Binomial.log
          vec1 = 1 / (exp (- linear_terms) - 1);
          g_Y = Y [, 1] - Y [, 2] * vec1;
          w   = rowSums (Y) * vec1;
        } else {                  # Binomial.nonlog
          vec1 = zeros_r;
          if (link_power == 0.5)  {
            vec1 = 1 / (1 - linear_terms ^ 2);
          } else { if (sum (linear_terms < 0) == 0) {
            vec1 = linear_terms ^ (- 2 + 1 / link_power) / (1 - linear_terms ^ (1 / link_power));
          } else {isNaN = 1;}}
          # We want a "zero-protected" version of
          #     vec2 = Y [, 1] / linear_terms;
          is_y_0 = (Y [, 1] == 0);
          vec2 = (Y [, 1] + is_y_0) / (linear_terms * (1 - is_y_0) + is_y_0) - is_y_0;
          g_Y =  (vec2 - Y [, 2] * vec1 * linear_terms) / link_power;
          w   =  rowSums (Y) * vec1 / link_power ^ 2;
        }
      } else {
        is_LT_pos_infinite = (linear_terms == 1.0/0.0);
        is_LT_neg_infinite = (linear_terms == -1.0/0.0);
        is_LT_infinite = is_LT_pos_infinite %*% one_zero + is_LT_neg_infinite %*% zero_one;
        finite_linear_terms = replace (target =        linear_terms, pattern =  1.0/0.0, replacement = 0);
        finite_linear_terms = replace (target = finite_linear_terms, pattern = -1.0/0.0, replacement = 0);
        if (link_type == 2)                           { # Binomial.logit
          Y_prob = exp (finite_linear_terms) %*% one_zero + ones_r %*% zero_one;
          Y_prob = Y_prob / (rowSums (Y_prob) %*% ones_2);
          Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
          g_Y = rowSums (Y * (Y_prob %*% flip_neg));           ### = y_residual;
          w   = rowSums (Y * (Y_prob %*% flip_pos) * Y_prob);  ### = y_variance;
        } else { if (link_type == 3)                  { # Binomial.probit
          is_lt_pos = (linear_terms >= 0);
          t_gp = 1.0 / (1.0 + abs (finite_linear_terms) * 0.231641888);  # 0.231641888 = 0.3275911 / sqrt (2.0)
          pt_gp = t_gp * ( 0.254829592 
                           + t_gp * (-0.284496736 # "Handbook of Mathematical Functions", ed. by M. Abramowitz and I.A. Stegun,
                                     + t_gp * ( 1.421413741 # U.S. Nat-l Bureau of Standards, 10th print (Dec 1972), Sec. 7.1.26, p. 299
                                                + t_gp * (-1.453152027 
                                                          + t_gp *   1.061405429))));
          the_gauss_exp = exp (- (linear_terms ^ 2) / 2.0);
          vec1 = 0.25 * pt_gp * (2 - the_gauss_exp * pt_gp);
          vec2 = Y [, 1] - rowSums (Y) * is_lt_pos + the_gauss_exp * pt_gp * rowSums (Y) * (is_lt_pos - 0.5);
          w   = the_gauss_exp * (one_over_sqrt_two_pi ^ 2) * rowSums (Y) / vec1;
          g_Y = one_over_sqrt_two_pi * vec2 / vec1;
        } else { if (link_type == 4)                  { # Binomial.cloglog
          the_exp = exp (linear_terms)
          the_exp_exp = exp (- the_exp);
          is_too_small = ((10000000 + the_exp) == 10000000);
          the_exp_ratio = (1 - is_too_small) * (1 - the_exp_exp) / (the_exp + is_too_small) + is_too_small * (1 - the_exp / 2);
          g_Y =  (rowSums (Y) * the_exp_exp - Y [, 2]) / the_exp_ratio;
          w   =  the_exp_exp * the_exp * rowSums (Y) / the_exp_ratio;
        } else { if (link_type == 5)                  { # Binomial.cauchit
          Y_prob = 0.5 + (atan (finite_linear_terms) %*% p_one_m_one) / 3.1415926535897932384626433832795;
          Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
          y_residual = Y [, 1] * Y_prob [, 2] - Y [, 2] * Y_prob [, 1];
          var_function = rowSums (Y) * Y_prob [, 1] * Y_prob [, 2];
          link_gradient_normalized = (1 + linear_terms ^ 2) * 3.1415926535897932384626433832795;
          g_Y =  rowSums (Y) * y_residual / (var_function * link_gradient_normalized);
          w   = (rowSums (Y) ^ 2) / (var_function * link_gradient_normalized ^ 2);
        }}}}   
      }
    }
  }


glm_log_likelihood_part = function (Matrix[double] linear_terms, Matrix[double] Y,
                                    int dist_type, double var_power, int link_type, double link_power)
  return (double log_l, int isNaN)
{
    isNaN = 0;
    log_l = 0.0;
    num_records = nrow (Y);
    zeros_r = matrix (0.0, rows = num_records, cols = 1);
    
    if (dist_type == 1 & link_type == 1)
    { # POWER DISTRIBUTION
      b_cumulant = zeros_r;
      natural_parameters = zeros_r;
      is_natural_parameter_log_zero = zeros_r;
      if          (var_power == 1.0 & link_power == 0)  { # Poisson.log
        b_cumulant = exp (linear_terms);
        is_natural_parameter_log_zero = (linear_terms == -1.0/0.0);
        natural_parameters = replace (target = linear_terms, pattern = -1.0/0.0, replacement = 0);
      } else { if (var_power == 1.0 & link_power == 1.0)  { # Poisson.id
        if (sum (linear_terms < 0) == 0)  {
          b_cumulant = linear_terms;
          is_natural_parameter_log_zero = (linear_terms == 0);
          natural_parameters = log (linear_terms + is_natural_parameter_log_zero);
        } else {isNaN = 1;}
      } else { if (var_power == 1.0 & link_power == 0.5)  { # Poisson.sqrt
        if (sum (linear_terms < 0) == 0)  {
          b_cumulant = linear_terms ^ 2;
          is_natural_parameter_log_zero = (linear_terms == 0);
          natural_parameters = 2.0 * log (linear_terms + is_natural_parameter_log_zero);
        } else {isNaN = 1;}
      } else { if (var_power == 1.0 & link_power  > 0)  { # Poisson.power_nonlog, pos
        if (sum (linear_terms < 0) == 0)  {
          is_natural_parameter_log_zero = (linear_terms == 0);
          b_cumulant = (linear_terms + is_natural_parameter_log_zero) ^ (1.0 / link_power) - is_natural_parameter_log_zero;
          natural_parameters = log (linear_terms + is_natural_parameter_log_zero) / link_power;
        } else {isNaN = 1;}
      } else { if (var_power == 1.0)                      { # Poisson.power_nonlog, neg
        if (sum (linear_terms <= 0) == 0) {
          b_cumulant = linear_terms ^ (1.0 / link_power);
          natural_parameters = log (linear_terms) / link_power;
        } else {isNaN = 1;}
      } else { if (var_power == 2.0 & link_power == -1.0) { # Gamma.inverse
        if (sum (linear_terms <= 0) == 0) {
          b_cumulant = - log (linear_terms);
          natural_parameters = - linear_terms;
        } else {isNaN = 1;}
      } else { if (var_power == 2.0 & link_power ==  1.0) { # Gamma.id
        if (sum (linear_terms <= 0) == 0) {
          b_cumulant = log (linear_terms);
          natural_parameters = - 1.0 / linear_terms;
        } else {isNaN = 1;}
      } else { if (var_power == 2.0 & link_power ==  0) { # Gamma.log
        b_cumulant = linear_terms;
        natural_parameters = - exp (- linear_terms);
      } else { if (var_power == 2.0)                      { # Gamma.power_nonlog
        if (sum (linear_terms <= 0) == 0) {
          b_cumulant = log (linear_terms) / link_power;
          natural_parameters = - linear_terms ^ (- 1.0 / link_power);
        } else {isNaN = 1;}
      } else { if                    (link_power ==  0) { # PowerDist.log
        natural_parameters = exp (linear_terms * (1.0 - var_power)) / (1.0 - var_power);
        b_cumulant = exp (linear_terms * (2.0 - var_power)) / (2.0 - var_power);
      } else {                                              # PowerDist.power_nonlog
        if          (-2 * link_power == 1.0 - var_power) {
          natural_parameters = 1.0 / (linear_terms ^ 2) / (1.0 - var_power);
        } else { if (-1 * link_power == 1.0 - var_power) {
          natural_parameters = 1.0 / linear_terms / (1.0 - var_power);
        } else { if (     link_power == 1.0 - var_power) {
          natural_parameters = linear_terms / (1.0 - var_power);
        } else { if ( 2 * link_power == 1.0 - var_power) {
          natural_parameters = linear_terms ^ 2 / (1.0 - var_power);
        } else {
          if (sum (linear_terms <= 0) == 0) {
            power = (1.0 - var_power) / link_power;
            natural_parameters = (linear_terms ^ power) / (1.0 - var_power);
          } else {isNaN = 1;}
        }}}}
        if          (-2 * link_power == 2.0 - var_power) {
          b_cumulant = 1.0 / (linear_terms ^ 2) / (2.0 - var_power);
        } else { if (-1 * link_power == 2.0 - var_power) {
          b_cumulant = 1.0 / linear_terms / (2.0 - var_power);
        } else { if (     link_power == 2.0 - var_power) {
          b_cumulant = linear_terms / (2.0 - var_power);
        } else { if ( 2 * link_power == 2.0 - var_power) {
          b_cumulant = linear_terms ^ 2 / (2.0 - var_power);
        } else {
          if (sum (linear_terms <= 0) == 0) {
            power = (2.0 - var_power) / link_power;
            b_cumulant = (linear_terms ^ power) / (2.0 - var_power);
          } else {isNaN = 1;}
        }}}}
      }}}}} }}}}}
      if (sum (is_natural_parameter_log_zero * abs (Y)) > 0) {
        log_l = -1.0 / 0.0;
        isNaN = 1;
      }
      if (isNaN == 0)
      {
        log_l = sum (Y * natural_parameters - b_cumulant);
        if (log_l != log_l | (log_l == log_l + 1.0 & log_l == log_l * 2.0)) {
          log_l = -1.0 / 0.0;
          isNaN = 1;
        }   }   }
    
    if (dist_type == 2 & link_type >= 1 & link_type <= 5)
    { # BINOMIAL/BERNOULLI DISTRIBUTION
      
      [Y_prob, isNaN] = binomial_probability_two_column (linear_terms, link_type, link_power);
      
      if (isNaN == 0) {            
        does_prob_contradict = (Y_prob <= 0);
        if (sum (does_prob_contradict * abs (Y)) == 0) {
          log_l = sum (Y * log (Y_prob * (1 - does_prob_contradict) + does_prob_contradict));
          if (log_l != log_l | (log_l == log_l + 1.0 & log_l == log_l * 2.0)) {
            isNaN = 1;
          }
        } else {
          log_l = -1.0 / 0.0;
          isNaN = 1;
        }   }   }
    
    if (isNaN == 1) {
      log_l = - 1.0 / 0.0; 
    }
  }


binomial_probability_two_column =
  function (Matrix[double] linear_terms, int link_type, double link_power)
    return   (Matrix[double] Y_prob, int isNaN)
{
      isNaN = 0;
      num_records = nrow (linear_terms);
      
      # Define some auxiliary matrices
      
      ones_2 = matrix (1.0, rows = 1, cols = 2);
      p_one_m_one = ones_2;
      p_one_m_one [1, 2] = -1.0;
      m_one_p_one = ones_2;
      m_one_p_one [1, 1] = -1.0;
      zero_one = ones_2;
      zero_one [1, 1] = 0.0;
      one_zero = ones_2;
      one_zero [1, 2] = 0.0;
      
      zeros_r = matrix (0.0, rows = num_records, cols = 1);
      ones_r = 1.0 + zeros_r;
      
      # Begin the function body
      
      Y_prob = zeros_r %*% ones_2;
      if (link_type == 1) { # Binomial.power
        if          (link_power == 0) { # Binomial.log
          Y_prob = exp (linear_terms) %*% p_one_m_one + ones_r %*% zero_one;    
        } else { if (link_power == 0.5) { # Binomial.sqrt
          Y_prob = (linear_terms ^ 2) %*% p_one_m_one + ones_r %*% zero_one;    
        } else {                          # Binomial.power_nonlog
          if (sum (linear_terms < 0) == 0) {
            Y_prob = (linear_terms ^ (1.0 / link_power)) %*% p_one_m_one + ones_r %*% zero_one;    
          } else {isNaN = 1;}
        }}
      } else {              # Binomial.non_power
        is_LT_pos_infinite = (linear_terms == 1.0/0.0);
        is_LT_neg_infinite = (linear_terms == -1.0/0.0);
        is_LT_infinite = is_LT_pos_infinite %*% one_zero + is_LT_neg_infinite %*% zero_one;
        finite_linear_terms = replace (target =        linear_terms, pattern =  1.0/0.0, replacement = 0);
        finite_linear_terms = replace (target = finite_linear_terms, pattern = -1.0/0.0, replacement = 0);
        if (link_type == 2)             { # Binomial.logit
          Y_prob = exp (finite_linear_terms) %*% one_zero + ones_r %*% zero_one;
          Y_prob = Y_prob / (rowSums (Y_prob) %*% ones_2);
        } else { if (link_type == 3)    { # Binomial.probit
          lt_pos_neg = (finite_linear_terms >= 0) %*% p_one_m_one + ones_r %*% zero_one;
          t_gp = 1.0 / (1.0 + abs (finite_linear_terms) * 0.231641888);  # 0.231641888 = 0.3275911 / sqrt (2.0)
          pt_gp = t_gp * ( 0.254829592 
                           + t_gp * (-0.284496736 # "Handbook of Mathematical Functions", ed. by M. Abramowitz and I.A. Stegun,
                                     + t_gp * ( 1.421413741 # U.S. Nat-l Bureau of Standards, 10th print (Dec 1972), Sec. 7.1.26, p. 299
                                                + t_gp * (-1.453152027 
                                                          + t_gp *   1.061405429))));
          the_gauss_exp = exp (- (finite_linear_terms ^ 2) / 2.0);
          Y_prob = lt_pos_neg + ((the_gauss_exp * pt_gp) %*% ones_2) * (0.5 - lt_pos_neg);
        } else { if (link_type == 4)    { # Binomial.cloglog
          the_exp = exp (finite_linear_terms);
          the_exp_exp = exp (- the_exp);
          is_too_small = ((10000000 + the_exp) == 10000000);
          Y_prob [, 1] = (1 - is_too_small) * (1 - the_exp_exp) + is_too_small * the_exp * (1 - the_exp / 2);
          Y_prob [, 2] = the_exp_exp;
        } else { if (link_type == 5)    { # Binomial.cauchit
          Y_prob = 0.5 + (atan (finite_linear_terms) %*% p_one_m_one) / 3.1415926535897932384626433832795;
        } else {
          isNaN = 1;
        }}}}
        Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
      }   }            


# THE CG-STEIHAUG PROCEDURE SCRIPT

# Apply Conjugate Gradient - Steihaug algorithm in order to approximately minimize
# 0.5 z^T (X^T diag(w) X + diag (lambda)) z + (g + lambda * beta)^T z
# under constraint:  ||z|| <= trust_delta.
# See Alg. 7.2 on p. 171 of "Numerical Optimization" 2nd ed. by Nocedal and Wright
# IN THE ABOVE, "X" IS UNDERSTOOD TO BE "X %*% (SHIFT/SCALE TRANSFORM)"; this transform
# is given separately because sparse "X" may become dense after applying the transform.
#
get_CG_Steihaug_point =
  function (Matrix[double] X, Matrix[double] scale_X, Matrix[double] shift_X, Matrix[double] w,
            Matrix[double] g, Matrix[double] beta, Matrix[double] lambda, double trust_delta, int max_iter_CG)
    return (Matrix[double] z, double neg_log_l_change, int i_CG, int reached_trust_boundary)
{
      trust_delta_sq = trust_delta ^ 2;
      size_CG = nrow (g);
      z = matrix (0.0, rows = size_CG, cols = 1);
      neg_log_l_change = 0.0;
      reached_trust_boundary = 0;
      g_reg = g + lambda * beta;
      r_CG = g_reg;
      p_CG = -r_CG;
      rr_CG = sum(r_CG * r_CG);
      eps_CG = rr_CG * min (0.25, sqrt (rr_CG));
      converged_CG = 0;
      if (rr_CG < eps_CG) {
        converged_CG = 1;
      }
      
      max_iteration_CG = max_iter_CG;
      if (max_iteration_CG <= 0) {
        max_iteration_CG = size_CG;
      }
      i_CG = 0;
      while (converged_CG == 0)
      {
        i_CG = i_CG + 1;
        ssX_p_CG = diag (scale_X) %*% p_CG;
        ssX_p_CG [size_CG, ] = ssX_p_CG [size_CG, ] + t(shift_X) %*% p_CG;
        temp_CG = t(X) %*% (w * (X %*% ssX_p_CG));
        q_CG = (lambda * p_CG) + diag (scale_X) %*% temp_CG + shift_X %*% temp_CG [size_CG, ];
        pq_CG = sum (p_CG * q_CG);
        if (pq_CG <= 0) {
          pp_CG = sum (p_CG * p_CG);  
          if (pp_CG > 0) {
            [z, neg_log_l_change] = 
              get_trust_boundary_point (g_reg, z, p_CG, q_CG, r_CG, pp_CG, pq_CG, trust_delta_sq);
            reached_trust_boundary = 1;
          } else {
            neg_log_l_change = 0.5 * sum (z * (r_CG + g_reg));
          }
          converged_CG = 1;
        }
        if (converged_CG == 0) {
          alpha_CG = rr_CG / pq_CG;
          new_z = z + alpha_CG * p_CG;
          if (sum(new_z * new_z) >= trust_delta_sq) {
            pp_CG = sum (p_CG * p_CG);  
            [z, neg_log_l_change] = 
              get_trust_boundary_point (g_reg, z, p_CG, q_CG, r_CG, pp_CG, pq_CG, trust_delta_sq);
            reached_trust_boundary = 1;
            converged_CG = 1;
          }
          if (converged_CG == 0) {
            z = new_z;
            old_rr_CG = rr_CG;
            r_CG = r_CG + alpha_CG * q_CG;
            rr_CG = sum(r_CG * r_CG);
            if (i_CG == max_iteration_CG | rr_CG < eps_CG) {
              neg_log_l_change = 0.5 * sum (z * (r_CG + g_reg));
              reached_trust_boundary = 0;
              converged_CG = 1;
            }
            if (converged_CG == 0) {
              p_CG = -r_CG + (rr_CG / old_rr_CG) * p_CG;
            }   }   }   }   }


# An auxiliary function used twice inside the CG-STEIHAUG loop:
get_trust_boundary_point = 
  function (Matrix[double] g, Matrix[double] z, Matrix[double] p, 
            Matrix[double] q, Matrix[double] r, double pp, double pq, 
            double trust_delta_sq)
    return (Matrix[double] new_z, double f_change)
{
      zz = sum (z * z);  pz = sum (p * z);
      sq_root_d = sqrt (pz * pz - pp * (zz - trust_delta_sq));
      tau_1 = (- pz + sq_root_d) / pp;
      tau_2 = (- pz - sq_root_d) / pp;
      zq = sum (z * q);  gp = sum (g * p);
      f_extra = 0.5 * sum (z * (r + g));
      f_change_1 = f_extra + (0.5 * tau_1 * pq + zq + gp) * tau_1;
      f_change_2 = f_extra + (0.5 * tau_2 * pq + zq + gp) * tau_2;
      ind1 = as.integer(f_change_1 < f_change_2);
      ind2 = as.integer(f_change_1 >= f_change_2);
      new_z = z + ((ind1 * tau_1 + ind2 * tau_2) * p);
      f_change = ind1 * f_change_1 + ind2 * f_change_2;
}


# Computes vector w such that  ||X %*% w - 1|| -> MIN  given  avg(X %*% w) = 1
# We find z_LS such that ||X %*% z_LS - 1|| -> MIN unconditionally, then scale
# it to compute  w = c * z_LS  such that  sum(X %*% w) = nrow(X).
straightenX =
  function (Matrix[double] X, double eps, int max_iter_CG)
    return   (Matrix[double] w)
{
      w_X = t(colSums(X));
      lambda_LS = 0.000001 * sum(X ^ 2) / ncol(X);
      eps_LS = eps * nrow(X);
      
      # BEGIN LEAST SQUARES
      
      r_LS = - w_X;
      z_LS = matrix (0.0, rows = ncol(X), cols = 1);
      p_LS = - r_LS;
      norm_r2_LS = sum (r_LS ^ 2);
      i_LS = 0;
      while (i_LS < max_iter_CG & i_LS < ncol(X) & norm_r2_LS >= eps_LS)
      {
        q_LS = t(X) %*% X %*% p_LS;
        q_LS = q_LS + lambda_LS * p_LS;
        alpha_LS = norm_r2_LS / sum (p_LS * q_LS);
        z_LS = z_LS + alpha_LS * p_LS;
        old_norm_r2_LS = norm_r2_LS;
        r_LS = r_LS + alpha_LS * q_LS;
        norm_r2_LS = sum (r_LS ^ 2);
        p_LS = -r_LS + (norm_r2_LS / old_norm_r2_LS) * p_LS;
        i_LS = i_LS + 1;
      }
      
      # END LEAST SQUARES
      
      w = (nrow(X) / sum (w_X * z_LS)) * z_LS;
    }