You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@systemml.apache.org by du...@apache.org on 2016/01/22 17:34:01 UTC
[25/51] [partial] incubator-systemml git commit: [SYSTEMML-482]
[SYSTEMML-480] Adding a Git attributes file to enfore Unix-styled line
endings, and normalizing all of the line endings.
http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/05d2c0a8/src/test/scripts/applications/impute/imputeGaussMCMC.nogradient.dml
----------------------------------------------------------------------
diff --git a/src/test/scripts/applications/impute/imputeGaussMCMC.nogradient.dml b/src/test/scripts/applications/impute/imputeGaussMCMC.nogradient.dml
index f30698c..00210c5 100644
--- a/src/test/scripts/applications/impute/imputeGaussMCMC.nogradient.dml
+++ b/src/test/scripts/applications/impute/imputeGaussMCMC.nogradient.dml
@@ -1,453 +1,453 @@
-#-------------------------------------------------------------
-#
-# 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.
-#
-#-------------------------------------------------------------
-
-# Implements the MCMC algorithm for imputation of missing data into a time-series of "reports".
-# Each report is a fixed-size vector of attribute values; reports come out each year/quarter/month ("term").
-# Hard linear equality constraints restrict values in/across the reports, e.g. total cost = sum of all costs.
-# Soft linear regression constraints define dependencies between values in/across the reports.
-# Linear regression parameters are unknown and sampled together with the missing values in the reports.
-#
-# INPUT 1: Initial reports matrix [1 : num_attrs, 1 : num_terms] with missing values usually set to zero,
-# but it MUST BE CONSISTENT with hard constraints! Set some missing values to nonzero if needed.
-# There are "num_terms" reports in the matrix, each having "num_attrs" attribute values.
-#
-# INPUT 2: Sparse matrix [1 : (num_terms * num_attrs), 1 : num_frees] that defines a linear map from
-# "free" variables to the reports matrix. A tensor of size (num_terms * num_attrs * num_frees)
-# where the reports matrix is stretched into a column-vector [1 : (num_terms * num_attrs)].
-# Term = t, attribute = i --> index = (t-1) * num_attrs + i
-#
-# INPUT 3: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : (num_terms * num_attrs)] that defines
-# a linear map from the stretched matrix of reports to the stretched matrix of regression factors.
-#
-# INPUT 4: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines "default regression factors"
-# (if nonzero) to be added to the regression factors before they are multiplied by parameters.
-#
-# INPUT 5: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : num_params] that defines a linear map
-# from the vector of parameters to the stretched matrix of regression factors.
-#
-# INPUT 6: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines default regression coefficients
-# (if nonzero) to be added to the parameters (if any) before being multiplied by regression factors.
-#
-# INPUT 7: A vector [1 : num_reg_eqs, 1] of scale multipliers, one per regression
-#
-# INPUT 8 : Number of factors in a regression equation, including the estimated value
-# INPUT 9 : Maximum number of burn-in full iterations (that sample each variable and each parameter once)
-# BUT the actual number of burn-in iterations may be smaller if "free fall" ends sooner
-# INPUT 10: Maximum number of observed full iterations (that sample each variable and each parameter once)
-#
-# INPUT 11: Output file name and path for the average MCMC reports table
-# INPUT 12: Output file for debugging (currently: the average parameters vector)
-#
-# Example:
-# hadoop jar SystemML.jar -f test/scripts/applications/impute/imputeGaussMCMC.dml -exec singlenode -args
-# test/scripts/applications/impute/initial_reports
-# test/scripts/applications/impute/CReps
-# test/scripts/applications/impute/RegresValueMap
-# test/scripts/applications/impute/RegresFactorDefault
-# test/scripts/applications/impute/RegresParamMap
-# test/scripts/applications/impute/RegresCoeffDefault
-# test/scripts/applications/impute/RegresScaleMult
-# 4 1000 100
-# test/scripts/applications/impute/output_reports
-# test/scripts/applications/impute/debug_info
-
-
-print ("START ImputeGaussMCMC");
-print ("Reading the input files...");
-
-initial_reports = read ($1);
-CReps = read ($2);
-
-num_terms = ncol (initial_reports); # Number of periods for which reports are produced
-num_attrs = nrow (initial_reports); # Number of attribute values per each term report
-num_frees = ncol (CReps); # Number of free variables used to describe all consistent reports
-
-dReps_size = num_terms * num_attrs;
-dReps = matrix (initial_reports, rows = dReps_size, cols = 1, byrow = FALSE);
-
-# We assume that all report-series consistent with hard constraints form an affine set:
-# reports = CReps %*% freeVars + dReps
-# Here "freeVars" is a vector of "free variables" (degrees of freedom), "CReps" is a linear mapping,
-# and "dReps" are the default values for the reports that correspond to the initial reports matrix.
-
-RegresValueMap = read ($3);
-RegresFactorDefault = read ($4);
-RegresParamMap = read ($5);
-RegresCoeffDefault = read ($6);
-RegresScaleMult = read ($7);
-
-num_factors = $8; # Number of factors in each regression equation, including the estimated value
-num_reg_eqs = nrow (RegresParamMap) / num_factors; # Number of regression equations
-num_params = ncol (RegresParamMap); # Number of parameters used in all regressions
-max_num_burnin_iterations = $9;
-max_num_observed_iterations = $10;
-
-num_opt_iter = 20;
-print ("Performing initial optimization (" + num_opt_iter + " alternating CG steps)...");
-
-freeVars = matrix (0.0, rows = num_frees, cols = 1);
-params = matrix (1.0, rows = num_params, cols = 1);
-reports = CReps %*% freeVars + dReps;
-
-regresValues = RegresValueMap %*% reports + RegresFactorDefault;
-regresParams = RegresParamMap %*% params + RegresCoeffDefault;
-bilinear_vector = regresValues * regresParams;
-
-### DELETE: bilinear_form = matricize (bilinear_vector, num_factors);
-bilinear_form = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
-bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
-
-opt_iter = 1;
-is_step_params = 1;
-is_opt_converged = 0;
-
-print ("Before optimization: Initial bilinear form value = " + bilinear_form_value);
-
-
-while (is_opt_converged == 0)
-{
- deg = is_step_params * num_params + (1 - is_step_params) * num_frees;
-
- # Compute gradient
-
- gradient = matrix (0.0, rows = deg, cols = 1);
- for (i in 1:deg)
- {
- if (is_step_params == 1) {
- bilinear_vector = regresValues * RegresParamMap [, i];
- } else {
- bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
- }
- ### DELETE: bilinear_updater = matricize (bilinear_vector, num_factors);
- bilinear_updater = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
-
- q_minus_1 = sum (RegresScaleMult * rowSums (bilinear_form - bilinear_updater) * rowSums (bilinear_form - bilinear_updater));
- q_plus_1 = sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater) * rowSums (bilinear_form + bilinear_updater));
- gradient [i, 1] = 0.5 * (q_plus_1 - q_minus_1) + gradient [i, 1];
- }
-
- # Make a few conjugate gradient steps
-
- shift_vector = matrix (0.0, rows = deg, cols = 1);
- residual = gradient;
- p = - residual;
- norm_r2 = sum (residual * residual);
-
- for (j in 1:3)
- {
- q = matrix (0.0, rows = deg, cols = 1);
- for (i in 1:deg) # Want: q = A %*% p;
- {
- if (is_step_params == 1) {
- bilinear_vector = regresValues * RegresParamMap [, i];
- } else {
- bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
- }
- ### DELETE: bilinear_updater_1 = matricize (bilinear_vector, num_factors);
- bilinear_updater_1 = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
-
- if (is_step_params == 1) {
- bilinear_vector = regresValues * (RegresParamMap %*% p);
- } else {
- bilinear_vector = (RegresValueMap %*% CReps %*% p) * regresParams;
- }
- ### DELETE: bilinear_updater_p = matricize (bilinear_vector, num_factors);
- bilinear_updater_p = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
-
- quadratic_plus_1 =
- sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_1) * rowSums (bilinear_form + bilinear_updater_1));
- quadratic_plus_p =
- sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_p) * rowSums (bilinear_form + bilinear_updater_p));
- quadratic_plus_both =
- sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p) * rowSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p));
- q [i, 1] = (quadratic_plus_both - quadratic_plus_1 - quadratic_plus_p + bilinear_form_value) + q [i, 1];
- }
-
- alpha = norm_r2 / castAsScalar (t(p) %*% q);
- shift_vector = shift_vector + alpha * p;
- old_norm_r2 = norm_r2;
- residual = residual + alpha * q;
- norm_r2 = sum (residual * residual);
- p = - residual + (norm_r2 / old_norm_r2) * p;
- }
-
- if (is_step_params == 1) {
- params = params + shift_vector;
- regresParams = RegresParamMap %*% params + RegresCoeffDefault;
- } else {
- freeVars = freeVars + shift_vector;
- reports = CReps %*% freeVars + dReps;
- regresValues = RegresValueMap %*% reports + RegresFactorDefault;
- }
-
- # Update the bilinear form and check convergence
-
- if (is_step_params == 1) {
- old_bilinear_form_value = bilinear_form_value;
- }
- bilinear_vector = regresValues * regresParams;
-
- ### DELETE: bilinear_form = matricize (bilinear_vector, num_factors);
- bilinear_form = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
- bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
-
- if (is_step_params == 1) {
- print ("Optimization step " + opt_iter + " (params) : bilinear form value = " + bilinear_form_value);
- } else {
- print ("Optimization step " + opt_iter + " (reports): bilinear form value = " + bilinear_form_value);
- }
-
- is_step_params = 1 - is_step_params;
- opt_iter = opt_iter + 1;
-
- if (is_step_params == 1 & opt_iter > num_opt_iter) {
- is_opt_converged = 1;
- }
-}
-
-print ("Performing MCMC initialization...");
-
-max_num_iter = max_num_burnin_iterations + max_num_observed_iterations;
-dim_sample = num_frees + num_params;
-sample_ones = matrix (1.0, rows = dim_sample, cols = 1);
-
-# Generate a random permutation matrix for the sampling order of freeVars and params
-
-SampleOrder = diag (sample_ones);
-num_swaps = 10 * dim_sample;
-rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
-left_swap = round (0.5 + dim_sample * rnd);
-rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
-right_swap = round (0.5 + dim_sample * rnd);
-for (swap_i in 1:num_swaps) {
- l = castAsScalar (left_swap [swap_i, 1]);
- r = castAsScalar (right_swap [swap_i, 1]);
- if (l != r) {
- tmp_row = SampleOrder [l, ];
- SampleOrder [l, ] = SampleOrder [r, ];
- SampleOrder [r, ] = tmp_row;
- }
-}
-
-pi = 3.1415926535897932384626433832795;
-zero = matrix (0.0, rows = 1, cols = 1);
-
-isVar = colSums (SampleOrder [1 : num_frees, ]);
-sum_of_observed_reports = matrix (0.0, rows = num_attrs, cols = num_terms);
-sum_of_observed_params = matrix (0.0, rows = num_params, cols = 1);
-num_of_observed_reports = 0;
-sum_of_observed_losses = 0.0;
-is_observed = 0;
-
-is_calculating_loss_change = 0;
-is_monitoring_loss_change = 0;
-avg_prob_of_loss_increase = 0;
-update_factor_for_avg_loss_change = 0.02;
-avg_loss_change = -50.0 * update_factor_for_avg_loss_change;
-old_bilinear_form_value = bilinear_form_value;
-
-# Starting MCMC iterations
-
-iter = 0;
-
-while ((iter < max_num_iter) & (num_of_observed_reports < max_num_observed_iterations))
-{
- iter = iter + 1;
-
- # Initialize (bi-)linear forms
-
- regresValues = RegresValueMap %*% reports + RegresFactorDefault;
- regresParams = RegresParamMap %*% params + RegresCoeffDefault;
- bilinear_form_vector = regresValues * regresParams;
-
- ### DELETE: bilinear_form = matricize (bilinear_form_vector, num_factors);
- bilinear_form = matrix (bilinear_form_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
- bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
-
- if (bilinear_form_value > old_bilinear_form_value) {
- avg_prob_of_loss_increase = avg_prob_of_loss_increase * (1 - update_factor_for_avg_loss_change) + 1 * update_factor_for_avg_loss_change;
- } else {
- avg_prob_of_loss_increase = avg_prob_of_loss_increase * (1 - update_factor_for_avg_loss_change);
- }
- if (is_calculating_loss_change == 0 & avg_prob_of_loss_increase > 0.4) {
- is_calculating_loss_change = 1;
- }
- if (is_monitoring_loss_change == 0 & avg_prob_of_loss_increase > 0.5) {
- is_calculating_loss_change = 1;
- is_monitoring_loss_change = 1;
- print ("Monitoring the average loss change is ON. ");
- }
- if (is_calculating_loss_change == 1) {
- avg_loss_change = avg_loss_change * (1 - update_factor_for_avg_loss_change)
- + (bilinear_form_value - old_bilinear_form_value) * update_factor_for_avg_loss_change;
- }
- if (is_observed == 0 & ((is_monitoring_loss_change == 1 & avg_loss_change > 0) | iter > max_num_burnin_iterations)) {
- print ("Burn-in ENDS, observation STARTS. ");
- is_observed = 1;
- }
-
- old_bilinear_form_value = bilinear_form_value;
-
- bilinear_form_value_to_print = bilinear_form_value;
- if (bilinear_form_value < 100000) {
- bilinear_form_value_to_print = round (10000 * bilinear_form_value) / 10000;
- } else {
- if (bilinear_form_value < 1000000000) {
- bilinear_form_value_to_print = round (bilinear_form_value);
- }}
-
- if (is_monitoring_loss_change == 0) {
- print ("MCMC iteration " + iter + ": Prob [loss_increase] = " + (round (10000 * avg_prob_of_loss_increase) / 10000)
- + ", bilinear form value = " + bilinear_form_value_to_print);
- } else {
- print ("MCMC iteration " + iter + ": Prob [loss_increase] = " + (round (10000 * avg_prob_of_loss_increase) / 10000)
- + ", bilinear form value = " + bilinear_form_value_to_print + ", avg_loss_change = " + (round (10000 * avg_loss_change) / 10000));
- }
-
- # Create a normally distributed random sample
-
- dim_half_sample = castAsScalar (round (dim_sample / 2 + 0.1 + zero));
- rnd1 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
- rnd2 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
- rnd_normal_1 = sqrt (- 2.0 * log (rnd1)) * sin (2 * pi * rnd2);
- rnd_normal_2 = sqrt (- 2.0 * log (rnd1)) * cos (2 * pi * rnd2);
- rnd_normal = matrix (0.0, rows = dim_sample, cols = 1);
- rnd_normal [1 : dim_half_sample, ] = rnd_normal_1;
- rnd_normal [(dim_sample - dim_half_sample + 1) : dim_sample, ] = rnd_normal_2;
-
- # Initialize updaters
-
- freeVars_updater = freeVars * 0.0;
- params_updater = params * 0.0;
- regresValues_updater = regresValues * 0.0;
- regresParams_updater = regresParams * 0.0;
- bilinear_updater_vector = bilinear_form_vector * 0.0;
-
- # Perform the sampling
-
- for (idx in 1:dim_sample)
- {
- # Generate the sample unit-vector and updaters
-
- if (castAsScalar (isVar [1, idx]) > 0.5) {
- freeVars_updater = SampleOrder [1 : num_frees, idx];
- regresValues_updater = RegresValueMap %*% CReps %*% freeVars_updater;
- bilinear_updater_vector = regresValues_updater * regresParams;
- } else {
- params_updater = SampleOrder [(num_frees + 1) : dim_sample, idx];
- regresParams_updater = RegresParamMap %*% params_updater;
- bilinear_updater_vector = regresValues * regresParams_updater;
- }
- ### DELETE: bilinear_updater = matricize (bilinear_updater_vector, num_factors);
- bilinear_updater = matrix (bilinear_updater_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
-
- # Compute the quadratic by three shift-points: -1, 0, +1
-
- bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
- q_minus_1 = sum (RegresScaleMult * rowSums (bilinear_form - bilinear_updater) * rowSums (bilinear_form - bilinear_updater));
- q_plus_1 = sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater) * rowSums (bilinear_form + bilinear_updater));
- coeff_b = (q_plus_1 - q_minus_1) / 2.0;
- coeff_a = (q_plus_1 + q_minus_1) / 2.0 - bilinear_form_value;
-
- # Find the mean and the sigma for f(x) ~ exp (- (ax^2 + bx + c)),
- # then compute the shift to get the new sample
-
- mean_shift = - coeff_b / (2.0 * coeff_a);
- sigma_shift = 1.0 / sqrt (2.0 * coeff_a);
- shift = mean_shift + sigma_shift * castAsScalar (rnd_normal [idx, 1]);
-
-# BEGIN DEBUG INSERT
-# mmm = 1;
-# if (castAsScalar (isVar [1, idx]) > 0.5 & # IT IS A FREE VARIABLE, NOT A PARAMETER
-# castAsScalar (freeVars_updater [mmm, 1]) > 0) # IT IS mmm-TH FREE VARIABLE
-# {
-# # print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a + ", coeff_b = " + coeff_b);
-# print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", mean_shift = " + mean_shift + ", sigma_shift = " + sigma_shift + ", shift = " + shift);
-# }
-# if (castAsScalar (isVar [1, idx]) <= 0.5 & # IT IS A PARAMETER, NOT A FREE VARIABLE
-# castAsScalar (params_updater [mmm, 1]) > 0) # IT IS mmm-TH PARAMETER
-# {
-# # print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a + ", coeff_b = " + coeff_b);
-# print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", mean_shift = " + mean_shift + ", sigma_shift = " + sigma_shift + ", shift = " + shift);
-# }
-# END DEBUG INSERT
-
- # Perform the updates
-
- bilinear_form = bilinear_form + shift * bilinear_updater;
- if (castAsScalar (isVar [1, idx]) > 0.5) {
- freeVars = freeVars + shift * freeVars_updater;
- regresValues = regresValues + shift * regresValues_updater;
- } else {
- params = params + shift * params_updater;
- regresParams = regresParams + shift * regresParams_updater;
- }
- }
-
- # Update / adjust the reports and the parameters
-
- reports = CReps %*% freeVars + dReps;
- ### DELETE: reports_matrix = matricize (reports, num_attrs);
- reports_matrix = matrix (reports, rows = num_attrs, cols = num_terms, byrow = FALSE);
-
- # Make an observation of the reports and/or the parameters
-
- if (is_observed > 0)
- {
- sum_of_observed_reports = sum_of_observed_reports + reports_matrix;
- num_of_observed_reports = num_of_observed_reports + 1;
-
- sum_of_observed_params = sum_of_observed_params + params;
- sum_of_observed_losses = sum_of_observed_losses + bilinear_form_value;
- }
-
-# v1 =castAsScalar(round(10000*reports[1 + (num_terms - 1) * num_attrs, 1])/10000);
-# v2 =castAsScalar(round(10000*reports[2 + (num_terms - 1) * num_attrs, 1])/10000);
-# v3 =castAsScalar(round(10000*reports[3 + (num_terms - 1) * num_attrs, 1])/10000);
-# v4 =castAsScalar(round(10000*reports[4 + (num_terms - 1) * num_attrs, 1])/10000);
-# w1 =castAsScalar(round(10000*reports_matrix[ 1,num_terms])/10000);
-# w2 =castAsScalar(round(10000*reports_matrix[ 2,num_terms])/10000);
-# w3 =castAsScalar(round(10000*reports_matrix[ 3,num_terms])/10000);
-# w4 =castAsScalar(round(10000*reports_matrix[ 4,num_terms])/10000);
-
-# v5 =castAsScalar(round(reports_matrix[ 5,num_terms]));
-# v8 =castAsScalar(round(reports_matrix[ 8,num_terms]));
-# v9 =castAsScalar(round(reports_matrix[ 9,num_terms]));
-# v10=castAsScalar(round(reports_matrix[10,num_terms]));
-# v16=castAsScalar(round(reports_matrix[16,num_terms]));
-# v19=castAsScalar(round(reports_matrix[19,num_terms]));
-
-#print (" Sample = 1:" + v1 + ", 2:" + v2 + ", 3:" + v3 + ", 4:" + v4);
-## + ", 5:" + v5 + ", 8:" + v8 + ", 9:" + v9 + ", 10:" + v10 + ", 16:" + v16 + ", 19:" + v19);
-#print (" Sample = 1:" + w1 + ", 2:" + w2 + ", 3:" + w3 + ", 4:" + w4);
-## + ", 5:" + w5 + ", 8:" + w8 + ", 9:" + w9 + ", 10:" + w10 + ", 16:" + w16 + ", 19:" + w19);
-
-}
-
-print ("Average observed loss = " + (sum_of_observed_losses / num_of_observed_reports));
-print ("Writing out the results...");
-
-avg_reports_matrix = sum_of_observed_reports / num_of_observed_reports;
-avg_params = sum_of_observed_params / num_of_observed_reports;
-write (avg_reports_matrix, $11, format="text");
-write (avg_params, $12, format="text");
-
-print ("END ImputeGaussMCMC");
+#-------------------------------------------------------------
+#
+# 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.
+#
+#-------------------------------------------------------------
+
+# Implements the MCMC algorithm for imputation of missing data into a time-series of "reports".
+# Each report is a fixed-size vector of attribute values; reports come out each year/quarter/month ("term").
+# Hard linear equality constraints restrict values in/across the reports, e.g. total cost = sum of all costs.
+# Soft linear regression constraints define dependencies between values in/across the reports.
+# Linear regression parameters are unknown and sampled together with the missing values in the reports.
+#
+# INPUT 1: Initial reports matrix [1 : num_attrs, 1 : num_terms] with missing values usually set to zero,
+# but it MUST BE CONSISTENT with hard constraints! Set some missing values to nonzero if needed.
+# There are "num_terms" reports in the matrix, each having "num_attrs" attribute values.
+#
+# INPUT 2: Sparse matrix [1 : (num_terms * num_attrs), 1 : num_frees] that defines a linear map from
+# "free" variables to the reports matrix. A tensor of size (num_terms * num_attrs * num_frees)
+# where the reports matrix is stretched into a column-vector [1 : (num_terms * num_attrs)].
+# Term = t, attribute = i --> index = (t-1) * num_attrs + i
+#
+# INPUT 3: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : (num_terms * num_attrs)] that defines
+# a linear map from the stretched matrix of reports to the stretched matrix of regression factors.
+#
+# INPUT 4: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines "default regression factors"
+# (if nonzero) to be added to the regression factors before they are multiplied by parameters.
+#
+# INPUT 5: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : num_params] that defines a linear map
+# from the vector of parameters to the stretched matrix of regression factors.
+#
+# INPUT 6: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines default regression coefficients
+# (if nonzero) to be added to the parameters (if any) before being multiplied by regression factors.
+#
+# INPUT 7: A vector [1 : num_reg_eqs, 1] of scale multipliers, one per regression
+#
+# INPUT 8 : Number of factors in a regression equation, including the estimated value
+# INPUT 9 : Maximum number of burn-in full iterations (that sample each variable and each parameter once)
+# BUT the actual number of burn-in iterations may be smaller if "free fall" ends sooner
+# INPUT 10: Maximum number of observed full iterations (that sample each variable and each parameter once)
+#
+# INPUT 11: Output file name and path for the average MCMC reports table
+# INPUT 12: Output file for debugging (currently: the average parameters vector)
+#
+# Example:
+# hadoop jar SystemML.jar -f test/scripts/applications/impute/imputeGaussMCMC.dml -exec singlenode -args
+# test/scripts/applications/impute/initial_reports
+# test/scripts/applications/impute/CReps
+# test/scripts/applications/impute/RegresValueMap
+# test/scripts/applications/impute/RegresFactorDefault
+# test/scripts/applications/impute/RegresParamMap
+# test/scripts/applications/impute/RegresCoeffDefault
+# test/scripts/applications/impute/RegresScaleMult
+# 4 1000 100
+# test/scripts/applications/impute/output_reports
+# test/scripts/applications/impute/debug_info
+
+
+print ("START ImputeGaussMCMC");
+print ("Reading the input files...");
+
+initial_reports = read ($1);
+CReps = read ($2);
+
+num_terms = ncol (initial_reports); # Number of periods for which reports are produced
+num_attrs = nrow (initial_reports); # Number of attribute values per each term report
+num_frees = ncol (CReps); # Number of free variables used to describe all consistent reports
+
+dReps_size = num_terms * num_attrs;
+dReps = matrix (initial_reports, rows = dReps_size, cols = 1, byrow = FALSE);
+
+# We assume that all report-series consistent with hard constraints form an affine set:
+# reports = CReps %*% freeVars + dReps
+# Here "freeVars" is a vector of "free variables" (degrees of freedom), "CReps" is a linear mapping,
+# and "dReps" are the default values for the reports that correspond to the initial reports matrix.
+
+RegresValueMap = read ($3);
+RegresFactorDefault = read ($4);
+RegresParamMap = read ($5);
+RegresCoeffDefault = read ($6);
+RegresScaleMult = read ($7);
+
+num_factors = $8; # Number of factors in each regression equation, including the estimated value
+num_reg_eqs = nrow (RegresParamMap) / num_factors; # Number of regression equations
+num_params = ncol (RegresParamMap); # Number of parameters used in all regressions
+max_num_burnin_iterations = $9;
+max_num_observed_iterations = $10;
+
+num_opt_iter = 20;
+print ("Performing initial optimization (" + num_opt_iter + " alternating CG steps)...");
+
+freeVars = matrix (0.0, rows = num_frees, cols = 1);
+params = matrix (1.0, rows = num_params, cols = 1);
+reports = CReps %*% freeVars + dReps;
+
+regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+bilinear_vector = regresValues * regresParams;
+
+### DELETE: bilinear_form = matricize (bilinear_vector, num_factors);
+bilinear_form = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
+
+opt_iter = 1;
+is_step_params = 1;
+is_opt_converged = 0;
+
+print ("Before optimization: Initial bilinear form value = " + bilinear_form_value);
+
+
+while (is_opt_converged == 0)
+{
+ deg = is_step_params * num_params + (1 - is_step_params) * num_frees;
+
+ # Compute gradient
+
+ gradient = matrix (0.0, rows = deg, cols = 1);
+ for (i in 1:deg)
+ {
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * RegresParamMap [, i];
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
+ }
+ ### DELETE: bilinear_updater = matricize (bilinear_vector, num_factors);
+ bilinear_updater = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+
+ q_minus_1 = sum (RegresScaleMult * rowSums (bilinear_form - bilinear_updater) * rowSums (bilinear_form - bilinear_updater));
+ q_plus_1 = sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater) * rowSums (bilinear_form + bilinear_updater));
+ gradient [i, 1] = 0.5 * (q_plus_1 - q_minus_1) + gradient [i, 1];
+ }
+
+ # Make a few conjugate gradient steps
+
+ shift_vector = matrix (0.0, rows = deg, cols = 1);
+ residual = gradient;
+ p = - residual;
+ norm_r2 = sum (residual * residual);
+
+ for (j in 1:3)
+ {
+ q = matrix (0.0, rows = deg, cols = 1);
+ for (i in 1:deg) # Want: q = A %*% p;
+ {
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * RegresParamMap [, i];
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
+ }
+ ### DELETE: bilinear_updater_1 = matricize (bilinear_vector, num_factors);
+ bilinear_updater_1 = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * (RegresParamMap %*% p);
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps %*% p) * regresParams;
+ }
+ ### DELETE: bilinear_updater_p = matricize (bilinear_vector, num_factors);
+ bilinear_updater_p = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+
+ quadratic_plus_1 =
+ sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_1) * rowSums (bilinear_form + bilinear_updater_1));
+ quadratic_plus_p =
+ sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_p) * rowSums (bilinear_form + bilinear_updater_p));
+ quadratic_plus_both =
+ sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p) * rowSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p));
+ q [i, 1] = (quadratic_plus_both - quadratic_plus_1 - quadratic_plus_p + bilinear_form_value) + q [i, 1];
+ }
+
+ alpha = norm_r2 / castAsScalar (t(p) %*% q);
+ shift_vector = shift_vector + alpha * p;
+ old_norm_r2 = norm_r2;
+ residual = residual + alpha * q;
+ norm_r2 = sum (residual * residual);
+ p = - residual + (norm_r2 / old_norm_r2) * p;
+ }
+
+ if (is_step_params == 1) {
+ params = params + shift_vector;
+ regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+ } else {
+ freeVars = freeVars + shift_vector;
+ reports = CReps %*% freeVars + dReps;
+ regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+ }
+
+ # Update the bilinear form and check convergence
+
+ if (is_step_params == 1) {
+ old_bilinear_form_value = bilinear_form_value;
+ }
+ bilinear_vector = regresValues * regresParams;
+
+ ### DELETE: bilinear_form = matricize (bilinear_vector, num_factors);
+ bilinear_form = matrix (bilinear_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+ bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
+
+ if (is_step_params == 1) {
+ print ("Optimization step " + opt_iter + " (params) : bilinear form value = " + bilinear_form_value);
+ } else {
+ print ("Optimization step " + opt_iter + " (reports): bilinear form value = " + bilinear_form_value);
+ }
+
+ is_step_params = 1 - is_step_params;
+ opt_iter = opt_iter + 1;
+
+ if (is_step_params == 1 & opt_iter > num_opt_iter) {
+ is_opt_converged = 1;
+ }
+}
+
+print ("Performing MCMC initialization...");
+
+max_num_iter = max_num_burnin_iterations + max_num_observed_iterations;
+dim_sample = num_frees + num_params;
+sample_ones = matrix (1.0, rows = dim_sample, cols = 1);
+
+# Generate a random permutation matrix for the sampling order of freeVars and params
+
+SampleOrder = diag (sample_ones);
+num_swaps = 10 * dim_sample;
+rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
+left_swap = round (0.5 + dim_sample * rnd);
+rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
+right_swap = round (0.5 + dim_sample * rnd);
+for (swap_i in 1:num_swaps) {
+ l = castAsScalar (left_swap [swap_i, 1]);
+ r = castAsScalar (right_swap [swap_i, 1]);
+ if (l != r) {
+ tmp_row = SampleOrder [l, ];
+ SampleOrder [l, ] = SampleOrder [r, ];
+ SampleOrder [r, ] = tmp_row;
+ }
+}
+
+pi = 3.1415926535897932384626433832795;
+zero = matrix (0.0, rows = 1, cols = 1);
+
+isVar = colSums (SampleOrder [1 : num_frees, ]);
+sum_of_observed_reports = matrix (0.0, rows = num_attrs, cols = num_terms);
+sum_of_observed_params = matrix (0.0, rows = num_params, cols = 1);
+num_of_observed_reports = 0;
+sum_of_observed_losses = 0.0;
+is_observed = 0;
+
+is_calculating_loss_change = 0;
+is_monitoring_loss_change = 0;
+avg_prob_of_loss_increase = 0;
+update_factor_for_avg_loss_change = 0.02;
+avg_loss_change = -50.0 * update_factor_for_avg_loss_change;
+old_bilinear_form_value = bilinear_form_value;
+
+# Starting MCMC iterations
+
+iter = 0;
+
+while ((iter < max_num_iter) & (num_of_observed_reports < max_num_observed_iterations))
+{
+ iter = iter + 1;
+
+ # Initialize (bi-)linear forms
+
+ regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+ regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+ bilinear_form_vector = regresValues * regresParams;
+
+ ### DELETE: bilinear_form = matricize (bilinear_form_vector, num_factors);
+ bilinear_form = matrix (bilinear_form_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+ bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
+
+ if (bilinear_form_value > old_bilinear_form_value) {
+ avg_prob_of_loss_increase = avg_prob_of_loss_increase * (1 - update_factor_for_avg_loss_change) + 1 * update_factor_for_avg_loss_change;
+ } else {
+ avg_prob_of_loss_increase = avg_prob_of_loss_increase * (1 - update_factor_for_avg_loss_change);
+ }
+ if (is_calculating_loss_change == 0 & avg_prob_of_loss_increase > 0.4) {
+ is_calculating_loss_change = 1;
+ }
+ if (is_monitoring_loss_change == 0 & avg_prob_of_loss_increase > 0.5) {
+ is_calculating_loss_change = 1;
+ is_monitoring_loss_change = 1;
+ print ("Monitoring the average loss change is ON. ");
+ }
+ if (is_calculating_loss_change == 1) {
+ avg_loss_change = avg_loss_change * (1 - update_factor_for_avg_loss_change)
+ + (bilinear_form_value - old_bilinear_form_value) * update_factor_for_avg_loss_change;
+ }
+ if (is_observed == 0 & ((is_monitoring_loss_change == 1 & avg_loss_change > 0) | iter > max_num_burnin_iterations)) {
+ print ("Burn-in ENDS, observation STARTS. ");
+ is_observed = 1;
+ }
+
+ old_bilinear_form_value = bilinear_form_value;
+
+ bilinear_form_value_to_print = bilinear_form_value;
+ if (bilinear_form_value < 100000) {
+ bilinear_form_value_to_print = round (10000 * bilinear_form_value) / 10000;
+ } else {
+ if (bilinear_form_value < 1000000000) {
+ bilinear_form_value_to_print = round (bilinear_form_value);
+ }}
+
+ if (is_monitoring_loss_change == 0) {
+ print ("MCMC iteration " + iter + ": Prob [loss_increase] = " + (round (10000 * avg_prob_of_loss_increase) / 10000)
+ + ", bilinear form value = " + bilinear_form_value_to_print);
+ } else {
+ print ("MCMC iteration " + iter + ": Prob [loss_increase] = " + (round (10000 * avg_prob_of_loss_increase) / 10000)
+ + ", bilinear form value = " + bilinear_form_value_to_print + ", avg_loss_change = " + (round (10000 * avg_loss_change) / 10000));
+ }
+
+ # Create a normally distributed random sample
+
+ dim_half_sample = castAsScalar (round (dim_sample / 2 + 0.1 + zero));
+ rnd1 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
+ rnd2 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
+ rnd_normal_1 = sqrt (- 2.0 * log (rnd1)) * sin (2 * pi * rnd2);
+ rnd_normal_2 = sqrt (- 2.0 * log (rnd1)) * cos (2 * pi * rnd2);
+ rnd_normal = matrix (0.0, rows = dim_sample, cols = 1);
+ rnd_normal [1 : dim_half_sample, ] = rnd_normal_1;
+ rnd_normal [(dim_sample - dim_half_sample + 1) : dim_sample, ] = rnd_normal_2;
+
+ # Initialize updaters
+
+ freeVars_updater = freeVars * 0.0;
+ params_updater = params * 0.0;
+ regresValues_updater = regresValues * 0.0;
+ regresParams_updater = regresParams * 0.0;
+ bilinear_updater_vector = bilinear_form_vector * 0.0;
+
+ # Perform the sampling
+
+ for (idx in 1:dim_sample)
+ {
+ # Generate the sample unit-vector and updaters
+
+ if (castAsScalar (isVar [1, idx]) > 0.5) {
+ freeVars_updater = SampleOrder [1 : num_frees, idx];
+ regresValues_updater = RegresValueMap %*% CReps %*% freeVars_updater;
+ bilinear_updater_vector = regresValues_updater * regresParams;
+ } else {
+ params_updater = SampleOrder [(num_frees + 1) : dim_sample, idx];
+ regresParams_updater = RegresParamMap %*% params_updater;
+ bilinear_updater_vector = regresValues * regresParams_updater;
+ }
+ ### DELETE: bilinear_updater = matricize (bilinear_updater_vector, num_factors);
+ bilinear_updater = matrix (bilinear_updater_vector, rows = num_reg_eqs, cols = num_factors, byrow = TRUE);
+
+ # Compute the quadratic by three shift-points: -1, 0, +1
+
+ bilinear_form_value = sum (RegresScaleMult * rowSums (bilinear_form) * rowSums (bilinear_form));
+ q_minus_1 = sum (RegresScaleMult * rowSums (bilinear_form - bilinear_updater) * rowSums (bilinear_form - bilinear_updater));
+ q_plus_1 = sum (RegresScaleMult * rowSums (bilinear_form + bilinear_updater) * rowSums (bilinear_form + bilinear_updater));
+ coeff_b = (q_plus_1 - q_minus_1) / 2.0;
+ coeff_a = (q_plus_1 + q_minus_1) / 2.0 - bilinear_form_value;
+
+ # Find the mean and the sigma for f(x) ~ exp (- (ax^2 + bx + c)),
+ # then compute the shift to get the new sample
+
+ mean_shift = - coeff_b / (2.0 * coeff_a);
+ sigma_shift = 1.0 / sqrt (2.0 * coeff_a);
+ shift = mean_shift + sigma_shift * castAsScalar (rnd_normal [idx, 1]);
+
+# BEGIN DEBUG INSERT
+# mmm = 1;
+# if (castAsScalar (isVar [1, idx]) > 0.5 & # IT IS A FREE VARIABLE, NOT A PARAMETER
+# castAsScalar (freeVars_updater [mmm, 1]) > 0) # IT IS mmm-TH FREE VARIABLE
+# {
+# # print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a + ", coeff_b = " + coeff_b);
+# print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", mean_shift = " + mean_shift + ", sigma_shift = " + sigma_shift + ", shift = " + shift);
+# }
+# if (castAsScalar (isVar [1, idx]) <= 0.5 & # IT IS A PARAMETER, NOT A FREE VARIABLE
+# castAsScalar (params_updater [mmm, 1]) > 0) # IT IS mmm-TH PARAMETER
+# {
+# # print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a + ", coeff_b = " + coeff_b);
+# print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", mean_shift = " + mean_shift + ", sigma_shift = " + sigma_shift + ", shift = " + shift);
+# }
+# END DEBUG INSERT
+
+ # Perform the updates
+
+ bilinear_form = bilinear_form + shift * bilinear_updater;
+ if (castAsScalar (isVar [1, idx]) > 0.5) {
+ freeVars = freeVars + shift * freeVars_updater;
+ regresValues = regresValues + shift * regresValues_updater;
+ } else {
+ params = params + shift * params_updater;
+ regresParams = regresParams + shift * regresParams_updater;
+ }
+ }
+
+ # Update / adjust the reports and the parameters
+
+ reports = CReps %*% freeVars + dReps;
+ ### DELETE: reports_matrix = matricize (reports, num_attrs);
+ reports_matrix = matrix (reports, rows = num_attrs, cols = num_terms, byrow = FALSE);
+
+ # Make an observation of the reports and/or the parameters
+
+ if (is_observed > 0)
+ {
+ sum_of_observed_reports = sum_of_observed_reports + reports_matrix;
+ num_of_observed_reports = num_of_observed_reports + 1;
+
+ sum_of_observed_params = sum_of_observed_params + params;
+ sum_of_observed_losses = sum_of_observed_losses + bilinear_form_value;
+ }
+
+# v1 =castAsScalar(round(10000*reports[1 + (num_terms - 1) * num_attrs, 1])/10000);
+# v2 =castAsScalar(round(10000*reports[2 + (num_terms - 1) * num_attrs, 1])/10000);
+# v3 =castAsScalar(round(10000*reports[3 + (num_terms - 1) * num_attrs, 1])/10000);
+# v4 =castAsScalar(round(10000*reports[4 + (num_terms - 1) * num_attrs, 1])/10000);
+# w1 =castAsScalar(round(10000*reports_matrix[ 1,num_terms])/10000);
+# w2 =castAsScalar(round(10000*reports_matrix[ 2,num_terms])/10000);
+# w3 =castAsScalar(round(10000*reports_matrix[ 3,num_terms])/10000);
+# w4 =castAsScalar(round(10000*reports_matrix[ 4,num_terms])/10000);
+
+# v5 =castAsScalar(round(reports_matrix[ 5,num_terms]));
+# v8 =castAsScalar(round(reports_matrix[ 8,num_terms]));
+# v9 =castAsScalar(round(reports_matrix[ 9,num_terms]));
+# v10=castAsScalar(round(reports_matrix[10,num_terms]));
+# v16=castAsScalar(round(reports_matrix[16,num_terms]));
+# v19=castAsScalar(round(reports_matrix[19,num_terms]));
+
+#print (" Sample = 1:" + v1 + ", 2:" + v2 + ", 3:" + v3 + ", 4:" + v4);
+## + ", 5:" + v5 + ", 8:" + v8 + ", 9:" + v9 + ", 10:" + v10 + ", 16:" + v16 + ", 19:" + v19);
+#print (" Sample = 1:" + w1 + ", 2:" + w2 + ", 3:" + w3 + ", 4:" + w4);
+## + ", 5:" + w5 + ", 8:" + w8 + ", 9:" + w9 + ", 10:" + w10 + ", 16:" + w16 + ", 19:" + w19);
+
+}
+
+print ("Average observed loss = " + (sum_of_observed_losses / num_of_observed_reports));
+print ("Writing out the results...");
+
+avg_reports_matrix = sum_of_observed_reports / num_of_observed_reports;
+avg_params = sum_of_observed_params / num_of_observed_reports;
+write (avg_reports_matrix, $11, format="text");
+write (avg_params, $12, format="text");
+
+print ("END ImputeGaussMCMC");
http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/05d2c0a8/src/test/scripts/applications/impute/old/imputeGaussMCMC.dml
----------------------------------------------------------------------
diff --git a/src/test/scripts/applications/impute/old/imputeGaussMCMC.dml b/src/test/scripts/applications/impute/old/imputeGaussMCMC.dml
index 3b50c7e..77bd21c 100644
--- a/src/test/scripts/applications/impute/old/imputeGaussMCMC.dml
+++ b/src/test/scripts/applications/impute/old/imputeGaussMCMC.dml
@@ -1,420 +1,420 @@
-#-------------------------------------------------------------
-#
-# 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.
-#
-#-------------------------------------------------------------
-
-# Implements the MCMC algorithm for imputation of missing data into a time-series of "reports".
-# Each report is a fixed-size vector of attribute values; reports come out each year/quarter/month ("term").
-# Hard linear equality constraints restrict values in/across the reports, e.g. total cost = sum of all costs.
-# Soft linear regression constraints define dependencies between values in/across the reports.
-# Linear regression parameters are unknown and sampled together with the missing values in the reports.
-#
-# INPUT 1: Initial reports matrix [1 : num_attrs, 1 : num_terms] with missing values usually set to zero,
-# but it MUST BE CONSISTENT with hard constraints! Set some missing values to nonzero if needed.
-# There are "num_terms" reports in the matrix, each having "num_attrs" attribute values.
-#
-# INPUT 2: Sparse matrix [1 : (num_terms * num_attrs), 1 : num_frees] that defines a linear map from
-# "free" variables to the reports matrix. A tensor of size (num_terms * num_attrs * num_frees)
-# where the reports matrix is stretched into a column-vector [1 : (num_terms * num_attrs)].
-# Term = t, attribute = i --> index = (t-1) * num_attrs + i
-#
-# INPUT 3: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : (num_terms * num_attrs)] that defines
-# a linear map from the stretched matrix of reports to the stretched matrix of regression factors.
-#
-# INPUT 4: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines "default regression factors"
-# (if nonzero) to be added to the regression factors before they are multiplied by parameters.
-#
-# INPUT 5: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : num_params] that defines a linear map
-# from the vector of parameters to the stretched matrix of regression factors.
-#
-# INPUT 6: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines default regression coefficients
-# (if nonzero) to be added to the parameters (if any) before being multiplied by regression factors.
-#
-# INPUT 7: A vector [1 : num_reg_eqs, 1] of scale multipliers, one per regression
-#
-# INPUT 8 : Number of factors in a regression equation, including the estimated value
-# INPUT 9 : Number of burn-in full iterations (that sample each variable and each parameter once)
-# INPUT 10: Number of observed full iterations (that sample each variable and each parameter once)
-#
-# INPUT 11: Output file name and path for the average MCMC reports table
-# INPUT 12: Output file for debugging (currently: the average parameters vector)
-#
-# Example:
-# hadoop jar SystemML.jar -f test/scripts/applications/impute/imputeGaussMCMC.dml -exec singlenode -args
-# test/scripts/applications/impute/initial_reports
-# test/scripts/applications/impute/CReps
-# test/scripts/applications/impute/RegresValueMap
-# test/scripts/applications/impute/RegresFactorDefault
-# test/scripts/applications/impute/RegresParamMap
-# test/scripts/applications/impute/RegresCoeffDefault
-# test/scripts/applications/impute/RegresScaleMult
-# 4 200 1000
-# test/scripts/applications/impute/output_reports
-# test/scripts/applications/impute/debug_info
-
-
-print ("START ImputeGaussMCMC");
-print ("Reading the input files...");
-
-initial_reports = read ($1);
-dReps = vectorize (initial_reports);
-CReps = read ($2);
-
-num_terms = ncol (initial_reports); # Number of periods for which reports are produced
-num_attrs = nrow (initial_reports); # Number of attribute values per each term report
-num_frees = ncol (CReps); # Number of free variables used to describe all consistent reports
-
-# We assume that all report-series consistent with hard constraints form an affine set:
-# reports = CReps %*% freeVars + dReps
-# Here "freeVars" is a vector of "free variables" (degrees of freedom), "CReps" is a linear mapping,
-# and "dReps" are the default values for the reports that correspond to the initial reports matrix.
-
-RegresValueMap = read ($3);
-RegresFactorDefault = read ($4);
-RegresParamMap = read ($5);
-RegresCoeffDefault = read ($6);
-RegresScaleMult = read ($7);
-tRegresScaleMult = t(RegresScaleMult);
-
-num_factors = $8; # Number of factors in each regression equation, including the estimated value
-num_reg_eqs = nrow (RegresParamMap) / num_factors; # Number of regression equations
-num_params = ncol (RegresParamMap); # Number of parameters used in all regressions
-num_burnin_iterations = $9;
-num_observed_iterations = $10;
-
-num_opt_iter = 20;
-print ("Performing initial optimization (" + num_opt_iter + " alternating CG steps)...");
-
-freeVars = matrix (0.0, rows = num_frees, cols = 1);
-params = matrix (1.0, rows = num_params, cols = 1);
-reports = CReps %*% freeVars + dReps;
-
-regresValues = RegresValueMap %*% reports + RegresFactorDefault;
-regresParams = RegresParamMap %*% params + RegresCoeffDefault;
-bilinear_vector = regresValues * regresParams;
-bilinear_form = matricize (bilinear_vector, num_factors);
-bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
-
-opt_iter = 1;
-is_step_params = 1;
-is_opt_converged = 0;
-
-print ("Before optimization: Initial bilinear form value = " + bilinear_form_value);
-
-
-while (is_opt_converged == 0)
-{
- deg = is_step_params * num_params + (1 - is_step_params) * num_frees;
-
- # Compute gradient
-
- gradient = matrix (0.0, rows = deg, cols = 1);
- for (i in 1:deg)
- {
- if (is_step_params == 1) {
- bilinear_vector = regresValues * RegresParamMap [, i];
- } else {
- bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
- }
- bilinear_updater = matricize (bilinear_vector, num_factors);
- q_minus_1 = sum (tRegresScaleMult * colSums (bilinear_form - bilinear_updater) * colSums (bilinear_form - bilinear_updater));
- q_plus_1 = sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater) * colSums (bilinear_form + bilinear_updater));
- gradient [i, 1] = 0.5 * (q_plus_1 - q_minus_1) + gradient [i, 1];
- }
-
- # Make a few conjugate gradient steps
-
- shift_vector = matrix (0.0, rows = deg, cols = 1);
- residual = gradient;
- p = - residual;
- norm_r2 = sum (residual * residual);
-
- for (j in 1:3)
- {
- q = matrix (0.0, rows = deg, cols = 1);
- for (i in 1:deg) # Want: q = A %*% p;
- {
- if (is_step_params == 1) {
- bilinear_vector = regresValues * RegresParamMap [, i];
- } else {
- bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
- }
- bilinear_updater_1 = matricize (bilinear_vector, num_factors);
-
- if (is_step_params == 1) {
- bilinear_vector = regresValues * (RegresParamMap %*% p);
- } else {
- bilinear_vector = (RegresValueMap %*% CReps %*% p) * regresParams;
- }
- bilinear_updater_p = matricize (bilinear_vector, num_factors);
-
- quadratic_plus_1 =
- sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_1) * colSums (bilinear_form + bilinear_updater_1));
- quadratic_plus_p =
- sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_p) * colSums (bilinear_form + bilinear_updater_p));
- quadratic_plus_both =
- sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p) * colSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p));
- q [i, 1] = (quadratic_plus_both - quadratic_plus_1 - quadratic_plus_p + bilinear_form_value) + q [i, 1];
- }
-
- alpha = norm_r2 / castAsScalar (t(p) %*% q);
- shift_vector = shift_vector + alpha * p;
- old_norm_r2 = norm_r2;
- residual = residual + alpha * q;
- norm_r2 = sum (residual * residual);
- p = - residual + (norm_r2 / old_norm_r2) * p;
- }
-
- if (is_step_params == 1) {
- params = params + shift_vector;
- regresParams = RegresParamMap %*% params + RegresCoeffDefault;
- } else {
- freeVars = freeVars + shift_vector;
- reports = CReps %*% freeVars + dReps;
- regresValues = RegresValueMap %*% reports + RegresFactorDefault;
- }
-
- # Update the bilinear form and check convergence
-
- if (is_step_params == 1) {
- old_bilinear_form_value = bilinear_form_value;
- }
- bilinear_vector = regresValues * regresParams;
- bilinear_form = matricize (bilinear_vector, num_factors);
- bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
-
- if (is_step_params == 1) {
- print ("Optimization step " + opt_iter + " (params) : bilinear form value = " + bilinear_form_value);
- } else {
- print ("Optimization step " + opt_iter + " (reports): bilinear form value = " + bilinear_form_value);
- }
-
- is_step_params = 1 - is_step_params;
- opt_iter = opt_iter + 1;
-
- if (is_step_params == 1 & opt_iter > num_opt_iter) {
- is_opt_converged = 1;
- }
-}
-
-print ("Performing MCMC initialization...");
-
-num_iter = num_burnin_iterations + num_observed_iterations;
-dim_sample = num_frees + num_params;
-sample_ones = matrix (1.0, rows = dim_sample, cols = 1);
-
-# Generate a random permutation matrix for the sampling order of freeVars and params
-
-SampleOrder = diag (sample_ones);
-num_swaps = 10 * dim_sample;
-rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
-left_swap = round (0.5 + dim_sample * rnd);
-rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
-right_swap = round (0.5 + dim_sample * rnd);
-for (swap_i in 1:num_swaps) {
- l = castAsScalar (left_swap [swap_i, 1]);
- r = castAsScalar (right_swap [swap_i, 1]);
- if (l != r) {
- tmp_row = SampleOrder [l, ];
- SampleOrder [l, ] = SampleOrder [r, ];
- SampleOrder [r, ] = tmp_row;
- }
-}
-
-isVar = colSums (SampleOrder [1 : num_frees, ]);
-sum_of_observed_reports = matrix (0.0, rows = num_attrs, cols = num_terms);
-sum_of_observed_params = matrix (0.0, rows = num_params, cols = 1);
-
-num_of_observed_reports = 0;
-pi = 3.1415926535897932384626433832795;
-zero = matrix (0.0, rows = 1, cols = 1);
-
-# Starting MCMC iterations
-
-for (iter in 1:num_iter)
-{
- # Initialize (bi-)linear forms
-
- regresValues = RegresValueMap %*% reports + RegresFactorDefault;
- regresParams = RegresParamMap %*% params + RegresCoeffDefault;
- bilinear_form_vector = regresValues * regresParams;
- bilinear_form = matricize (bilinear_form_vector, num_factors);
- bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
-
- if (iter <= num_burnin_iterations) {
- print ("MCMC iteration " + iter + " (burn-in) : bilinear form value = " + bilinear_form_value);
- } else {
- print ("MCMC iteration " + iter + " (observed): bilinear form value = " + bilinear_form_value);
- }
-
- # Create a normally distributed random sample
-
- dim_half_sample = castAsScalar (round (dim_sample / 2 + 0.1 + zero));
- rnd1 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
- rnd2 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
- rnd_normal_1 = sqrt (- 2.0 * log (rnd1)) * sin (2 * pi * rnd2);
- rnd_normal_2 = sqrt (- 2.0 * log (rnd1)) * cos (2 * pi * rnd2);
- rnd_normal = matrix (0.0, rows = dim_sample, cols = 1);
- rnd_normal [1 : dim_half_sample, ] = rnd_normal_1;
- rnd_normal [(dim_sample - dim_half_sample + 1) : dim_sample, ] = rnd_normal_2;
-
- # Initialize updaters
-
- freeVars_updater = freeVars * 0.0;
- params_updater = params * 0.0;
- regresValues_updater = regresValues * 0.0;
- regresParams_updater = regresParams * 0.0;
- bilinear_updater_vector = bilinear_form_vector * 0.0;
-
- # Perform the sampling
-
- for (idx in 1:dim_sample)
- {
- # Generate the sample unit-vector and updaters
-
- if (castAsScalar (isVar [1, idx]) > 0.5) {
- freeVars_updater = SampleOrder [1 : num_frees, idx];
- regresValues_updater = RegresValueMap %*% CReps %*% freeVars_updater;
- bilinear_updater_vector = regresValues_updater * regresParams;
- } else {
- params_updater = SampleOrder [(num_frees + 1) : dim_sample, idx];
- regresParams_updater = RegresParamMap %*% params_updater;
- bilinear_updater_vector = regresValues * regresParams_updater;
- }
- bilinear_updater = matricize (bilinear_updater_vector, num_factors);
-
- # Compute the quadratic by three shift-points: -1, 0, +1
-
- bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
- q_minus_1 = sum (tRegresScaleMult * colSums (bilinear_form - bilinear_updater) * colSums (bilinear_form - bilinear_updater));
- q_plus_1 = sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater) * colSums (bilinear_form + bilinear_updater));
- coeff_b = (q_plus_1 - q_minus_1) / 2.0;
- coeff_a = (q_plus_1 + q_minus_1) / 2.0 - bilinear_form_value;
-
-# BEGIN DEBUG INSERT
-# mmm = 1;
-# if (castAsScalar (isVar [1, idx]) > 0.5) {
-# for (iii in 2:num_frees) {if (castAsScalar (freeVars_updater [iii, 1] - freeVars_updater [mmm, 1]) > 0) {mmm = iii;}}
-# print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a);
-# } else {
-# for (iii in 2:num_params) {if (castAsScalar (params_updater [iii, 1] - params_updater [mmm, 1]) > 0) {mmm = iii;}}
-# print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a);
-# }
-# END DEBUG INSERT
-
- # Find the mean and the sigma for f(x) ~ exp (- (ax^2 + bx + c)),
- # then compute the shift to get the new sample
-
- mean_shift = - coeff_b / (2.0 * coeff_a);
- sigma_shift = 1.0 / sqrt (2.0 * coeff_a);
- shift = mean_shift + sigma_shift * castAsScalar (rnd_normal [idx, 1]);
-
- # Perform the updates
-
- bilinear_form = bilinear_form + shift * bilinear_updater;
- if (castAsScalar (isVar [1, idx]) > 0.5) {
- freeVars = freeVars + shift * freeVars_updater;
- regresValues = regresValues + shift * regresValues_updater;
- } else {
- params = params + shift * params_updater;
- regresParams = regresParams + shift * regresParams_updater;
- }
- }
-
- # Update / adjust the reports and the parameters
-
- reports = CReps %*% freeVars + dReps;
- reports_matrix = matricize (reports, num_attrs);
-
- # Make an observation of the reports and/or the parameters
-
- if (iter > num_burnin_iterations)
- {
- sum_of_observed_reports = sum_of_observed_reports + reports_matrix;
- num_of_observed_reports = num_of_observed_reports + 1;
-
- sum_of_observed_params = sum_of_observed_params + params;
- }
-
-
-v1 =castAsScalar(round(reports_matrix[ 1,num_terms]));
-v2 =castAsScalar(round(reports_matrix[ 2,num_terms]));
-v3 =castAsScalar(round(reports_matrix[ 3,num_terms]));
-v4 =castAsScalar(round(reports_matrix[ 4,num_terms]));
-v5 =castAsScalar(round(reports_matrix[ 5,num_terms]));
-v8 =castAsScalar(round(reports_matrix[ 8,num_terms]));
-v9 =castAsScalar(round(reports_matrix[ 9,num_terms]));
-v10=castAsScalar(round(reports_matrix[10,num_terms]));
-v16=castAsScalar(round(reports_matrix[16,num_terms]));
-v19=castAsScalar(round(reports_matrix[19,num_terms]));
-print (
-" Sample = 1:" + v1 + ", 2:" + v2 + ", 3:" + v3 + ", 4:" + v4 + ", 5:" + v5 +
-", 8:" + v8 + ", 9:" + v9 + ", 10:" + v10 + ", 16:" + v16 + ", 19:" + v19);
-
-}
-
-print ("Writing out the results...");
-
-avg_reports_matrix = sum_of_observed_reports / num_of_observed_reports;
-avg_params = sum_of_observed_params / num_of_observed_reports;
-write (avg_reports_matrix, $11, format="text");
-write (avg_params, $12, format="text");
-
-print ("END ImputeGaussMCMC");
-
-
-
-
-# Outputs a column vector that consists of the columns of the input matrix in sequential order
-# NEEDS TO BE PARALLELIZED
-vectorize = function (Matrix[double] M) return (Matrix[double] v)
-{
- n_rows = nrow (M);
- n_cols = ncol (M);
- n = n_rows * n_cols;
- v = matrix (0.0, rows = n, cols = 1);
- for (i in 1:n_cols) {
- left_row = n_rows * (i-1) + 1;
- right_row = n_rows * i;
- v [left_row : right_row, 1] = M [, i];
- }
-}
-
-# Takes a column vector, splits it into columns of "n_rows" in each, and combines into a matrix
-# NEEDS TO BE PARALLELIZED
-matricize = function (Matrix[double] v, int n_rows) return (Matrix[double] M)
-{
- zero = matrix (0.0, rows = 1, cols = 1);
- n = nrow (v);
- n_cols = castAsScalar (round (zero + (n / n_rows)));
- if (n_cols * n_rows < n) {
- n_cols = n_cols + 1;
- }
- M = matrix (0.0, rows = n_rows, cols = n_cols);
- for (i in 1:n_cols) {
- left_row = n_rows * (i-1) + 1;
- right_row = n_rows * i;
- if (right_row <= n) {
- M [, i] = v [left_row : right_row, 1];
- } else {
- new_right = n - left_row + 1;
- M [1 : new_right, i] = v [left_row : n, 1];
- }
- }
-}
+#-------------------------------------------------------------
+#
+# 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.
+#
+#-------------------------------------------------------------
+
+# Implements the MCMC algorithm for imputation of missing data into a time-series of "reports".
+# Each report is a fixed-size vector of attribute values; reports come out each year/quarter/month ("term").
+# Hard linear equality constraints restrict values in/across the reports, e.g. total cost = sum of all costs.
+# Soft linear regression constraints define dependencies between values in/across the reports.
+# Linear regression parameters are unknown and sampled together with the missing values in the reports.
+#
+# INPUT 1: Initial reports matrix [1 : num_attrs, 1 : num_terms] with missing values usually set to zero,
+# but it MUST BE CONSISTENT with hard constraints! Set some missing values to nonzero if needed.
+# There are "num_terms" reports in the matrix, each having "num_attrs" attribute values.
+#
+# INPUT 2: Sparse matrix [1 : (num_terms * num_attrs), 1 : num_frees] that defines a linear map from
+# "free" variables to the reports matrix. A tensor of size (num_terms * num_attrs * num_frees)
+# where the reports matrix is stretched into a column-vector [1 : (num_terms * num_attrs)].
+# Term = t, attribute = i --> index = (t-1) * num_attrs + i
+#
+# INPUT 3: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : (num_terms * num_attrs)] that defines
+# a linear map from the stretched matrix of reports to the stretched matrix of regression factors.
+#
+# INPUT 4: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines "default regression factors"
+# (if nonzero) to be added to the regression factors before they are multiplied by parameters.
+#
+# INPUT 5: Sparse matrix [1 : (num_reg_eqs * num_factors), 1 : num_params] that defines a linear map
+# from the vector of parameters to the stretched matrix of regression factors.
+#
+# INPUT 6: Sparse vector [1 : (num_reg_eqs * num_factors), 1] that defines default regression coefficients
+# (if nonzero) to be added to the parameters (if any) before being multiplied by regression factors.
+#
+# INPUT 7: A vector [1 : num_reg_eqs, 1] of scale multipliers, one per regression
+#
+# INPUT 8 : Number of factors in a regression equation, including the estimated value
+# INPUT 9 : Number of burn-in full iterations (that sample each variable and each parameter once)
+# INPUT 10: Number of observed full iterations (that sample each variable and each parameter once)
+#
+# INPUT 11: Output file name and path for the average MCMC reports table
+# INPUT 12: Output file for debugging (currently: the average parameters vector)
+#
+# Example:
+# hadoop jar SystemML.jar -f test/scripts/applications/impute/imputeGaussMCMC.dml -exec singlenode -args
+# test/scripts/applications/impute/initial_reports
+# test/scripts/applications/impute/CReps
+# test/scripts/applications/impute/RegresValueMap
+# test/scripts/applications/impute/RegresFactorDefault
+# test/scripts/applications/impute/RegresParamMap
+# test/scripts/applications/impute/RegresCoeffDefault
+# test/scripts/applications/impute/RegresScaleMult
+# 4 200 1000
+# test/scripts/applications/impute/output_reports
+# test/scripts/applications/impute/debug_info
+
+
+print ("START ImputeGaussMCMC");
+print ("Reading the input files...");
+
+initial_reports = read ($1);
+dReps = vectorize (initial_reports);
+CReps = read ($2);
+
+num_terms = ncol (initial_reports); # Number of periods for which reports are produced
+num_attrs = nrow (initial_reports); # Number of attribute values per each term report
+num_frees = ncol (CReps); # Number of free variables used to describe all consistent reports
+
+# We assume that all report-series consistent with hard constraints form an affine set:
+# reports = CReps %*% freeVars + dReps
+# Here "freeVars" is a vector of "free variables" (degrees of freedom), "CReps" is a linear mapping,
+# and "dReps" are the default values for the reports that correspond to the initial reports matrix.
+
+RegresValueMap = read ($3);
+RegresFactorDefault = read ($4);
+RegresParamMap = read ($5);
+RegresCoeffDefault = read ($6);
+RegresScaleMult = read ($7);
+tRegresScaleMult = t(RegresScaleMult);
+
+num_factors = $8; # Number of factors in each regression equation, including the estimated value
+num_reg_eqs = nrow (RegresParamMap) / num_factors; # Number of regression equations
+num_params = ncol (RegresParamMap); # Number of parameters used in all regressions
+num_burnin_iterations = $9;
+num_observed_iterations = $10;
+
+num_opt_iter = 20;
+print ("Performing initial optimization (" + num_opt_iter + " alternating CG steps)...");
+
+freeVars = matrix (0.0, rows = num_frees, cols = 1);
+params = matrix (1.0, rows = num_params, cols = 1);
+reports = CReps %*% freeVars + dReps;
+
+regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+bilinear_vector = regresValues * regresParams;
+bilinear_form = matricize (bilinear_vector, num_factors);
+bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
+
+opt_iter = 1;
+is_step_params = 1;
+is_opt_converged = 0;
+
+print ("Before optimization: Initial bilinear form value = " + bilinear_form_value);
+
+
+while (is_opt_converged == 0)
+{
+ deg = is_step_params * num_params + (1 - is_step_params) * num_frees;
+
+ # Compute gradient
+
+ gradient = matrix (0.0, rows = deg, cols = 1);
+ for (i in 1:deg)
+ {
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * RegresParamMap [, i];
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
+ }
+ bilinear_updater = matricize (bilinear_vector, num_factors);
+ q_minus_1 = sum (tRegresScaleMult * colSums (bilinear_form - bilinear_updater) * colSums (bilinear_form - bilinear_updater));
+ q_plus_1 = sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater) * colSums (bilinear_form + bilinear_updater));
+ gradient [i, 1] = 0.5 * (q_plus_1 - q_minus_1) + gradient [i, 1];
+ }
+
+ # Make a few conjugate gradient steps
+
+ shift_vector = matrix (0.0, rows = deg, cols = 1);
+ residual = gradient;
+ p = - residual;
+ norm_r2 = sum (residual * residual);
+
+ for (j in 1:3)
+ {
+ q = matrix (0.0, rows = deg, cols = 1);
+ for (i in 1:deg) # Want: q = A %*% p;
+ {
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * RegresParamMap [, i];
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps [, i]) * regresParams;
+ }
+ bilinear_updater_1 = matricize (bilinear_vector, num_factors);
+
+ if (is_step_params == 1) {
+ bilinear_vector = regresValues * (RegresParamMap %*% p);
+ } else {
+ bilinear_vector = (RegresValueMap %*% CReps %*% p) * regresParams;
+ }
+ bilinear_updater_p = matricize (bilinear_vector, num_factors);
+
+ quadratic_plus_1 =
+ sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_1) * colSums (bilinear_form + bilinear_updater_1));
+ quadratic_plus_p =
+ sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_p) * colSums (bilinear_form + bilinear_updater_p));
+ quadratic_plus_both =
+ sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p) * colSums (bilinear_form + bilinear_updater_1 + bilinear_updater_p));
+ q [i, 1] = (quadratic_plus_both - quadratic_plus_1 - quadratic_plus_p + bilinear_form_value) + q [i, 1];
+ }
+
+ alpha = norm_r2 / castAsScalar (t(p) %*% q);
+ shift_vector = shift_vector + alpha * p;
+ old_norm_r2 = norm_r2;
+ residual = residual + alpha * q;
+ norm_r2 = sum (residual * residual);
+ p = - residual + (norm_r2 / old_norm_r2) * p;
+ }
+
+ if (is_step_params == 1) {
+ params = params + shift_vector;
+ regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+ } else {
+ freeVars = freeVars + shift_vector;
+ reports = CReps %*% freeVars + dReps;
+ regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+ }
+
+ # Update the bilinear form and check convergence
+
+ if (is_step_params == 1) {
+ old_bilinear_form_value = bilinear_form_value;
+ }
+ bilinear_vector = regresValues * regresParams;
+ bilinear_form = matricize (bilinear_vector, num_factors);
+ bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
+
+ if (is_step_params == 1) {
+ print ("Optimization step " + opt_iter + " (params) : bilinear form value = " + bilinear_form_value);
+ } else {
+ print ("Optimization step " + opt_iter + " (reports): bilinear form value = " + bilinear_form_value);
+ }
+
+ is_step_params = 1 - is_step_params;
+ opt_iter = opt_iter + 1;
+
+ if (is_step_params == 1 & opt_iter > num_opt_iter) {
+ is_opt_converged = 1;
+ }
+}
+
+print ("Performing MCMC initialization...");
+
+num_iter = num_burnin_iterations + num_observed_iterations;
+dim_sample = num_frees + num_params;
+sample_ones = matrix (1.0, rows = dim_sample, cols = 1);
+
+# Generate a random permutation matrix for the sampling order of freeVars and params
+
+SampleOrder = diag (sample_ones);
+num_swaps = 10 * dim_sample;
+rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
+left_swap = round (0.5 + dim_sample * rnd);
+rnd = Rand (rows = num_swaps, cols = 1, min = 0.0, max = 1.0);
+right_swap = round (0.5 + dim_sample * rnd);
+for (swap_i in 1:num_swaps) {
+ l = castAsScalar (left_swap [swap_i, 1]);
+ r = castAsScalar (right_swap [swap_i, 1]);
+ if (l != r) {
+ tmp_row = SampleOrder [l, ];
+ SampleOrder [l, ] = SampleOrder [r, ];
+ SampleOrder [r, ] = tmp_row;
+ }
+}
+
+isVar = colSums (SampleOrder [1 : num_frees, ]);
+sum_of_observed_reports = matrix (0.0, rows = num_attrs, cols = num_terms);
+sum_of_observed_params = matrix (0.0, rows = num_params, cols = 1);
+
+num_of_observed_reports = 0;
+pi = 3.1415926535897932384626433832795;
+zero = matrix (0.0, rows = 1, cols = 1);
+
+# Starting MCMC iterations
+
+for (iter in 1:num_iter)
+{
+ # Initialize (bi-)linear forms
+
+ regresValues = RegresValueMap %*% reports + RegresFactorDefault;
+ regresParams = RegresParamMap %*% params + RegresCoeffDefault;
+ bilinear_form_vector = regresValues * regresParams;
+ bilinear_form = matricize (bilinear_form_vector, num_factors);
+ bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
+
+ if (iter <= num_burnin_iterations) {
+ print ("MCMC iteration " + iter + " (burn-in) : bilinear form value = " + bilinear_form_value);
+ } else {
+ print ("MCMC iteration " + iter + " (observed): bilinear form value = " + bilinear_form_value);
+ }
+
+ # Create a normally distributed random sample
+
+ dim_half_sample = castAsScalar (round (dim_sample / 2 + 0.1 + zero));
+ rnd1 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
+ rnd2 = Rand (rows = dim_half_sample, cols = 1, min = 0.0, max = 1.0);
+ rnd_normal_1 = sqrt (- 2.0 * log (rnd1)) * sin (2 * pi * rnd2);
+ rnd_normal_2 = sqrt (- 2.0 * log (rnd1)) * cos (2 * pi * rnd2);
+ rnd_normal = matrix (0.0, rows = dim_sample, cols = 1);
+ rnd_normal [1 : dim_half_sample, ] = rnd_normal_1;
+ rnd_normal [(dim_sample - dim_half_sample + 1) : dim_sample, ] = rnd_normal_2;
+
+ # Initialize updaters
+
+ freeVars_updater = freeVars * 0.0;
+ params_updater = params * 0.0;
+ regresValues_updater = regresValues * 0.0;
+ regresParams_updater = regresParams * 0.0;
+ bilinear_updater_vector = bilinear_form_vector * 0.0;
+
+ # Perform the sampling
+
+ for (idx in 1:dim_sample)
+ {
+ # Generate the sample unit-vector and updaters
+
+ if (castAsScalar (isVar [1, idx]) > 0.5) {
+ freeVars_updater = SampleOrder [1 : num_frees, idx];
+ regresValues_updater = RegresValueMap %*% CReps %*% freeVars_updater;
+ bilinear_updater_vector = regresValues_updater * regresParams;
+ } else {
+ params_updater = SampleOrder [(num_frees + 1) : dim_sample, idx];
+ regresParams_updater = RegresParamMap %*% params_updater;
+ bilinear_updater_vector = regresValues * regresParams_updater;
+ }
+ bilinear_updater = matricize (bilinear_updater_vector, num_factors);
+
+ # Compute the quadratic by three shift-points: -1, 0, +1
+
+ bilinear_form_value = sum (tRegresScaleMult * colSums (bilinear_form) * colSums (bilinear_form));
+ q_minus_1 = sum (tRegresScaleMult * colSums (bilinear_form - bilinear_updater) * colSums (bilinear_form - bilinear_updater));
+ q_plus_1 = sum (tRegresScaleMult * colSums (bilinear_form + bilinear_updater) * colSums (bilinear_form + bilinear_updater));
+ coeff_b = (q_plus_1 - q_minus_1) / 2.0;
+ coeff_a = (q_plus_1 + q_minus_1) / 2.0 - bilinear_form_value;
+
+# BEGIN DEBUG INSERT
+# mmm = 1;
+# if (castAsScalar (isVar [1, idx]) > 0.5) {
+# for (iii in 2:num_frees) {if (castAsScalar (freeVars_updater [iii, 1] - freeVars_updater [mmm, 1]) > 0) {mmm = iii;}}
+# print ("freeVars[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a);
+# } else {
+# for (iii in 2:num_params) {if (castAsScalar (params_updater [iii, 1] - params_updater [mmm, 1]) > 0) {mmm = iii;}}
+# print (" params[" + mmm + "]: q_minus_1 = " + q_minus_1 + ", q_plus_1 = " + q_plus_1 + ", coeff_a = " + coeff_a);
+# }
+# END DEBUG INSERT
+
+ # Find the mean and the sigma for f(x) ~ exp (- (ax^2 + bx + c)),
+ # then compute the shift to get the new sample
+
+ mean_shift = - coeff_b / (2.0 * coeff_a);
+ sigma_shift = 1.0 / sqrt (2.0 * coeff_a);
+ shift = mean_shift + sigma_shift * castAsScalar (rnd_normal [idx, 1]);
+
+ # Perform the updates
+
+ bilinear_form = bilinear_form + shift * bilinear_updater;
+ if (castAsScalar (isVar [1, idx]) > 0.5) {
+ freeVars = freeVars + shift * freeVars_updater;
+ regresValues = regresValues + shift * regresValues_updater;
+ } else {
+ params = params + shift * params_updater;
+ regresParams = regresParams + shift * regresParams_updater;
+ }
+ }
+
+ # Update / adjust the reports and the parameters
+
+ reports = CReps %*% freeVars + dReps;
+ reports_matrix = matricize (reports, num_attrs);
+
+ # Make an observation of the reports and/or the parameters
+
+ if (iter > num_burnin_iterations)
+ {
+ sum_of_observed_reports = sum_of_observed_reports + reports_matrix;
+ num_of_observed_reports = num_of_observed_reports + 1;
+
+ sum_of_observed_params = sum_of_observed_params + params;
+ }
+
+
+v1 =castAsScalar(round(reports_matrix[ 1,num_terms]));
+v2 =castAsScalar(round(reports_matrix[ 2,num_terms]));
+v3 =castAsScalar(round(reports_matrix[ 3,num_terms]));
+v4 =castAsScalar(round(reports_matrix[ 4,num_terms]));
+v5 =castAsScalar(round(reports_matrix[ 5,num_terms]));
+v8 =castAsScalar(round(reports_matrix[ 8,num_terms]));
+v9 =castAsScalar(round(reports_matrix[ 9,num_terms]));
+v10=castAsScalar(round(reports_matrix[10,num_terms]));
+v16=castAsScalar(round(reports_matrix[16,num_terms]));
+v19=castAsScalar(round(reports_matrix[19,num_terms]));
+print (
+" Sample = 1:" + v1 + ", 2:" + v2 + ", 3:" + v3 + ", 4:" + v4 + ", 5:" + v5 +
+", 8:" + v8 + ", 9:" + v9 + ", 10:" + v10 + ", 16:" + v16 + ", 19:" + v19);
+
+}
+
+print ("Writing out the results...");
+
+avg_reports_matrix = sum_of_observed_reports / num_of_observed_reports;
+avg_params = sum_of_observed_params / num_of_observed_reports;
+write (avg_reports_matrix, $11, format="text");
+write (avg_params, $12, format="text");
+
+print ("END ImputeGaussMCMC");
+
+
+
+
+# Outputs a column vector that consists of the columns of the input matrix in sequential order
+# NEEDS TO BE PARALLELIZED
+vectorize = function (Matrix[double] M) return (Matrix[double] v)
+{
+ n_rows = nrow (M);
+ n_cols = ncol (M);
+ n = n_rows * n_cols;
+ v = matrix (0.0, rows = n, cols = 1);
+ for (i in 1:n_cols) {
+ left_row = n_rows * (i-1) + 1;
+ right_row = n_rows * i;
+ v [left_row : right_row, 1] = M [, i];
+ }
+}
+
+# Takes a column vector, splits it into columns of "n_rows" in each, and combines into a matrix
+# NEEDS TO BE PARALLELIZED
+matricize = function (Matrix[double] v, int n_rows) return (Matrix[double] M)
+{
+ zero = matrix (0.0, rows = 1, cols = 1);
+ n = nrow (v);
+ n_cols = castAsScalar (round (zero + (n / n_rows)));
+ if (n_cols * n_rows < n) {
+ n_cols = n_cols + 1;
+ }
+ M = matrix (0.0, rows = n_rows, cols = n_cols);
+ for (i in 1:n_cols) {
+ left_row = n_rows * (i-1) + 1;
+ right_row = n_rows * i;
+ if (right_row <= n) {
+ M [, i] = v [left_row : right_row, 1];
+ } else {
+ new_right = n - left_row + 1;
+ M [1 : new_right, i] = v [left_row : n, 1];
+ }
+ }
+}