Mercurial > hg > camir-aes2014

annotate toolboxes/SVM-light/src/svm_hideo.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /***********************************************************************/
wolffd@0 2 /* */
wolffd@0 3 /* svm_hideo.c */
wolffd@0 4 /* */
wolffd@0 5 /* The Hildreth and D'Espo solver specialized for SVMs. */
wolffd@0 6 /* */
wolffd@0 7 /* Author: Thorsten Joachims */
wolffd@0 8 /* Date: 02.07.02 */
wolffd@0 9 /* */
wolffd@0 10 /* Copyright (c) 2002 Thorsten Joachims - All rights reserved */
wolffd@0 11 /* */
wolffd@0 12 /* This software is available for non-commercial use only. It must */
wolffd@0 13 /* not be modified and distributed without prior permission of the */
wolffd@0 14 /* author. The author is not responsible for implications from the */
wolffd@0 15 /* use of this software. */
wolffd@0 16 /* */
wolffd@0 17 /***********************************************************************/
wolffd@0 18
wolffd@0 19 # include <math.h>
wolffd@0 20 # include "svm_common.h"
wolffd@0 21
wolffd@0 22 /*
wolffd@0 23 solve the quadratic programming problem
wolffd@0 24
wolffd@0 25 minimize g0 * x + 1/2 x' * G * x
wolffd@0 26 subject to ce*x = ce0
wolffd@0 27 l <= x <= u
wolffd@0 28
wolffd@0 29 The linear constraint vector ce can only have -1/+1 as entries
wolffd@0 30 */
wolffd@0 31
wolffd@0 32 /* Common Block Declarations */
wolffd@0 33
wolffd@0 34 long verbosity;
wolffd@0 35
wolffd@0 36 # define PRIMAL_OPTIMAL 1
wolffd@0 37 # define DUAL_OPTIMAL 2
wolffd@0 38 # define MAXITER_EXCEEDED 3
wolffd@0 39 # define NAN_SOLUTION 4
wolffd@0 40 # define ONLY_ONE_VARIABLE 5
wolffd@0 41
wolffd@0 42 # define LARGEROUND 0
wolffd@0 43 # define SMALLROUND 1
wolffd@0 44
wolffd@0 45 /* /////////////////////////////////////////////////////////////// */
wolffd@0 46
wolffd@0 47 # define DEF_PRECISION 1E-5
wolffd@0 48 # define DEF_MAX_ITERATIONS 200
wolffd@0 49 # define DEF_LINDEP_SENSITIVITY 1E-8
wolffd@0 50 # define EPSILON_HIDEO 1E-20
wolffd@0 51 # define EPSILON_EQ 1E-5
wolffd@0 52
wolffd@0 53 double *optimize_qp(QP *, double *, long, double *, LEARN_PARM *);
wolffd@0 54 double *primal=0,*dual=0;
wolffd@0 55 long precision_violations=0;
wolffd@0 56 double opt_precision=DEF_PRECISION;
wolffd@0 57 long maxiter=DEF_MAX_ITERATIONS;
wolffd@0 58 double lindep_sensitivity=DEF_LINDEP_SENSITIVITY;
wolffd@0 59 double *buffer;
wolffd@0 60 long *nonoptimal;
wolffd@0 61
wolffd@0 62 long smallroundcount=0;
wolffd@0 63 long roundnumber=0;
wolffd@0 64
wolffd@0 65 /* /////////////////////////////////////////////////////////////// */
wolffd@0 66
wolffd@0 67 void *my_malloc();
wolffd@0 68
wolffd@0 69 int optimize_hildreth_despo(long,long,double,double,double,long,long,long,double,double *,
wolffd@0 70 double *,double *,double *,double *,double *,
wolffd@0 71 double *,double *,double *,long *,double *,double *);
wolffd@0 72 int solve_dual(long,long,double,double,long,double *,double *,double *,
wolffd@0 73 double *,double *,double *,double *,double *,double *,
wolffd@0 74 double *,double *,double *,double *,long);
wolffd@0 75
wolffd@0 76 void linvert_matrix(double *, long, double *, double, long *);
wolffd@0 77 void lprint_matrix(double *, long);
wolffd@0 78 void ladd_matrix(double *, long, double);
wolffd@0 79 void lcopy_matrix(double *, long, double *);
wolffd@0 80 void lswitch_rows_matrix(double *, long, long, long);
wolffd@0 81 void lswitchrk_matrix(double *, long, long, long);
wolffd@0 82
wolffd@0 83 double calculate_qp_objective(long, double *, double *, double *);
wolffd@0 84
wolffd@0 85
wolffd@0 86
wolffd@0 87 double *optimize_qp(qp,epsilon_crit,nx,threshold,learn_parm)
wolffd@0 88 QP *qp;
wolffd@0 89 double *epsilon_crit;
wolffd@0 90 long nx; /* Maximum number of variables in QP */
wolffd@0 91 double *threshold;
wolffd@0 92 LEARN_PARM *learn_parm;
wolffd@0 93 /* start the optimizer and return the optimal values */
wolffd@0 94 /* The HIDEO optimizer does not necessarily fully solve the problem. */
wolffd@0 95 /* Since it requires a strictly positive definite hessian, the solution */
wolffd@0 96 /* is restricted to a linear independent subset in case the matrix is */
wolffd@0 97 /* only semi-definite. */
wolffd@0 98 {
wolffd@0 99 long i,j;
wolffd@0 100 int result;
wolffd@0 101 double eq,progress;
wolffd@0 102
wolffd@0 103 roundnumber++;
wolffd@0 104
wolffd@0 105 if(!primal) { /* allocate memory at first call */
wolffd@0 106 primal=(double *)my_malloc(sizeof(double)*nx);
wolffd@0 107 dual=(double *)my_malloc(sizeof(double)*((nx+1)*2));
wolffd@0 108 nonoptimal=(long *)my_malloc(sizeof(long)*(nx));
wolffd@0 109 buffer=(double *)my_malloc(sizeof(double)*((nx+1)*2*(nx+1)*2+
wolffd@0 110 nx*nx+2*(nx+1)*2+2*nx+1+2*nx+
wolffd@0 111 nx+nx+nx*nx));
wolffd@0 112 (*threshold)=0;
wolffd@0 113 for(i=0;i<nx;i++) {
wolffd@0 114 primal[i]=0;
wolffd@0 115 }
wolffd@0 116 }
wolffd@0 117
wolffd@0 118 if(verbosity>=4) { /* really verbose */
wolffd@0 119 printf("\n\n");
wolffd@0 120 eq=qp->opt_ce0[0];
wolffd@0 121 for(i=0;i<qp->opt_n;i++) {
wolffd@0 122 eq+=qp->opt_xinit[i]*qp->opt_ce[i];
wolffd@0 123 printf("%f: ",qp->opt_g0[i]);
wolffd@0 124 for(j=0;j<qp->opt_n;j++) {
wolffd@0 125 printf("%f ",qp->opt_g[i*qp->opt_n+j]);
wolffd@0 126 }
wolffd@0 127 printf(": a=%.10f < %f",qp->opt_xinit[i],qp->opt_up[i]);
wolffd@0 128 printf(": y=%f\n",qp->opt_ce[i]);
wolffd@0 129 }
wolffd@0 130 if(qp->opt_m) {
wolffd@0 131 printf("EQ: %f*x0",qp->opt_ce[0]);
wolffd@0 132 for(i=1;i<qp->opt_n;i++) {
wolffd@0 133 printf(" + %f*x%ld",qp->opt_ce[i],i);
wolffd@0 134 }
wolffd@0 135 printf(" = %f\n\n",-qp->opt_ce0[0]);
wolffd@0 136 }
wolffd@0 137 }
wolffd@0 138
wolffd@0 139 result=optimize_hildreth_despo(qp->opt_n,qp->opt_m,
wolffd@0 140 opt_precision,(*epsilon_crit),
wolffd@0 141 learn_parm->epsilon_a,maxiter,
wolffd@0 142 /* (long)PRIMAL_OPTIMAL, */
wolffd@0 143 (long)0, (long)0,
wolffd@0 144 lindep_sensitivity,
wolffd@0 145 qp->opt_g,qp->opt_g0,qp->opt_ce,qp->opt_ce0,
wolffd@0 146 qp->opt_low,qp->opt_up,primal,qp->opt_xinit,
wolffd@0 147 dual,nonoptimal,buffer,&progress);
wolffd@0 148 if(verbosity>=3) {
wolffd@0 149 printf("return(%d)...",result);
wolffd@0 150 }
wolffd@0 151
wolffd@0 152 if(learn_parm->totwords < learn_parm->svm_maxqpsize) {
wolffd@0 153 /* larger working sets will be linear dependent anyway */
wolffd@0 154 learn_parm->svm_maxqpsize=maxl(learn_parm->totwords,(long)2);
wolffd@0 155 }
wolffd@0 156
wolffd@0 157 if(result == NAN_SOLUTION) {
wolffd@0 158 lindep_sensitivity*=2; /* throw out linear dependent examples more */
wolffd@0 159 /* generously */
wolffd@0 160 if(learn_parm->svm_maxqpsize>2) {
wolffd@0 161 learn_parm->svm_maxqpsize--; /* decrease size of qp-subproblems */
wolffd@0 162 }
wolffd@0 163 precision_violations++;
wolffd@0 164 }
wolffd@0 165
wolffd@0 166 /* take one round of only two variable to get unstuck */
wolffd@0 167 if((result != PRIMAL_OPTIMAL) || (!(roundnumber % 31)) || (progress <= 0)) {
wolffd@0 168
wolffd@0 169 smallroundcount++;
wolffd@0 170
wolffd@0 171 result=optimize_hildreth_despo(qp->opt_n,qp->opt_m,
wolffd@0 172 opt_precision,(*epsilon_crit),
wolffd@0 173 learn_parm->epsilon_a,(long)maxiter,
wolffd@0 174 (long)PRIMAL_OPTIMAL,(long)SMALLROUND,
wolffd@0 175 lindep_sensitivity,
wolffd@0 176 qp->opt_g,qp->opt_g0,qp->opt_ce,qp->opt_ce0,
wolffd@0 177 qp->opt_low,qp->opt_up,primal,qp->opt_xinit,
wolffd@0 178 dual,nonoptimal,buffer,&progress);
wolffd@0 179 if(verbosity>=3) {
wolffd@0 180 printf("return_srd(%d)...",result);
wolffd@0 181 }
wolffd@0 182
wolffd@0 183 if(result != PRIMAL_OPTIMAL) {
wolffd@0 184 if(result != ONLY_ONE_VARIABLE)
wolffd@0 185 precision_violations++;
wolffd@0 186 if(result == MAXITER_EXCEEDED)
wolffd@0 187 maxiter+=100;
wolffd@0 188 if(result == NAN_SOLUTION) {
wolffd@0 189 lindep_sensitivity*=2; /* throw out linear dependent examples more */
wolffd@0 190 /* generously */
wolffd@0 191 /* results not valid, so return inital values */
wolffd@0 192 for(i=0;i<qp->opt_n;i++) {
wolffd@0 193 primal[i]=qp->opt_xinit[i];
wolffd@0 194 }
wolffd@0 195 }
wolffd@0 196 }
wolffd@0 197 }
wolffd@0 198
wolffd@0 199
wolffd@0 200 if(precision_violations > 50) {
wolffd@0 201 precision_violations=0;
wolffd@0 202 (*epsilon_crit)*=10.0;
wolffd@0 203 if(verbosity>=1) {
wolffd@0 204 printf("\nWARNING: Relaxing epsilon on KT-Conditions (%f).\n",
wolffd@0 205 (*epsilon_crit));
wolffd@0 206 }
wolffd@0 207 }
wolffd@0 208
wolffd@0 209 if((qp->opt_m>0) && (result != NAN_SOLUTION) && (!isnan(dual[1]-dual[0])))
wolffd@0 210 (*threshold)=dual[1]-dual[0];
wolffd@0 211 else
wolffd@0 212 (*threshold)=0;
wolffd@0 213
wolffd@0 214 if(verbosity>=4) { /* really verbose */
wolffd@0 215 printf("\n\n");
wolffd@0 216 eq=qp->opt_ce0[0];
wolffd@0 217 for(i=0;i<qp->opt_n;i++) {
wolffd@0 218 eq+=primal[i]*qp->opt_ce[i];
wolffd@0 219 printf("%f: ",qp->opt_g0[i]);
wolffd@0 220 for(j=0;j<qp->opt_n;j++) {
wolffd@0 221 printf("%f ",qp->opt_g[i*qp->opt_n+j]);
wolffd@0 222 }
wolffd@0 223 printf(": a=%.30f",primal[i]);
wolffd@0 224 printf(": nonopti=%ld",nonoptimal[i]);
wolffd@0 225 printf(": y=%f\n",qp->opt_ce[i]);
wolffd@0 226 }
wolffd@0 227 printf("eq-constraint=%.30f\n",eq);
wolffd@0 228 printf("b=%f\n",(*threshold));
wolffd@0 229 printf(" smallroundcount=%ld ",smallroundcount);
wolffd@0 230 }
wolffd@0 231
wolffd@0 232 return(primal);
wolffd@0 233 }
wolffd@0 234
wolffd@0 235
wolffd@0 236
wolffd@0 237 int optimize_hildreth_despo(n,m,precision,epsilon_crit,epsilon_a,maxiter,goal,
wolffd@0 238 smallround,lindep_sensitivity,g,g0,ce,ce0,low,up,
wolffd@0 239 primal,init,dual,lin_dependent,buffer,progress)
wolffd@0 240 long n; /* number of variables */
wolffd@0 241 long m; /* number of linear equality constraints [0,1] */
wolffd@0 242 double precision; /* solve at least to this dual precision */
wolffd@0 243 double epsilon_crit; /* stop, if KT-Conditions approx fulfilled */
wolffd@0 244 double epsilon_a; /* precision of alphas at bounds */
wolffd@0 245 long maxiter; /* stop after this many iterations */
wolffd@0 246 long goal; /* keep going until goal fulfilled */
wolffd@0 247 long smallround; /* use only two variables of steepest descent */
wolffd@0 248 double lindep_sensitivity; /* epsilon for detecting linear dependent ex */
wolffd@0 249 double *g; /* hessian of objective */
wolffd@0 250 double *g0; /* linear part of objective */
wolffd@0 251 double *ce,*ce0; /* linear equality constraints */
wolffd@0 252 double *low,*up; /* box constraints */
wolffd@0 253 double *primal; /* primal variables */
wolffd@0 254 double *init; /* initial values of primal */
wolffd@0 255 double *dual; /* dual variables */
wolffd@0 256 long *lin_dependent;
wolffd@0 257 double *buffer;
wolffd@0 258 double *progress; /* delta in the objective function between
wolffd@0 259 before and after */
wolffd@0 260 {
wolffd@0 261 long i,j,k,from,to,n_indep,changed;
wolffd@0 262 double sum,bmin=0,bmax=0;
wolffd@0 263 double *d,*d0,*ig,*dual_old,*temp,*start;
wolffd@0 264 double *g0_new,*g_new,*ce_new,*ce0_new,*low_new,*up_new;
wolffd@0 265 double add,t;
wolffd@0 266 int result;
wolffd@0 267 double obj_before,obj_after;
wolffd@0 268 long b1,b2;
wolffd@0 269 double g0_b1,g0_b2,ce0_b;
wolffd@0 270
wolffd@0 271 g0_new=&(buffer[0]); /* claim regions of buffer */
wolffd@0 272 d=&(buffer[n]);
wolffd@0 273 d0=&(buffer[n+(n+m)*2*(n+m)*2]);
wolffd@0 274 ce_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2]);
wolffd@0 275 ce0_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n]);
wolffd@0 276 ig=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m]);
wolffd@0 277 dual_old=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n]);
wolffd@0 278 low_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2]);
wolffd@0 279 up_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n]);
wolffd@0 280 start=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n]);
wolffd@0 281 g_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n+n]);
wolffd@0 282 temp=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n+n+n*n]);
wolffd@0 283
wolffd@0 284 b1=-1;
wolffd@0 285 b2=-1;
wolffd@0 286 for(i=0;i<n;i++) { /* get variables with steepest feasible descent */
wolffd@0 287 sum=g0[i];
wolffd@0 288 for(j=0;j<n;j++)
wolffd@0 289 sum+=init[j]*g[i*n+j];
wolffd@0 290 sum=sum*ce[i];
wolffd@0 291 if(((b1==-1) || (sum<bmin))
wolffd@0 292 && (!((init[i]<=(low[i]+epsilon_a)) && (ce[i]<0.0)))
wolffd@0 293 && (!((init[i]>=( up[i]-epsilon_a)) && (ce[i]>0.0)))
wolffd@0 294 ) {
wolffd@0 295 bmin=sum;
wolffd@0 296 b1=i;
wolffd@0 297 }
wolffd@0 298 if(((b2==-1) || (sum>=bmax))
wolffd@0 299 && (!((init[i]<=(low[i]+epsilon_a)) && (ce[i]>0.0)))
wolffd@0 300 && (!((init[i]>=( up[i]-epsilon_a)) && (ce[i]<0.0)))
wolffd@0 301 ) {
wolffd@0 302 bmax=sum;
wolffd@0 303 b2=i;
wolffd@0 304 }
wolffd@0 305 }
wolffd@0 306 /* in case of unbiased hyperplane, the previous projection on */
wolffd@0 307 /* equality constraint can lead to b1 or b2 being -1. */
wolffd@0 308 if((b1 == -1) || (b2 == -1)) {
wolffd@0 309 b1=maxl(b1,b2);
wolffd@0 310 b2=maxl(b1,b2);
wolffd@0 311 }
wolffd@0 312
wolffd@0 313 for(i=0;i<n;i++) {
wolffd@0 314 start[i]=init[i];
wolffd@0 315 }
wolffd@0 316
wolffd@0 317 /* in case both example vectors are linearly dependent */
wolffd@0 318 /* WARNING: Assumes that ce[] in {-1,1} */
wolffd@0 319 add=0;
wolffd@0 320 changed=0;
wolffd@0 321 if((b1 != b2) && (m==1)) {
wolffd@0 322 for(i=0;i<n;i++) { /* fix other vectors */
wolffd@0 323 if(i==b1)
wolffd@0 324 g0_b1=g0[i];
wolffd@0 325 if(i==b2)
wolffd@0 326 g0_b2=g0[i];
wolffd@0 327 }
wolffd@0 328 ce0_b=ce0[0];
wolffd@0 329 for(i=0;i<n;i++) {
wolffd@0 330 if((i!=b1) && (i!=b2)) {
wolffd@0 331 for(j=0;j<n;j++) {
wolffd@0 332 if(j==b1)
wolffd@0 333 g0_b1+=start[i]*g[i*n+j];
wolffd@0 334 if(j==b2)
wolffd@0 335 g0_b2+=start[i]*g[i*n+j];
wolffd@0 336 }
wolffd@0 337 ce0_b-=(start[i]*ce[i]);
wolffd@0 338 }
wolffd@0 339 }
wolffd@0 340 if((g[b1*n+b2] == g[b1*n+b1]) && (g[b1*n+b2] == g[b2*n+b2])) {
wolffd@0 341 /* printf("euqal\n"); */
wolffd@0 342 if(ce[b1] == ce[b2]) {
wolffd@0 343 if(g0_b1 <= g0_b2) { /* set b1 to upper bound */
wolffd@0 344 /* printf("case +=<\n"); */
wolffd@0 345 changed=1;
wolffd@0 346 t=up[b1]-init[b1];
wolffd@0 347 if((init[b2]-low[b2]) < t) {
wolffd@0 348 t=init[b2]-low[b2];
wolffd@0 349 }
wolffd@0 350 start[b1]=init[b1]+t;
wolffd@0 351 start[b2]=init[b2]-t;
wolffd@0 352 }
wolffd@0 353 else if(g0_b1 > g0_b2) { /* set b2 to upper bound */
wolffd@0 354 /* printf("case +=>\n"); */
wolffd@0 355 changed=1;
wolffd@0 356 t=up[b2]-init[b2];
wolffd@0 357 if((init[b1]-low[b1]) < t) {
wolffd@0 358 t=init[b1]-low[b1];
wolffd@0 359 }
wolffd@0 360 start[b1]=init[b1]-t;
wolffd@0 361 start[b2]=init[b2]+t;
wolffd@0 362 }
wolffd@0 363 }
wolffd@0 364 else if(((g[b1*n+b1]>0) || (g[b2*n+b2]>0))) { /* (ce[b1] != ce[b2]) */
wolffd@0 365 /* printf("case +!\n"); */
wolffd@0 366 t=((ce[b2]/ce[b1])*g0[b1]-g0[b2]+ce0[0]*(g[b1*n+b1]*ce[b2]/ce[b1]-g[b1*n+b2]/ce[b1]))/((ce[b2]*ce[b2]/(ce[b1]*ce[b1]))*g[b1*n+b1]+g[b2*n+b2]-2*(g[b1*n+b2]*ce[b2]/ce[b1]))-init[b2];
wolffd@0 367 changed=1;
wolffd@0 368 if((up[b2]-init[b2]) < t) {
wolffd@0 369 t=up[b2]-init[b2];
wolffd@0 370 }
wolffd@0 371 if((init[b2]-low[b2]) < -t) {
wolffd@0 372 t=-(init[b2]-low[b2]);
wolffd@0 373 }
wolffd@0 374 if((up[b1]-init[b1]) < t) {
wolffd@0 375 t=(up[b1]-init[b1]);
wolffd@0 376 }
wolffd@0 377 if((init[b1]-low[b1]) < -t) {
wolffd@0 378 t=-(init[b1]-low[b1]);
wolffd@0 379 }
wolffd@0 380 start[b1]=init[b1]+t;
wolffd@0 381 start[b2]=init[b2]+t;
wolffd@0 382 }
wolffd@0 383 }
wolffd@0 384 if((-g[b1*n+b2] == g[b1*n+b1]) && (-g[b1*n+b2] == g[b2*n+b2])) {
wolffd@0 385 /* printf("diffeuqal\n"); */
wolffd@0 386 if(ce[b1] != ce[b2]) {
wolffd@0 387 if((g0_b1+g0_b2) < 0) { /* set b1 and b2 to upper bound */
wolffd@0 388 /* printf("case -!<\n"); */
wolffd@0 389 changed=1;
wolffd@0 390 t=up[b1]-init[b1];
wolffd@0 391 if((up[b2]-init[b2]) < t) {
wolffd@0 392 t=up[b2]-init[b2];
wolffd@0 393 }
wolffd@0 394 start[b1]=init[b1]+t;
wolffd@0 395 start[b2]=init[b2]+t;
wolffd@0 396 }
wolffd@0 397 else if((g0_b1+g0_b2) >= 0) { /* set b1 and b2 to lower bound */
wolffd@0 398 /* printf("case -!>\n"); */
wolffd@0 399 changed=1;
wolffd@0 400 t=init[b1]-low[b1];
wolffd@0 401 if((init[b2]-low[b2]) < t) {
wolffd@0 402 t=init[b2]-low[b2];
wolffd@0 403 }
wolffd@0 404 start[b1]=init[b1]-t;
wolffd@0 405 start[b2]=init[b2]-t;
wolffd@0 406 }
wolffd@0 407 }
wolffd@0 408 else if(((g[b1*n+b1]>0) || (g[b2*n+b2]>0))) { /* (ce[b1]==ce[b2]) */
wolffd@0 409 /* printf("case -=\n"); */
wolffd@0 410 t=((ce[b2]/ce[b1])*g0[b1]-g0[b2]+ce0[0]*(g[b1*n+b1]*ce[b2]/ce[b1]-g[b1*n+b2]/ce[b1]))/((ce[b2]*ce[b2]/(ce[b1]*ce[b1]))*g[b1*n+b1]+g[b2*n+b2]-2*(g[b1*n+b2]*ce[b2]/ce[b1]))-init[b2];
wolffd@0 411 changed=1;
wolffd@0 412 if((up[b2]-init[b2]) < t) {
wolffd@0 413 t=up[b2]-init[b2];
wolffd@0 414 }
wolffd@0 415 if((init[b2]-low[b2]) < -t) {
wolffd@0 416 t=-(init[b2]-low[b2]);
wolffd@0 417 }
wolffd@0 418 if((up[b1]-init[b1]) < -t) {
wolffd@0 419 t=-(up[b1]-init[b1]);
wolffd@0 420 }
wolffd@0 421 if((init[b1]-low[b1]) < t) {
wolffd@0 422 t=init[b1]-low[b1];
wolffd@0 423 }
wolffd@0 424 start[b1]=init[b1]-t;
wolffd@0 425 start[b2]=init[b2]+t;
wolffd@0 426 }
wolffd@0 427 }
wolffd@0 428 }
wolffd@0 429 /* if we have a biased hyperplane, then adding a constant to the */
wolffd@0 430 /* hessian does not change the solution. So that is done for examples */
wolffd@0 431 /* with zero diagonal entry, since HIDEO cannot handle them. */
wolffd@0 432 if((m>0)
wolffd@0 433 && ((fabs(g[b1*n+b1]) < lindep_sensitivity)
wolffd@0 434 || (fabs(g[b2*n+b2]) < lindep_sensitivity))) {
wolffd@0 435 /* printf("Case 0\n"); */
wolffd@0 436 add+=0.093274;
wolffd@0 437 }
wolffd@0 438 /* in case both examples are linear dependent */
wolffd@0 439 else if((m>0)
wolffd@0 440 && (g[b1*n+b2] != 0 && g[b2*n+b2] != 0)
wolffd@0 441 && (fabs(g[b1*n+b1]/g[b1*n+b2] - g[b1*n+b2]/g[b2*n+b2])
wolffd@0 442 < lindep_sensitivity)) {
wolffd@0 443 /* printf("Case lindep\n"); */
wolffd@0 444 add+=0.078274;
wolffd@0 445 }
wolffd@0 446
wolffd@0 447 /* special case for zero diagonal entry on unbiased hyperplane */
wolffd@0 448 if((m==0) && (b1>=0)) {
wolffd@0 449 if(fabs(g[b1*n+b1]) < lindep_sensitivity) {
wolffd@0 450 /* printf("Case 0b1\n"); */
wolffd@0 451 for(i=0;i<n;i++) { /* fix other vectors */
wolffd@0 452 if(i==b1)
wolffd@0 453 g0_b1=g0[i];
wolffd@0 454 }
wolffd@0 455 for(i=0;i<n;i++) {
wolffd@0 456 if(i!=b1) {
wolffd@0 457 for(j=0;j<n;j++) {
wolffd@0 458 if(j==b1)
wolffd@0 459 g0_b1+=start[i]*g[i*n+j];
wolffd@0 460 }
wolffd@0 461 }
wolffd@0 462 }
wolffd@0 463 if(g0_b1<0)
wolffd@0 464 start[b1]=up[b1];
wolffd@0 465 if(g0_b1>=0)
wolffd@0 466 start[b1]=low[b1];
wolffd@0 467 }
wolffd@0 468 }
wolffd@0 469 if((m==0) && (b2>=0)) {
wolffd@0 470 if(fabs(g[b2*n+b2]) < lindep_sensitivity) {
wolffd@0 471 /* printf("Case 0b2\n"); */
wolffd@0 472 for(i=0;i<n;i++) { /* fix other vectors */
wolffd@0 473 if(i==b2)
wolffd@0 474 g0_b2=g0[i];
wolffd@0 475 }
wolffd@0 476 for(i=0;i<n;i++) {
wolffd@0 477 if(i!=b2) {
wolffd@0 478 for(j=0;j<n;j++) {
wolffd@0 479 if(j==b2)
wolffd@0 480 g0_b2+=start[i]*g[i*n+j];
wolffd@0 481 }
wolffd@0 482 }
wolffd@0 483 }
wolffd@0 484 if(g0_b2<0)
wolffd@0 485 start[b2]=up[b2];
wolffd@0 486 if(g0_b2>=0)
wolffd@0 487 start[b2]=low[b2];
wolffd@0 488 }
wolffd@0 489 }
wolffd@0 490
wolffd@0 491 /* printf("b1=%ld,b2=%ld\n",b1,b2); */
wolffd@0 492
wolffd@0 493 lcopy_matrix(g,n,d);
wolffd@0 494 if((m==1) && (add>0.0)) {
wolffd@0 495 for(j=0;j<n;j++) {
wolffd@0 496 for(k=0;k<n;k++) {
wolffd@0 497 d[j*n+k]+=add*ce[j]*ce[k];
wolffd@0 498 }
wolffd@0 499 }
wolffd@0 500 }
wolffd@0 501 else {
wolffd@0 502 add=0.0;
wolffd@0 503 }
wolffd@0 504
wolffd@0 505 if(n>2) { /* switch, so that variables are better mixed */
wolffd@0 506 lswitchrk_matrix(d,n,b1,(long)0);
wolffd@0 507 if(b2 == 0)
wolffd@0 508 lswitchrk_matrix(d,n,b1,(long)1);
wolffd@0 509 else
wolffd@0 510 lswitchrk_matrix(d,n,b2,(long)1);
wolffd@0 511 }
wolffd@0 512 if(smallround == SMALLROUND) {
wolffd@0 513 for(i=2;i<n;i++) {
wolffd@0 514 lin_dependent[i]=1;
wolffd@0 515 }
wolffd@0 516 if(m>0) { /* for biased hyperplane, pick two variables */
wolffd@0 517 lin_dependent[0]=0;
wolffd@0 518 lin_dependent[1]=0;
wolffd@0 519 }
wolffd@0 520 else { /* for unbiased hyperplane, pick only one variable */
wolffd@0 521 lin_dependent[0]=smallroundcount % 2;
wolffd@0 522 lin_dependent[1]=(smallroundcount+1) % 2;
wolffd@0 523 }
wolffd@0 524 }
wolffd@0 525 else {
wolffd@0 526 for(i=0;i<n;i++) {
wolffd@0 527 lin_dependent[i]=0;
wolffd@0 528 }
wolffd@0 529 }
wolffd@0 530 linvert_matrix(d,n,ig,lindep_sensitivity,lin_dependent);
wolffd@0 531 if(n>2) { /* now switch back */
wolffd@0 532 if(b2 == 0) {
wolffd@0 533 lswitchrk_matrix(ig,n,b1,(long)1);
wolffd@0 534 i=lin_dependent[1];
wolffd@0 535 lin_dependent[1]=lin_dependent[b1];
wolffd@0 536 lin_dependent[b1]=i;
wolffd@0 537 }
wolffd@0 538 else {
wolffd@0 539 lswitchrk_matrix(ig,n,b2,(long)1);
wolffd@0 540 i=lin_dependent[1];
wolffd@0 541 lin_dependent[1]=lin_dependent[b2];
wolffd@0 542 lin_dependent[b2]=i;
wolffd@0 543 }
wolffd@0 544 lswitchrk_matrix(ig,n,b1,(long)0);
wolffd@0 545 i=lin_dependent[0];
wolffd@0 546 lin_dependent[0]=lin_dependent[b1];
wolffd@0 547 lin_dependent[b1]=i;
wolffd@0 548 }
wolffd@0 549 /* lprint_matrix(d,n); */
wolffd@0 550 /* lprint_matrix(ig,n); */
wolffd@0 551
wolffd@0 552 lcopy_matrix(g,n,g_new); /* restore g_new matrix */
wolffd@0 553 if(add>0)
wolffd@0 554 for(j=0;j<n;j++) {
wolffd@0 555 for(k=0;k<n;k++) {
wolffd@0 556 g_new[j*n+k]+=add*ce[j]*ce[k];
wolffd@0 557 }
wolffd@0 558 }
wolffd@0 559
wolffd@0 560 for(i=0;i<n;i++) { /* fix linear dependent vectors */
wolffd@0 561 g0_new[i]=g0[i]+add*ce0[0]*ce[i];
wolffd@0 562 }
wolffd@0 563 if(m>0) ce0_new[0]=-ce0[0];
wolffd@0 564 for(i=0;i<n;i++) { /* fix linear dependent vectors */
wolffd@0 565 if(lin_dependent[i]) {
wolffd@0 566 for(j=0;j<n;j++) {
wolffd@0 567 if(!lin_dependent[j]) {
wolffd@0 568 g0_new[j]+=start[i]*g_new[i*n+j];
wolffd@0 569 }
wolffd@0 570 }
wolffd@0 571 if(m>0) ce0_new[0]-=(start[i]*ce[i]);
wolffd@0 572 }
wolffd@0 573 }
wolffd@0 574 from=0; /* remove linear dependent vectors */
wolffd@0 575 to=0;
wolffd@0 576 n_indep=0;
wolffd@0 577 for(i=0;i<n;i++) {
wolffd@0 578 if(!lin_dependent[i]) {
wolffd@0 579 g0_new[n_indep]=g0_new[i];
wolffd@0 580 ce_new[n_indep]=ce[i];
wolffd@0 581 low_new[n_indep]=low[i];
wolffd@0 582 up_new[n_indep]=up[i];
wolffd@0 583 primal[n_indep]=start[i];
wolffd@0 584 n_indep++;
wolffd@0 585 }
wolffd@0 586 for(j=0;j<n;j++) {
wolffd@0 587 if((!lin_dependent[i]) && (!lin_dependent[j])) {
wolffd@0 588 ig[to]=ig[from];
wolffd@0 589 g_new[to]=g_new[from];
wolffd@0 590 to++;
wolffd@0 591 }
wolffd@0 592 from++;
wolffd@0 593 }
wolffd@0 594 }
wolffd@0 595
wolffd@0 596 if(verbosity>=3) {
wolffd@0 597 printf("real_qp_size(%ld)...",n_indep);
wolffd@0 598 }
wolffd@0 599
wolffd@0 600 /* cannot optimize with only one variable */
wolffd@0 601 if((n_indep<=1) && (m>0) && (!changed)) {
wolffd@0 602 for(i=n-1;i>=0;i--) {
wolffd@0 603 primal[i]=init[i];
wolffd@0 604 }
wolffd@0 605 return((int)ONLY_ONE_VARIABLE);
wolffd@0 606 }
wolffd@0 607
wolffd@0 608 if((!changed) || (n_indep>1)) {
wolffd@0 609 result=solve_dual(n_indep,m,precision,epsilon_crit,maxiter,g_new,g0_new,
wolffd@0 610 ce_new,ce0_new,low_new,up_new,primal,d,d0,ig,
wolffd@0 611 dual,dual_old,temp,goal);
wolffd@0 612 }
wolffd@0 613 else {
wolffd@0 614 result=PRIMAL_OPTIMAL;
wolffd@0 615 }
wolffd@0 616
wolffd@0 617 j=n_indep;
wolffd@0 618 for(i=n-1;i>=0;i--) {
wolffd@0 619 if(!lin_dependent[i]) {
wolffd@0 620 j--;
wolffd@0 621 primal[i]=primal[j];
wolffd@0 622 }
wolffd@0 623 else {
wolffd@0 624 primal[i]=start[i]; /* leave as is */
wolffd@0 625 }
wolffd@0 626 temp[i]=primal[i];
wolffd@0 627 }
wolffd@0 628
wolffd@0 629 obj_before=calculate_qp_objective(n,g,g0,init);
wolffd@0 630 obj_after=calculate_qp_objective(n,g,g0,primal);
wolffd@0 631 (*progress)=obj_before-obj_after;
wolffd@0 632 if(verbosity>=3) {
wolffd@0 633 printf("before(%.30f)...after(%.30f)...result_sd(%d)...",
wolffd@0 634 obj_before,obj_after,result);
wolffd@0 635 }
wolffd@0 636
wolffd@0 637 return((int)result);
wolffd@0 638 }
wolffd@0 639
wolffd@0 640
wolffd@0 641 int solve_dual(n,m,precision,epsilon_crit,maxiter,g,g0,ce,ce0,low,up,primal,
wolffd@0 642 d,d0,ig,dual,dual_old,temp,goal)
wolffd@0 643 /* Solves the dual using the method of Hildreth and D'Espo. */
wolffd@0 644 /* Can only handle problems with zero or exactly one */
wolffd@0 645 /* equality constraints. */
wolffd@0 646
wolffd@0 647 long n; /* number of variables */
wolffd@0 648 long m; /* number of linear equality constraints */
wolffd@0 649 double precision; /* solve at least to this dual precision */
wolffd@0 650 double epsilon_crit; /* stop, if KT-Conditions approx fulfilled */
wolffd@0 651 long maxiter; /* stop after that many iterations */
wolffd@0 652 double *g;
wolffd@0 653 double *g0; /* linear part of objective */
wolffd@0 654 double *ce,*ce0; /* linear equality constraints */
wolffd@0 655 double *low,*up; /* box constraints */
wolffd@0 656 double *primal; /* variables (with initial values) */
wolffd@0 657 double *d,*d0,*ig,*dual,*dual_old,*temp; /* buffer */
wolffd@0 658 long goal;
wolffd@0 659 {
wolffd@0 660 long i,j,k,iter;
wolffd@0 661 double sum,w,maxviol,viol,temp1,temp2,isnantest;
wolffd@0 662 double model_b,dist;
wolffd@0 663 long retrain,maxfaktor,primal_optimal=0,at_bound,scalemaxiter;
wolffd@0 664 double epsilon_a=1E-15,epsilon_hideo;
wolffd@0 665 double eq;
wolffd@0 666
wolffd@0 667 if((m<0) || (m>1))
wolffd@0 668 perror("SOLVE DUAL: inappropriate number of eq-constrains!");
wolffd@0 669
wolffd@0 670 /*
wolffd@0 671 printf("\n");
wolffd@0 672 for(i=0;i<n;i++) {
wolffd@0 673 printf("%f: ",g0[i]);
wolffd@0 674 for(j=0;j<n;j++) {
wolffd@0 675 printf("%f ",g[i*n+j]);
wolffd@0 676 }
wolffd@0 677 printf(": a=%.30f",primal[i]);
wolffd@0 678 printf(": y=%f\n",ce[i]);
wolffd@0 679 }
wolffd@0 680 */
wolffd@0 681
wolffd@0 682 for(i=0;i<2*(n+m);i++) {
wolffd@0 683 dual[i]=0;
wolffd@0 684 dual_old[i]=0;
wolffd@0 685 }
wolffd@0 686 for(i=0;i<n;i++) {
wolffd@0 687 for(j=0;j<n;j++) { /* dual hessian for box constraints */
wolffd@0 688 d[i*2*(n+m)+j]=ig[i*n+j];
wolffd@0 689 d[(i+n)*2*(n+m)+j]=-ig[i*n+j];
wolffd@0 690 d[i*2*(n+m)+j+n]=-ig[i*n+j];
wolffd@0 691 d[(i+n)*2*(n+m)+j+n]=ig[i*n+j];
wolffd@0 692 }
wolffd@0 693 if(m>0) {
wolffd@0 694 sum=0; /* dual hessian for eq constraints */
wolffd@0 695 for(j=0;j<n;j++) {
wolffd@0 696 sum+=(ce[j]*ig[i*n+j]);
wolffd@0 697 }
wolffd@0 698 d[i*2*(n+m)+2*n]=sum;
wolffd@0 699 d[i*2*(n+m)+2*n+1]=-sum;
wolffd@0 700 d[(n+i)*2*(n+m)+2*n]=-sum;
wolffd@0 701 d[(n+i)*2*(n+m)+2*n+1]=sum;
wolffd@0 702 d[(n+n)*2*(n+m)+i]=sum;
wolffd@0 703 d[(n+n+1)*2*(n+m)+i]=-sum;
wolffd@0 704 d[(n+n)*2*(n+m)+(n+i)]=-sum;
wolffd@0 705 d[(n+n+1)*2*(n+m)+(n+i)]=sum;
wolffd@0 706
wolffd@0 707 sum=0;
wolffd@0 708 for(j=0;j<n;j++) {
wolffd@0 709 for(k=0;k<n;k++) {
wolffd@0 710 sum+=(ce[k]*ce[j]*ig[j*n+k]);
wolffd@0 711 }
wolffd@0 712 }
wolffd@0 713 d[(n+n)*2*(n+m)+2*n]=sum;
wolffd@0 714 d[(n+n)*2*(n+m)+2*n+1]=-sum;
wolffd@0 715 d[(n+n+1)*2*(n+m)+2*n]=-sum;
wolffd@0 716 d[(n+n+1)*2*(n+m)+2*n+1]=sum;
wolffd@0 717 }
wolffd@0 718 }
wolffd@0 719
wolffd@0 720 for(i=0;i<n;i++) { /* dual linear component for the box constraints */
wolffd@0 721 w=0;
wolffd@0 722 for(j=0;j<n;j++) {
wolffd@0 723 w+=(ig[i*n+j]*g0[j]);
wolffd@0 724 }
wolffd@0 725 d0[i]=up[i]+w;
wolffd@0 726 d0[i+n]=-low[i]-w;
wolffd@0 727 }
wolffd@0 728
wolffd@0 729 if(m>0) {
wolffd@0 730 sum=0; /* dual linear component for eq constraints */
wolffd@0 731 for(j=0;j<n;j++) {
wolffd@0 732 for(k=0;k<n;k++) {
wolffd@0 733 sum+=(ce[k]*ig[k*n+j]*g0[j]);
wolffd@0 734 }
wolffd@0 735 }
wolffd@0 736 d0[2*n]=ce0[0]+sum;
wolffd@0 737 d0[2*n+1]=-ce0[0]-sum;
wolffd@0 738 }
wolffd@0 739
wolffd@0 740 maxviol=999999;
wolffd@0 741 iter=0;
wolffd@0 742 retrain=1;
wolffd@0 743 maxfaktor=1;
wolffd@0 744 scalemaxiter=maxiter/5;
wolffd@0 745 while((retrain) && (maxviol > 0) && (iter < (scalemaxiter*maxfaktor))) {
wolffd@0 746 iter++;
wolffd@0 747
wolffd@0 748 while((maxviol > precision) && (iter < (scalemaxiter*maxfaktor))) {
wolffd@0 749 iter++;
wolffd@0 750 maxviol=0;
wolffd@0 751 for(i=0;i<2*(n+m);i++) {
wolffd@0 752 sum=d0[i];
wolffd@0 753 for(j=0;j<2*(n+m);j++) {
wolffd@0 754 sum+=d[i*2*(n+m)+j]*dual_old[j];
wolffd@0 755 }
wolffd@0 756 sum-=d[i*2*(n+m)+i]*dual_old[i];
wolffd@0 757 dual[i]=-sum/d[i*2*(n+m)+i];
wolffd@0 758 if(dual[i]<0) dual[i]=0;
wolffd@0 759
wolffd@0 760 viol=fabs(dual[i]-dual_old[i]);
wolffd@0 761 if(viol>maxviol)
wolffd@0 762 maxviol=viol;
wolffd@0 763 dual_old[i]=dual[i];
wolffd@0 764 }
wolffd@0 765 /*
wolffd@0 766 printf("%d) maxviol=%20f precision=%f\n",iter,maxviol,precision);
wolffd@0 767 */
wolffd@0 768 }
wolffd@0 769
wolffd@0 770 if(m>0) {
wolffd@0 771 for(i=0;i<n;i++) {
wolffd@0 772 temp[i]=dual[i]-dual[i+n]+ce[i]*(dual[n+n]-dual[n+n+1])+g0[i];
wolffd@0 773 }
wolffd@0 774 }
wolffd@0 775 else {
wolffd@0 776 for(i=0;i<n;i++) {
wolffd@0 777 temp[i]=dual[i]-dual[i+n]+g0[i];
wolffd@0 778 }
wolffd@0 779 }
wolffd@0 780 for(i=0;i<n;i++) {
wolffd@0 781 primal[i]=0; /* calc value of primal variables */
wolffd@0 782 for(j=0;j<n;j++) {
wolffd@0 783 primal[i]+=ig[i*n+j]*temp[j];
wolffd@0 784 }
wolffd@0 785 primal[i]*=-1.0;
wolffd@0 786 if(primal[i]<=(low[i])) { /* clip conservatively */
wolffd@0 787 primal[i]=low[i];
wolffd@0 788 }
wolffd@0 789 else if(primal[i]>=(up[i])) {
wolffd@0 790 primal[i]=up[i];
wolffd@0 791 }
wolffd@0 792 }
wolffd@0 793
wolffd@0 794 if(m>0)
wolffd@0 795 model_b=dual[n+n+1]-dual[n+n];
wolffd@0 796 else
wolffd@0 797 model_b=0;
wolffd@0 798
wolffd@0 799 epsilon_hideo=EPSILON_HIDEO;
wolffd@0 800 for(i=0;i<n;i++) { /* check precision of alphas */
wolffd@0 801 dist=-model_b*ce[i];
wolffd@0 802 dist+=(g0[i]+1.0);
wolffd@0 803 for(j=0;j<i;j++) {
wolffd@0 804 dist+=(primal[j]*g[j*n+i]);
wolffd@0 805 }
wolffd@0 806 for(j=i;j<n;j++) {
wolffd@0 807 dist+=(primal[j]*g[i*n+j]);
wolffd@0 808 }
wolffd@0 809 if((primal[i]<(up[i]-epsilon_hideo)) && (dist < (1.0-epsilon_crit))) {
wolffd@0 810 epsilon_hideo=(up[i]-primal[i])*2.0;
wolffd@0 811 }
wolffd@0 812 else if((primal[i]>(low[i]+epsilon_hideo)) &&(dist>(1.0+epsilon_crit))) {
wolffd@0 813 epsilon_hideo=(primal[i]-low[i])*2.0;
wolffd@0 814 }
wolffd@0 815 }
wolffd@0 816 /* printf("\nEPSILON_HIDEO=%.30f\n",epsilon_hideo); */
wolffd@0 817
wolffd@0 818 for(i=0;i<n;i++) { /* clip alphas to bounds */
wolffd@0 819 if(primal[i]<=(low[i]+epsilon_hideo)) {
wolffd@0 820 primal[i]=low[i];
wolffd@0 821 }
wolffd@0 822 else if(primal[i]>=(up[i]-epsilon_hideo)) {
wolffd@0 823 primal[i]=up[i];
wolffd@0 824 }
wolffd@0 825 }
wolffd@0 826
wolffd@0 827 retrain=0;
wolffd@0 828 primal_optimal=1;
wolffd@0 829 at_bound=0;
wolffd@0 830 for(i=0;(i<n);i++) { /* check primal KT-Conditions */
wolffd@0 831 dist=-model_b*ce[i];
wolffd@0 832 dist+=(g0[i]+1.0);
wolffd@0 833 for(j=0;j<i;j++) {
wolffd@0 834 dist+=(primal[j]*g[j*n+i]);
wolffd@0 835 }
wolffd@0 836 for(j=i;j<n;j++) {
wolffd@0 837 dist+=(primal[j]*g[i*n+j]);
wolffd@0 838 }
wolffd@0 839 if((primal[i]<(up[i]-epsilon_a)) && (dist < (1.0-epsilon_crit))) {
wolffd@0 840 retrain=1;
wolffd@0 841 primal_optimal=0;
wolffd@0 842 }
wolffd@0 843 else if((primal[i]>(low[i]+epsilon_a)) && (dist > (1.0+epsilon_crit))) {
wolffd@0 844 retrain=1;
wolffd@0 845 primal_optimal=0;
wolffd@0 846 }
wolffd@0 847 if((primal[i]<=(low[i]+epsilon_a)) || (primal[i]>=(up[i]-epsilon_a))) {
wolffd@0 848 at_bound++;
wolffd@0 849 }
wolffd@0 850 /* printf("HIDEOtemp: a[%ld]=%.30f, dist=%.6f, b=%f, at_bound=%ld\n",i,primal[i],dist,model_b,at_bound); */
wolffd@0 851 }
wolffd@0 852 if(m>0) {
wolffd@0 853 eq=-ce0[0]; /* check precision of eq-constraint */
wolffd@0 854 for(i=0;i<n;i++) {
wolffd@0 855 eq+=(ce[i]*primal[i]);
wolffd@0 856 }
wolffd@0 857 if((EPSILON_EQ < fabs(eq))
wolffd@0 858 /*
wolffd@0 859 && !((goal==PRIMAL_OPTIMAL)
wolffd@0 860 && (at_bound==n)) */
wolffd@0 861 ) {
wolffd@0 862 retrain=1;
wolffd@0 863 primal_optimal=0;
wolffd@0 864 }
wolffd@0 865 /* printf("\n eq=%.30f ce0=%f at-bound=%ld\n",eq,ce0[0],at_bound); */
wolffd@0 866 }
wolffd@0 867
wolffd@0 868 if(retrain) {
wolffd@0 869 precision/=10;
wolffd@0 870 if(((goal == PRIMAL_OPTIMAL) && (maxfaktor < 50000))
wolffd@0 871 || (maxfaktor < 5)) {
wolffd@0 872 maxfaktor++;
wolffd@0 873 }
wolffd@0 874 }
wolffd@0 875 }
wolffd@0 876
wolffd@0 877 if(!primal_optimal) {
wolffd@0 878 for(i=0;i<n;i++) {
wolffd@0 879 primal[i]=0; /* calc value of primal variables */
wolffd@0 880 for(j=0;j<n;j++) {
wolffd@0 881 primal[i]+=ig[i*n+j]*temp[j];
wolffd@0 882 }
wolffd@0 883 primal[i]*=-1.0;
wolffd@0 884 if(primal[i]<=(low[i]+epsilon_a)) { /* clip conservatively */
wolffd@0 885 primal[i]=low[i];
wolffd@0 886 }
wolffd@0 887 else if(primal[i]>=(up[i]-epsilon_a)) {
wolffd@0 888 primal[i]=up[i];
wolffd@0 889 }
wolffd@0 890 }
wolffd@0 891 }
wolffd@0 892
wolffd@0 893 isnantest=0;
wolffd@0 894 for(i=0;i<n;i++) { /* check for isnan */
wolffd@0 895 isnantest+=primal[i];
wolffd@0 896 }
wolffd@0 897
wolffd@0 898 if(m>0) {
wolffd@0 899 temp1=dual[n+n+1]; /* copy the dual variables for the eq */
wolffd@0 900 temp2=dual[n+n]; /* constraints to a handier location */
wolffd@0 901 for(i=n+n+1;i>=2;i--) {
wolffd@0 902 dual[i]=dual[i-2];
wolffd@0 903 }
wolffd@0 904 dual[0]=temp2;
wolffd@0 905 dual[1]=temp1;
wolffd@0 906 isnantest+=temp1+temp2;
wolffd@0 907 }
wolffd@0 908
wolffd@0 909 if(isnan(isnantest)) {
wolffd@0 910 return((int)NAN_SOLUTION);
wolffd@0 911 }
wolffd@0 912 else if(primal_optimal) {
wolffd@0 913 return((int)PRIMAL_OPTIMAL);
wolffd@0 914 }
wolffd@0 915 else if(maxviol == 0.0) {
wolffd@0 916 return((int)DUAL_OPTIMAL);
wolffd@0 917 }
wolffd@0 918 else {
wolffd@0 919 return((int)MAXITER_EXCEEDED);
wolffd@0 920 }
wolffd@0 921 }
wolffd@0 922
wolffd@0 923
wolffd@0 924 void linvert_matrix(matrix,depth,inverse,lindep_sensitivity,lin_dependent)
wolffd@0 925 double *matrix;
wolffd@0 926 long depth;
wolffd@0 927 double *inverse,lindep_sensitivity;
wolffd@0 928 long *lin_dependent; /* indicates the active parts of matrix on
wolffd@0 929 input and output*/
wolffd@0 930 {
wolffd@0 931 long i,j,k;
wolffd@0 932 double factor;
wolffd@0 933
wolffd@0 934 for(i=0;i<depth;i++) {
wolffd@0 935 /* lin_dependent[i]=0; */
wolffd@0 936 for(j=0;j<depth;j++) {
wolffd@0 937 inverse[i*depth+j]=0.0;
wolffd@0 938 }
wolffd@0 939 inverse[i*depth+i]=1.0;
wolffd@0 940 }
wolffd@0 941 for(i=0;i<depth;i++) {
wolffd@0 942 if(lin_dependent[i] || (fabs(matrix[i*depth+i])<lindep_sensitivity)) {
wolffd@0 943 lin_dependent[i]=1;
wolffd@0 944 }
wolffd@0 945 else {
wolffd@0 946 for(j=i+1;j<depth;j++) {
wolffd@0 947 factor=matrix[j*depth+i]/matrix[i*depth+i];
wolffd@0 948 for(k=i;k<depth;k++) {
wolffd@0 949 matrix[j*depth+k]-=(factor*matrix[i*depth+k]);
wolffd@0 950 }
wolffd@0 951 for(k=0;k<depth;k++) {
wolffd@0 952 inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
wolffd@0 953 }
wolffd@0 954 }
wolffd@0 955 }
wolffd@0 956 }
wolffd@0 957 for(i=depth-1;i>=0;i--) {
wolffd@0 958 if(!lin_dependent[i]) {
wolffd@0 959 factor=1/matrix[i*depth+i];
wolffd@0 960 for(k=0;k<depth;k++) {
wolffd@0 961 inverse[i*depth+k]*=factor;
wolffd@0 962 }
wolffd@0 963 matrix[i*depth+i]=1;
wolffd@0 964 for(j=i-1;j>=0;j--) {
wolffd@0 965 factor=matrix[j*depth+i];
wolffd@0 966 matrix[j*depth+i]=0;
wolffd@0 967 for(k=0;k<depth;k++) {
wolffd@0 968 inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
wolffd@0 969 }
wolffd@0 970 }
wolffd@0 971 }
wolffd@0 972 }
wolffd@0 973 }
wolffd@0 974
wolffd@0 975 void lprint_matrix(matrix,depth)
wolffd@0 976 double *matrix;
wolffd@0 977 long depth;
wolffd@0 978 {
wolffd@0 979 long i,j;
wolffd@0 980 for(i=0;i<depth;i++) {
wolffd@0 981 for(j=0;j<depth;j++) {
wolffd@0 982 printf("%5.2f ",(double)(matrix[i*depth+j]));
wolffd@0 983 }
wolffd@0 984 printf("\n");
wolffd@0 985 }
wolffd@0 986 printf("\n");
wolffd@0 987 }
wolffd@0 988
wolffd@0 989 void ladd_matrix(matrix,depth,scalar)
wolffd@0 990 double *matrix;
wolffd@0 991 long depth;
wolffd@0 992 double scalar;
wolffd@0 993 {
wolffd@0 994 long i,j;
wolffd@0 995 for(i=0;i<depth;i++) {
wolffd@0 996 for(j=0;j<depth;j++) {
wolffd@0 997 matrix[i*depth+j]+=scalar;
wolffd@0 998 }
wolffd@0 999 }
wolffd@0 1000 }
wolffd@0 1001
wolffd@0 1002 void lcopy_matrix(matrix,depth,matrix2)
wolffd@0 1003 double *matrix;
wolffd@0 1004 long depth;
wolffd@0 1005 double *matrix2;
wolffd@0 1006 {
wolffd@0 1007 long i;
wolffd@0 1008
wolffd@0 1009 for(i=0;i<(depth)*(depth);i++) {
wolffd@0 1010 matrix2[i]=matrix[i];
wolffd@0 1011 }
wolffd@0 1012 }
wolffd@0 1013
wolffd@0 1014 void lswitch_rows_matrix(matrix,depth,r1,r2)
wolffd@0 1015 double *matrix;
wolffd@0 1016 long depth,r1,r2;
wolffd@0 1017 {
wolffd@0 1018 long i;
wolffd@0 1019 double temp;
wolffd@0 1020
wolffd@0 1021 for(i=0;i<depth;i++) {
wolffd@0 1022 temp=matrix[r1*depth+i];
wolffd@0 1023 matrix[r1*depth+i]=matrix[r2*depth+i];
wolffd@0 1024 matrix[r2*depth+i]=temp;
wolffd@0 1025 }
wolffd@0 1026 }
wolffd@0 1027
wolffd@0 1028 void lswitchrk_matrix(matrix,depth,rk1,rk2)
wolffd@0 1029 double *matrix;
wolffd@0 1030 long depth,rk1,rk2;
wolffd@0 1031 {
wolffd@0 1032 long i;
wolffd@0 1033 double temp;
wolffd@0 1034
wolffd@0 1035 for(i=0;i<depth;i++) {
wolffd@0 1036 temp=matrix[rk1*depth+i];
wolffd@0 1037 matrix[rk1*depth+i]=matrix[rk2*depth+i];
wolffd@0 1038 matrix[rk2*depth+i]=temp;
wolffd@0 1039 }
wolffd@0 1040 for(i=0;i<depth;i++) {
wolffd@0 1041 temp=matrix[i*depth+rk1];
wolffd@0 1042 matrix[i*depth+rk1]=matrix[i*depth+rk2];
wolffd@0 1043 matrix[i*depth+rk2]=temp;
wolffd@0 1044 }
wolffd@0 1045 }
wolffd@0 1046
wolffd@0 1047 double calculate_qp_objective(opt_n,opt_g,opt_g0,alpha)
wolffd@0 1048 long opt_n;
wolffd@0 1049 double *opt_g,*opt_g0,*alpha;
wolffd@0 1050 {
wolffd@0 1051 double obj;
wolffd@0 1052 long i,j;
wolffd@0 1053 obj=0; /* calculate objective */
wolffd@0 1054 for(i=0;i<opt_n;i++) {
wolffd@0 1055 obj+=(opt_g0[i]*alpha[i]);
wolffd@0 1056 obj+=(0.5*alpha[i]*alpha[i]*opt_g[i*opt_n+i]);
wolffd@0 1057 for(j=0;j<i;j++) {
wolffd@0 1058 obj+=(alpha[j]*alpha[i]*opt_g[j*opt_n+i]);
wolffd@0 1059 }
wolffd@0 1060 }
wolffd@0 1061 return(obj);
wolffd@0 1062 }