Actual source code: ex3.c
1: #include <petscregressor.h>
3: static char help[] = "Tests some linear PetscRegressor types with different regularizers.\n\n";
5: typedef struct _AppCtx {
6: Mat X; /* Training data */
7: Vec y; /* Target data */
8: Vec y_predicted; /* Target data */
9: Vec coefficients;
10: PetscInt N; /* Data size */
11: PetscBool flg_string;
12: PetscBool flg_ascii;
13: PetscBool flg_view_sol;
14: PetscBool test_prefix;
15: } *AppCtx;
17: static PetscErrorCode DestroyCtx(AppCtx *ctx)
18: {
19: PetscFunctionBegin;
20: PetscCall(MatDestroy(&(*ctx)->X));
21: PetscCall(VecDestroy(&(*ctx)->y));
22: PetscCall(VecDestroy(&(*ctx)->y_predicted));
23: PetscCall(PetscFree(*ctx));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode TestRegressorViews(PetscRegressor regressor, AppCtx ctx)
28: {
29: PetscRegressorType check_type;
30: PetscBool match;
32: PetscFunctionBegin;
33: if (ctx->flg_view_sol) {
34: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Training target vector is\n"));
35: PetscCall(VecView(ctx->y, PETSC_VIEWER_STDOUT_WORLD));
36: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Predicted values are\n"));
37: PetscCall(VecView(ctx->y_predicted, PETSC_VIEWER_STDOUT_WORLD));
38: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coefficients are\n"));
39: PetscCall(VecView(ctx->coefficients, PETSC_VIEWER_STDOUT_WORLD));
40: }
42: if (ctx->flg_string) {
43: PetscViewer stringviewer;
44: char string[512];
45: const char *outstring;
47: PetscCall(PetscViewerStringOpen(PETSC_COMM_WORLD, string, sizeof(string), &stringviewer));
48: PetscCall(PetscRegressorView(regressor, stringviewer));
49: PetscCall(PetscViewerStringGetStringRead(stringviewer, &outstring, NULL));
50: PetscCheck((char *)outstring == (char *)string, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "String returned from viewer does not equal original string");
51: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output from string viewer:%s\n", outstring));
52: PetscCall(PetscViewerDestroy(&stringviewer));
53: } else if (ctx->flg_ascii) PetscCall(PetscRegressorView(regressor, PETSC_VIEWER_STDOUT_WORLD));
55: PetscCall(PetscRegressorGetType(regressor, &check_type));
56: PetscCall(PetscStrcmp(check_type, PETSCREGRESSORLINEAR, &match));
57: PetscCheck(match, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "Regressor type is not Linear");
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode TestPrefixRegressor(PetscRegressor regressor, AppCtx ctx)
62: {
63: PetscFunctionBegin;
64: if (ctx->test_prefix) {
65: PetscCall(PetscRegressorSetOptionsPrefix(regressor, "sys1_"));
66: PetscCall(PetscRegressorAppendOptionsPrefix(regressor, "sys2_"));
67: }
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: static PetscErrorCode CreateData(AppCtx ctx)
72: {
73: PetscMPIInt rank;
74: PetscScalar mean;
76: PetscFunctionBegin;
77: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
78: PetscCall(VecCreate(PETSC_COMM_WORLD, &ctx->y));
79: PetscCall(VecSetSizes(ctx->y, PETSC_DECIDE, ctx->N));
80: PetscCall(VecSetFromOptions(ctx->y));
81: PetscCall(VecDuplicate(ctx->y, &ctx->y_predicted));
82: PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx->X));
83: PetscCall(MatSetSizes(ctx->X, PETSC_DECIDE, PETSC_DECIDE, ctx->N, ctx->N));
84: PetscCall(MatSetFromOptions(ctx->X));
85: PetscCall(MatSetUp(ctx->X));
87: if (!rank) {
88: for (PetscInt i = 0; i < ctx->N; i++) {
89: PetscCall(VecSetValue(ctx->y, i, (PetscScalar)i, INSERT_VALUES));
90: PetscCall(MatSetValue(ctx->X, i, i, 1.0, INSERT_VALUES));
91: }
92: }
93: /* Set up a training data matrix that is the identity.
94: * We do this because this gives us a special case in which we can analytically determine what the regression
95: * coefficients should be for ordinary least squares, LASSO (L1 regularized), and ridge (L2 regularized) regression.
96: * See details in section 6.2 of James et al.'s An Introduction to Statistical Learning (ISLR), in the subsection
97: * titled "A Simple Special Case for Ridge Regression and the Lasso".
98: * Note that the coefficients we generate with ridge regression (-regressor_linear_type ridge -regressor_regularizer_weight <lambda>, or, equivalently,
99: * -tao_brgn_regularization_type l2pure -tao_brgn_regularizer_weight <lambda>) match those of the ISLR formula exactly.
100: * For LASSO it does not match the ISLR formula: where they use lambda/2, we need to use lambda.
101: * It also doesn't match what Scikit-learn does; in that case their lambda is 1/n_samples of our lambda. Apparently everyone is scaling
102: * their loss function by a different value, hence the need to change what "lambda" is. But it's clear that ISLR, Scikit-learn, and we
103: * are basically doing the same thing otherwise. */
104: PetscCall(VecAssemblyBegin(ctx->y));
105: PetscCall(VecAssemblyEnd(ctx->y));
106: PetscCall(MatAssemblyBegin(ctx->X, MAT_FINAL_ASSEMBLY));
107: PetscCall(MatAssemblyEnd(ctx->X, MAT_FINAL_ASSEMBLY));
108: /* Center the target vector we will train with. */
109: PetscCall(VecMean(ctx->y, &mean));
110: PetscCall(VecShift(ctx->y, -1.0 * mean));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode ConfigureContext(AppCtx ctx)
115: {
116: PetscFunctionBegin;
117: ctx->flg_string = PETSC_FALSE;
118: ctx->flg_ascii = PETSC_FALSE;
119: ctx->flg_view_sol = PETSC_FALSE;
120: ctx->test_prefix = PETSC_FALSE;
121: ctx->N = 10;
122: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Options for PetscRegressor ex3:", "");
123: PetscCall(PetscOptionsInt("-N", "Dimension of the N x N data matrix", "ex3.c", ctx->N, &ctx->N, NULL));
124: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_string_viewer", &ctx->flg_string, NULL));
125: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_ascii_viewer", &ctx->flg_ascii, NULL));
126: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_sols", &ctx->flg_view_sol, NULL));
127: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_prefix", &ctx->test_prefix, NULL));
128: PetscOptionsEnd();
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: int main(int argc, char **args)
133: {
134: AppCtx ctx;
135: PetscRegressor regressor;
136: PetscScalar intercept;
138: /* Initialize PETSc */
139: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
141: /* Initialize problem parameters and data */
142: PetscCall(PetscNew(&ctx));
143: PetscCall(ConfigureContext(ctx));
144: PetscCall(CreateData(ctx));
146: /* Create Regressor solver with desired type and options */
147: PetscCall(PetscRegressorCreate(PETSC_COMM_WORLD, ®ressor));
148: PetscCall(PetscRegressorSetType(regressor, PETSCREGRESSORLINEAR));
149: PetscCall(PetscRegressorLinearSetType(regressor, REGRESSOR_LINEAR_OLS));
150: PetscCall(PetscRegressorLinearSetFitIntercept(regressor, PETSC_FALSE));
151: /* Testing prefix functions for Regressor */
152: PetscCall(TestPrefixRegressor(regressor, ctx));
153: /* Check for command line options */
154: PetscCall(PetscRegressorSetFromOptions(regressor));
155: /* Fit the regressor */
156: PetscCall(PetscRegressorFit(regressor, ctx->X, ctx->y));
157: /* Predict data with fitted regressor */
158: PetscCall(PetscRegressorPredict(regressor, ctx->X, ctx->y_predicted));
159: /* Get other desired output data */
160: PetscCall(PetscRegressorLinearGetIntercept(regressor, &intercept));
161: PetscCall(PetscRegressorLinearGetCoefficients(regressor, &ctx->coefficients));
163: /* Testing Views, and GetTypes */
164: PetscCall(TestRegressorViews(regressor, ctx));
165: PetscCall(PetscRegressorDestroy(®ressor));
166: PetscCall(DestroyCtx(&ctx));
167: PetscCall(PetscFinalize());
168: return 0;
169: }
171: /*TEST
173: build:
174: requires: !complex !single !__float128 !defined(PETSC_USE_64BIT_INDICES)
176: test:
177: suffix: prefix_tao
178: args: -sys1_sys2_regressor_view ::ascii_info_detail -test_prefix
180: test:
181: suffix: prefix_ksp
182: args: -sys1_sys2_regressor_view -test_prefix -sys1_sys2_regressor_linear_use_ksp -sys1_sys2_regressor_linear_ksp_monitor
184: test:
185: suffix: prefix_ksp_cholesky
186: args: -sys1_sys2_regressor_view -test_prefix -sys1_sys2_regressor_linear_use_ksp -sys1_sys2_regressor_linear_pc_type cholesky
187: TODO: Could not locate a solver type for factorization type CHOLESKY and matrix type normal
189: test:
190: suffix: prefix_ksp_suitesparse
191: requires: suitesparse
192: args: -sys1_sys2_regressor_view -test_prefix -sys1_sys2_regressor_linear_use_ksp -sys1_sys2_regressor_linear_pc_type qr -sys1_sys2_regressor_linear_pc_factor_mat_solver_type spqr -sys1_sys2_regressor_linear_ksp_monitor
194: test:
195: suffix: asciiview
196: args: -test_ascii_viewer
198: test:
199: suffix: stringview
200: args: -test_string_viewer
202: test:
203: suffix: ksp_intercept
204: args: -regressor_linear_use_ksp -regressor_linear_fit_intercept -regressor_view
206: test:
207: suffix: ksp_no_intercept
208: args: -regressor_linear_use_ksp -regressor_view
210: test:
211: suffix: lasso_1
212: nsize: 1
213: args: -regressor_type linear -regressor_linear_type lasso -regressor_regularizer_weight 2 -regressor_linear_fit_intercept -view_sols
215: test:
216: suffix: lasso_2
217: nsize: 2
218: args: -regressor_type linear -regressor_linear_type lasso -regressor_regularizer_weight 2 -regressor_linear_fit_intercept -view_sols
220: test:
221: suffix: ridge_1
222: nsize: 1
223: args: -regressor_type linear -regressor_linear_type ridge -regressor_regularizer_weight 2 -regressor_linear_fit_intercept -view_sols
225: test:
226: suffix: ridge_2
227: nsize: 2
228: args: -regressor_type linear -regressor_linear_type ridge -regressor_regularizer_weight 2 -regressor_linear_fit_intercept -view_sols
230: test:
231: suffix: ols_1
232: nsize: 1
233: args: -regressor_type linear -regressor_linear_type ols -regressor_linear_fit_intercept -view_sols
235: test:
236: suffix: ols_2
237: nsize: 2
238: args: -regressor_type linear -regressor_linear_type ols -regressor_linear_fit_intercept -view_sols
240: TEST*/