From c6688ca94f700c75762559a698c7d74ee5a41b78 Mon Sep 17 00:00:00 2001 From: Jacob Smith <jake@laptop.(none)> Date: Fri, 23 Apr 2010 09:57:40 -0400 Subject: [PATCH] Added NPSOL --- CFiles/optpackage.c | 383 ++++++++++++++++++++++++++++++++++++++++++++ CFiles/optpackage.h | 79 +++++++++ 2 files changed, 462 insertions(+) diff --git a/CFiles/optpackage.c b/CFiles/optpackage.c index efd8521..1da1647 100644 --- a/CFiles/optpackage.c +++ b/CFiles/optpackage.c @@ -25,6 +25,17 @@ static void gradcd_imslconlin(int n, double *x, double *g); static double ObjFuncForModel_congrad(double *x0_p, int d_x0); +// NPSOL related +static TSdvector *XNPSOL_DV = NULL; //To save the minimized value in case the IMSL quits with a higher value. +void confun(int *mode, int *ncnln, int *n, int *ldJ, int *needc, double *x, double *c, double *cJac, int *nstate); +void npsol(int *n, int *nclin, int *ncnln, int *ldA, int *ldJ, int *ldR, double *A, double *bl, double *bu, + void (*funcon)(int *mode, int *ncnln, int *n, int *ldJ, int *needc, double *x, double *c, double *cJac, int *nstate), + void (*funobj)(int *mode, int *n, double *x, double *f, double *g, int *nstate), + int *inform, int *iter, int *istate, double *c, double *cJac, double *clamda, double *f, double *g, double *R, double *x, + int *iw, int *leniw, double *w, int *lenw); +static void ObjFuncForModel_npsolconlin(int *mode, int *N, double *x, double *f, double *g, int *nstate); +static struct TStateModel_tag *SetModelGlobalForNPSOLconlin(struct TStateModel_tag *smodel_ps); +static void npsolconlin_SetPrintFile(char *filename); ////TSminpack *CreateTSminpack(TFminpackage *minfinder_func, TFminobj *minobj_func, TFmingrad *mingrad_func, TFSetPrintFile *printinterresults_func, const int n, const int package) //, const int indxAnag) ////TSminpack *CreateTSminpack(TFminfinder *minfinder_func, TFminobj *minobj_func, TFmingrad *mingrad_func, char *filename_printout, const int n, const int package) //, const int indxAnag) @@ -1111,3 +1122,375 @@ static struct TSetc_imslconlin_tag *CreateTSetc_imslconlin(struct TSminpack_tag } /**/ + + +//----------------------------------------------------------------------- +// Using Linearly-constrained NPSOL minimization package. +// Added: 4/21/10 - Jake Smith +//----------------------------------------------------------------------- +void minfinder_noblocknpsolconlin(struct TSpackage_npsolconlin_tag *package_npsolconlin_ps, struct TSminpack_tag *minpack_ps, char *filename_printout, int ntheta) +{ +#ifdef _NPSOL_ + + + //ntheta: number of free model parameters (NOT including free transition matrix Q parameters). + //filename_printout: the file that stores the intermediate results. + + int i, j; + + //--- Model or project specific structure. + struct TStateModel_tag *smodel_ps = ((struct TSetc_minproj_tag *)minpack_ps->etc_project_ps)->smodel_ps; + //--- + + + TSdvector *x_dv = minpack_ps->x_dv; + TSdvector *g_dv = minpack_ps->g_dv; + double *x1_pd=NULL, *x2_pd=NULL; + //=== + TSdvector *xguess_dv = CreateVector_lf(x_dv->n); + + x1_pd = x_dv->v; //In the constant parameter model, this will point to invalid, + x2_pd = x_dv->v + ntheta; // but will be taken care of automatically by DW's function ConvertFreeParametersToQ(). + + CopyVector0(xguess_dv, x_dv); + + //======= NPSOL linearly-constrained optimization, which makes sure that the boundary condition is met. + npsolconlin_SetPrintFile(filename_printout); //Set the print-out file outputsp_min_tag.prn. + printf("\n\n======= Starting the NPSOL constrained optimization======= \n\n"); + fflush(stdout); + //====== The following linearly-constrained minimization works well for this kind of model but has a bugger of returning a higher value of the objective function. + CopyVector0(XNPSOL_DV, x_dv); //This is absolutely necessary because once imsl_d_min_con_gen_lin() is called, x_dv will be + // changed before ObjFuncForModel_imslconlin() is evaluated. It is possible that x_dv is changed + // so much that bad objective is returned and thus XIMSL_DV would be bad from the start, thus + // giving 1.eE+3000 from beginning to end. + GLB_FVALMIN = -LogPosterior_StatesIntegratedOut(smodel_ps); + CopyVector0(package_npsolconlin_ps->xsaved_dv, XNPSOL_DV); + //+ + SetModelGlobalForNPSOLconlin(smodel_ps); + + //////////////////////////////////////////////////////////////// + // Set up NPSOL variables + //////////////////////////////////////////////////////////////// + int N = package_npsolconlin_ps->n; + int nlin = package_npsolconlin_ps->num_lin; + int nclin = package_npsolconlin_ps->num_nlin; + int nbnd = N + nlin + nclin; + + double x[N]; + double g[N]; + //double *x = malloc(sizeof(double)*N); + //double *g = malloc(sizeof(double)*N); + + int ldJ = nclin; + if (nclin <= 0) + ldJ = 1; + + int ldR = N; + int ldA = nlin; + + double A[50][50]; + + // convert row major definition to column major + if (nlin > 0) { + TSdmatrix *m = package_npsolconlin_ps->A; + for (i=0; i<N; i++) { + for (j=0; j<ldA; j++) { + //TODO: Talk to Tao about how he wants this done + //A[i][j] = m->M[j][i]; + } + } + } + else + { + ldA = 1; + for (i=0; i<N; i++) + A[i][0] = 0; + } + + double lb[150]; + double ub[150]; + + + // the first part of ub/lb are the general bound constraints + for (i=0; i<N;i++) { + lb[i] = package_npsolconlin_ps->lb->v[i]; + ub[i] = package_npsolconlin_ps->ub->v[i]; + } + // the second part of the ub/lb are the linear constraints; + for (i=N; i<N+nlin; i++) { + lb[i] = package_npsolconlin_ps->lin_lb->v[i-N]; + ub[i] = package_npsolconlin_ps->lin_ub->v[i-N]; + } + // the third part of the ub/lb are the nonlinear constraints; + for (i=N+nlin; i<nbnd; i++) { + lb[i] = package_npsolconlin_ps->nlin_lb->v[i-N-nlin]; + ub[i] = package_npsolconlin_ps->nlin_ub->v[i-N-nlin]; + } + + double result; + double *c = NULL; + c = malloc(sizeof(double)*nclin); + + int inform; + + // workspace + double *R = malloc(ldR*N*sizeof(double)); + int itera; + int istate[nbnd]; + int leniw = 3*N + nlin + 2*nclin; + int lenw = 2*pow(N,2) + N*nlin + 2*N*nclin + 20*N + 11*nlin + 21*nclin; + int *iw = malloc(sizeof(int)*leniw); + double *w = malloc(sizeof(double)*lenw); + double *cJac = malloc(sizeof(double)*1); + double *clambda = malloc(sizeof(double)*nbnd); + + //////////////////////////////////////////////////////////////// + // CALL NPSOL + //////////////////////////////////////////////////////////////// + + npsol((int*) &N, &nlin, &nclin, &ldA, &ldJ, &ldR, *A, lb, ub, &confun, &ObjFuncForModel_npsolconlin, &inform, + &itera, istate, c, cJac, clambda, &result, g_dv->v, R, x_dv->v, iw, &leniw, w, &lenw); + + minpack_ps->fret = result; + + if (inform == 0) { + // The function completed successfullly + printf("\n===Ending the NPSOL constrained optimization===\n"); + } + else if (inform == 1) { + // The final iterate of x satisfies the optimality condition to eh accuracy, but has yet to converge + printf("\n===Ending the NPSOL constrained optimization, DID NOT CONVERGE===\n"); + } else { + // There was a problem + printf("\n-------NPSOL had a problem estimating ----------\n"); + printf("Error Code: %d\n\n", (int) inform); + } + + + //=== Printing out messages indicating that IMSL has bugs. + if (minpack_ps->fret > GLB_FVALMIN) + { + //NPSOL linearly-constrained optimization returns a higher obj. func. (a bug). + printf("\n----------NPSOL linearly-constrained minimization finished but with a higher objective function value!----------\n"); + printf("The improperly-returned value is %.10f while the lowest value of the objective function is %.16e.\n\n", minpack_ps->fret, GLB_FVALMIN); + fflush(stdout); + fprintf(FPTR_DEBUG, "\n----------NPSOL linearly-constrained minimization finished but with a higher objective function value!----------\n"); + fprintf(FPTR_DEBUG, "The improperly-returned value is %.16e while the lowest value of the objective function is %.16e.\n\n", minpack_ps->fret, GLB_FVALMIN); + fflush(FPTR_DEBUG); + } + + ConvertFreeParametersToQ(smodel_ps,x2_pd); + //DW's function, which takes care of the degenerate case where x2_pd points to an + // invalid place as in the constant parameter case. + ConvertFreeParametersToTheta(smodel_ps,x1_pd); //DW's function, which calls TZ's function. So essentially it's TZ's function. + + //Saved the last best results in case the IMSL quits with a bug. + CopyVector0(package_npsolconlin_ps->xsaved_dv, XNPSOL_DV); + + + //=== + DestroyVector_lf(xguess_dv); + +#endif +} + +//----------------------------------------------------------------------- +// ObjFuncForModel_npsolconlin(): +// Objective Function for linear constrained NPSOL problem +// +// mode: negative to exit NPSOL +// N: dimensions of X +// x0: parameters to maximize over +// f: value of objective function evaluated at x0 +// g: gradient vector evaluated at x0: g(j) -> df/dx_j +// nstate: can be ignored (allows for if the objective function needs +// time to initiate itself, then nstate=1 the first time it is called +//----------------------------------------------------------------------- +static void ObjFuncForModel_npsolconlin(int *mode, int *N, double *x0, double *f, double *g, int *nstate) +{ + TSdvector x0_sdv; + // + FILE *fptr_startingpoint_vec = NULL; + static int ncnt_fevals = -1; + + // printf("\n----- Entering the objective function. ------"); + // fflush(stdout); + x0_sdv.v = x0; + x0_sdv.n = *N; + x0_sdv.flag = V_DEF; + + *f = -opt_logOverallPosteriorKernal(SMODEL_PS, &x0_sdv); + if ( GLB_DISPLAY) { + printf("\nValue of objective function at the %dth evaluation: %.16e\n", ++ncnt_fevals, *f); + fflush(stdout); + } + if (*f < GLB_FVALMIN) { + //=== Resets GLB_FVALMIN at *fret_p and then prints the intermediate point to a file. + fptr_startingpoint_vec = tzFopen(filename_sp_vec_minproj,"w"); + fprintf(fptr_startingpoint_vec, "================= Output from IMSC linear constrained optimization ====================\n"); + fprintf(fptr_startingpoint_vec, "NPSOL: Value of objective miminization function at the %dth iteration: %.15f\n", ncnt_fevals, GLB_FVALMIN=*f); + fprintf(fptr_startingpoint_vec, "--------Restarting point---------\n"); + WriteVector(fptr_startingpoint_vec, &x0_sdv, " %0.16e "); + CopyVector0(XNPSOL_DV, &x0_sdv); //Saved in case the IMSL quits with a bug. + //=== Must print this results because imsl_d_min_con_gen_lin() has a bug and quits with a higher value. The printed-out results in the debug file may be used for imsl_d_min_con_nonlin() to continue. + fprintf(FPTR_DEBUG, "\nNPSOL: Value of objective miminization function at the %dth iteration: %.15f\n", ncnt_fevals, GLB_FVALMIN=*f); + fprintf(FPTR_DEBUG, "--------Restarting point---------\n"); + WriteVector(FPTR_DEBUG, &x0_sdv, " %0.16e "); + fflush(FPTR_DEBUG); + } + + // printf("\n----- Leaving the objective function. ------\n"); + // fflush(stdout); + + + // Deal with the gradient + + // printf("\n=== Entering the gradient function. ===\n"); + // fflush(stdout); + // TODO: JAKE FIX + //GLB_DISPLAY = 0; //This guarantees that the objective function printouts will not show when the gradient is computed. + gradcd_gen(g, x0, *N, ObjFuncForModel_congrad, (double *)NULL, ObjFuncForModel_congrad(x0, *N)); + + //=== Prints the intermediate gradient to a file. + // fptr_startingpoint_vec = tzFopen(filename_sp_vec_minproj,"r"); + // fprintf(fptr_startingpoint_vec, "--------Numerical gradient---------\n"); + // for (_i=0; _i<n; _i++) fprintf(fptr_startingpoint_vec, " %0.16e ", g[_i]); + // tzFclose(fptr_startingpoint_vec); + + GLB_DISPLAY = 1; + // printf("\n=== Leaving the gradient function. ===\n"); + // fflush(stdout); + + tzFclose(fptr_startingpoint_vec); + +} + +// Holder for eventual nonlinear constraints +void confun(int *mode, int *ncnln, int *n, int *ldJ, int *needc, double *x, double *c, double *cJac, int *nstate) +{ + *mode=1; +} + +//----------------------------------------------------------------------- +// Linearly-constrained IMSL minimization package. +//----------------------------------------------------------------------- +struct TSpackage_npsolconlin_tag *CreateTSpackagae_npsolconlin(const int npars_tot, const int num_lin, const int num_nlin) +{ + // n: total number of variables + // num_lin: the number of linear constraints + // num_nlin: the number of non-linear constraints + // + // General bounds again for equality lb(i) = ub(i) + // lb: n x 1 lower bound on paramters + // ub: n x 1 upper bound on paramters + // + // LINEAR BOUNDS if you want equality on an equation, set + // lin_lb(i) == lin_ub(i) + // lin_lb: (num_lin x 1) the lower bound of linear constraints + // lin_ub: (num_lin x 1) the upper bound of linear constraints + // + // NON-LINEAR BOUNDS + // nlin_lb: (num_nlin x 1) the lower bound for nonlinear constraint + // nlin_ub: (num_nlin x 1) the upper bound for nonlinear constraint + // + // A: (num_nlin x n) linear constraint matrix, with each row representing an equation + // + // + // EX: + // n = 4; + // x_1 + x_2 + x_3 >= 5 + // x_2 + 2*x_4 = 4 + // ==> + // A[0][0] = 1; + // A[0][1] = 1; + // A[0][2] = 1; + // A[1][1] = 1; + // A[1][3] = 2; + // nlin_lb[0] = 5; + // nlin_lb[1] = 4; + // nlin_ub[1] = 4; + + //=== + struct TSpackage_npsolconlin_tag *package_npsolconlin_ps = tzMalloc(1, struct TSpackage_npsolconlin_tag); + + if (npars_tot<=0) + fn_DisplayError("CreateTSpackage_npsolconlin(): make sure (1) # of equality constraints no greater than total # of constraints" + "\t and (2) number of free parameters must be greater than 0"); + package_npsolconlin_ps->n = npars_tot; + package_npsolconlin_ps->num_lin = num_lin; + package_npsolconlin_ps->num_nlin = num_nlin; + + if (num_lin<=0) + { + package_npsolconlin_ps->lin_lb = NULL; + package_npsolconlin_ps->lin_ub = NULL; + package_npsolconlin_ps->A = NULL; + } + else + { + package_npsolconlin_ps->lin_lb = CreateConstantVector_lf(num_lin, -NPSOLINF); + package_npsolconlin_ps->lin_ub = CreateConstantVector_lf(num_lin, NPSOLINF); + package_npsolconlin_ps->A = CreateConstantMatrix_lf(num_lin, npars_tot, 0); + } + if (num_nlin<=0) + { + package_npsolconlin_ps->nlin_lb = NULL; + package_npsolconlin_ps->nlin_ub = NULL; + } + else + { + package_npsolconlin_ps->nlin_lb = CreateConstantVector_lf(num_nlin, -NPSOLINF); + package_npsolconlin_ps->nlin_ub = CreateConstantVector_lf(num_nlin, NPSOLINF); + } + + // default to blank non linear constraint function + package_npsolconlin_ps->nlin_func = &confun; + + + package_npsolconlin_ps->lb = CreateConstantVector_lf(npars_tot, -NPSOLINF); + package_npsolconlin_ps->ub = CreateConstantVector_lf(npars_tot, NPSOLINF); + //- + package_npsolconlin_ps->xsaved_dv = CreateVector_lf(package_npsolconlin_ps->n); + XNPSOL_DV = CreateVector_lf(package_npsolconlin_ps->n); //Used in ObjFuncForModel_npsolconlin() to save the minimized value in case the IMSL quits with a higher value. + + package_npsolconlin_ps->crit = CRIT_NPSOLCONLIN; + package_npsolconlin_ps->itmax = ITMAX_NPSOLCONLIN; + + + return (package_npsolconlin_ps); +} +//--- +struct TSpackage_npsolconlin_tag *DestroyTSpackagae_npsolconlin(struct TSpackage_npsolconlin_tag *package_npsolconlin_ps) +{ + if (package_npsolconlin_ps) + { + DestroyVector_lf(package_npsolconlin_ps->nlin_lb); + DestroyVector_lf(package_npsolconlin_ps->nlin_ub); + DestroyVector_lf(package_npsolconlin_ps->lin_lb); + DestroyVector_lf(package_npsolconlin_ps->lin_ub); + DestroyMatrix_lf(package_npsolconlin_ps->A); + + // + DestroyVector_lf(package_npsolconlin_ps->xsaved_dv); + DestroyVector_lf(XNPSOL_DV); + + //=== + free(package_npsolconlin_ps); + return ((struct TSpackage_npsolconlin_tag *)NULL); + } + else return (package_npsolconlin_ps); +} + +static struct TStateModel_tag *SetModelGlobalForNPSOLconlin(struct TStateModel_tag *smodel_ps) +{ + //Returns the old pointer in order to preserve the previous value. + struct TStateModel_tag *tmp_ps = SMODEL_PS; + SMODEL_PS = smodel_ps; + return (tmp_ps); +} + +static void npsolconlin_SetPrintFile(char *filename) { + if (!filename) sprintf(filename_sp_vec_minproj, "outdata5npsolonlin.txt"); //Default filename. + else if (STRLEN-1 < strlen(filename)) fn_DisplayError(".../optpackage.c: the allocated length STRLEN for filename_sp_vec_minproj is too short. Must increase the string length"); + else strcpy(filename_sp_vec_minproj, filename); +} diff --git a/CFiles/optpackage.h b/CFiles/optpackage.h index 088d99c..50d187f 100644 --- a/CFiles/optpackage.h +++ b/CFiles/optpackage.h @@ -23,6 +23,12 @@ #include "dw_switch_opt.h" //DW's optimization routines for Markov-switching models. #include "cstz.h" //Used for gradcd_gen() only in the IMSL linear constrainted problem. +// ARE WE GOING TO BE USING NPSOL? +#ifdef _NPSOL_ + #define npsol npsol_ + #define npoptn npoptn_ +#endif + //-------------- Attributes for selecting optimization packages. -------------- #define MIN_DEFAULT 0 //0 or NULL: default or no minimization package. #define MIN_CSMINWEL 0x0001 //1: csminwel unconstrained minimization package. @@ -53,6 +59,15 @@ #define CRIT_IMSLCONLIN 1.0e-09 //Overall convergence criterion on the first-order conditions. #define ITMAX_IMSLCONLIN 100000 //Maximum number of iterations. + +//-------------- Minimization package: linearly-nconstrained NPSOL. -------------- +#define CRIT_NPSOLCONLIN 1.0e-09 //Overall convergence criterion on the first-order conditions. +#define ITMAX_NPSOLCONLIN 100000 //Maximum number of iterations. +#define NPSOLINF 1.0e21 +#define MAXPARAMETERS 50; +#define MAXLINCONSTRAINTS 20; + + //-------------- Minimization package: conjugate gradient method I. -------------- #define CRIT_CONGRAD1 1.0e-09 //Overall convergence criterion on the first-order conditions. #define ITMAX_CONGRAD1 100000 //Maximum number of iterations. @@ -338,6 +353,70 @@ typedef struct TSminpack_csminwel_tag { +//================ For NPSOL multivariate linearly-constrained minimizaiton package only. ================// +//----------------------------------------------------------------------- +// Using Linearly-constrained NPSOL minimization package. +// Added: 4/21/10 - Jake Smith +//----------------------------------------------------------------------- + +typedef struct TSpackage_npsolconlin_tag +{ + // n: total number of variables + // nlin: the number of linear constraints + // nclin: the number of non-linear constraints + // + // General bounds again for equality lb(i) = ub(i) + // lb: n x 1 lower bound on paramters + // ub: n x 1 upper bound on paramters + // + // LINEAR BOUNDS if you want equality on an equation, set + // nlin_lb(i) == nlin_ub(i) + // nlin_lb: (nlin x 1) the lower bound of linear constraints + // nlin_ub: (nlin x 1) the upper bound of linear constraints + // + // NON-LINEAR BOUNDS + // nclin_lb: (nclin x 1) the lower bound for nonlinear constraint + // nclin_ub: (nclin x 1) the upper bound for nonlinear constraint + // + // A: (nclin x n) linear constraint matrix, with each row representing an equation + + int n; //Total number of free parameters for the optimaization. + //For example, model free parameters + free transition matrix parameters in the regime-switching case. + int num_lin; //the number of linear constraints + int num_nlin; //the number of non-linear constraints + + TSdmatrix *A; //(nclin x n) linear constraint matrix, with each row representing an equation + + TSdvector *lb; // lb: n x 1 lower bound on paramters + TSdvector *ub; // ub: n x 1 upper bound on paramters + // + // LINEAR BOUNDS if you want equality on an equation, set + // lin_lb(i) == lin_ub(i) + TSdvector *lin_lb; // nlin_lb: (nlin x 1) the lower bound of linear constraints + TSdvector *lin_ub; // nlin_ub: (nlin x 1) the upper bound of linear constraints + // + // NON-LINEAR BOUNDS + TSdvector *nlin_lb; // nclin_lb: (nclin x 1) the lower bound for nonlinear constraint + TSdvector *nlin_ub; // nclin_ub: (nclin x 1) the upper bound for nonlinear constraint + + //=== Other output. + TSdvector *xsaved_dv; //npars_tot-by-1. Saved the parameters that give the minimal value of the objective function. + + // NONLINEAR CONSTRAINT FUNCTION + void *nlin_func; + + //=== Other inputs. + double crit; //Overall convergence criterion on the first-order conditions. + int itmax; //Maximum number of iterations. +} TSpackage_npsolconlin; +//+ +struct TSpackage_npsolconlin_tag *CreateTSpackagae_npsolconlin(const int npars_tot, const int neqs, const int ncons); +struct TSpackage_npsolconlin_tag *DestroyTSpackagae_npsolconlin(struct TSpackage_npsolconlin_tag *package_npsolconlin_ps); + +void minfinder_noblocknpsolconlin(struct TSpackage_npsolconlin_tag *package_npsolconlin_ps, struct TSminpack_tag *minpack_ps, char *filename_printout, int ntheta); +static void ObjFuncForModel_npsolconlin(int *mode, int *N, double *x, double *f, double *g, int *nstate); + + #endif -- GitLab