share/simplex/rtest_simplex.mac

   1 (kill(all), load(simplex), 0);
   2 0;
   3
   4 linear_program(matrix([1,1,-1,0],[2,-3,0,-1],[4,-5,0,0]), [1,1,6], [1,-2,0,0]);
   5 [[13/2, 4, 19/2, 0], -3/2];
   6
   7 minimize_lp(x, [y>=x-1, y>=-x-1, y<=x+1, y<=1-x, y=x/2]);
   8 [-(2/3), [x=-(2/3), y=-(1/3)]];
   9
  10 maximize_lp(x+y, [y<=-x/2+3, y<=-x+4], [x, y]);
  11 [4, [x = 2, y = 2]];
  12
  13 maximize_lp(x+y, [y<=-x/2+3, y<=-x+4], all);
  14 [4, [x = 2, y = 2]];
  15
  16 maximize_lp(x+y, [y<=-x/2+3, y<=-x+4]), nonnegative_lp=true;
  17 [4, [x = 2, y = 2]];
  18
  19 maximize_lp(10*x1 - 57*x2 - 9*x3 - 24*x4,
  20   [1/2*x1 - 11/2*x2 - 5/2*x3 + 9*x4 < 0,
  21    1/2*x1 - 3+2*x2  - 1/2*x3 + x4   < 0,
  22    x1 < 1],
  23   all)$
  24 [41/5,[x4=0,x3=1/5,x2=0,x1=1]];
  25
  26 minimize_lp(x, [x<1, y>0]);
  27 "Problem not bounded!";
  28
  29 minimize_lp(x, [x<1, y>0], all);
  30 [0,[y=0,x=0]];
  31
  32 minimize_lp(x, [x<1, y<-1], all);
  33 "Problem not feasible!";
  34
  35 block([A,b,c],
  36   A : read_matrix(file_search("Tests/afiro_A.csv"), 'csv),
  37   b : read_list(file_search("Tests/afiro_b.csv"), 'csv),
  38   c : read_list(file_search("Tests/afiro_c.csv"), 'csv),
  39   is (abs(second(linear_program(A, b, c) + 464.7531428571429))<10^-4)
  40 );
  41 true;
  42
  43 block([A,b,c],
  44   A : read_matrix(file_search("Tests/sc50a_A.csv"), 'csv),
  45   b : read_list(file_search("Tests/sc50a_b.csv"), 'csv),
  46   c : read_list(file_search("Tests/sc50a_c.csv"), 'csv),
  47   is (abs(second(linear_program(A, b, c) + 64.5750770585645))<10^-4)
  48 );
  49 true;
  50
  51 minimize_lp(-x, [1e-9*x + y <= 1], [x,y])$
  52 "Problem not bounded!"$
  53
  54 is(abs(first(minimize_lp(-x, [1e-9*x + y <= 1], [x,y])) + 1e9) < 1e-5), epsilon_lp=0 $
  55 true;
  56
  57 minimize_lp(-x, [10^-9*x + y <= 1], [x,y])$
  58 [- 1000000000, [y = 0, x = 1000000000]]$
  59
  60 /* Thanks to Michael Soegtrop
  61 * https://sourceforge.net/p/maxima/mailman/message/36414030/
  62 *
  63 * Do max-norm interpolation of a polynomial (x^3) on an interval
  64 * ([3/4,1]) with a quintic and sextic at 10 points. Using rational
  65 * arithmetic, we should get x^3; using floating point arithmetic and the
  66 * default value of epsilon_lp, we do for the quintic (1.0*x^3) but not
  67 * the sextic.
  68 *
  69 */
  70
  71 block([solutionpoly_lmad,f,sol5,sol6,sol5f,sol6f,checkerr],
  72   local(solutionpoly_lmad,f,checkerr),
  73   solutionpoly_lmad(f,degree,npoints,xmin,xmax,float_coerce):=block([xp,gep,gem,ineq,sol,yp,ydiff],
  74     xp:float_coerce(makelist(xmin+(xmax-xmin)*(i-1)/(npoints-1),i,1,npoints)),
  75     yp:makelist(apply(f,[xp[i]]),i,1,npoints),
  76     ydiff:makelist(yp[i]-sum(a[k]*xp[i]^k,k,0,degree),i,1,npoints),
  77     gep:makelist(umax>=ydiff[i],i,1,npoints),
  78     gem:makelist(umax>=-ydiff[i],i,1,npoints),
  79     ineq:append(gep,gem),
  80     sol:second(minimize_lp(umax,ineq)),
  81     /*[at(sum(a[i]*'x^i,i,0,degree),sol),at(umax,sol),sol]*/
  82     at(umax,sol)
  83     ),
  84   f(x):=x^3,
  85   sol5:solutionpoly_lmad(f,5,10,3/4,1,identity),
  86   sol6:solutionpoly_lmad(f,6,10,3/4,1,identity),
  87   sol5f:solutionpoly_lmad(f,5,10,3/4,1,'float),
  88   sol6f:solutionpoly_lmad(f,6,10,3/4,1,'float),
  89   checkerr(u,v) := is(abs(u)<=v),
  90   map(checkerr, [sol5,sol6,sol5f,sol6f,sol6f,sol6f],[0,0,1e-16,1e-7,1e-10,sol5f])
  91   );
  92 [true,true,true,true,false,false] $
  93
  94 (assume(c>4/3), declare(c,constant),
  95   ratsimp(maximize_lp(x+y, [y<=-x/c+3, y<=-x+4], [x, y]) - [4,[x = c/(c-1),y = (3*c-4)/(c-1)]])), epsilon_lp=0 $
  96 [0,[0 = 0,0 = 0]] $
  97
  98 (forget(c>4/3), remove(c,constant), kill(nonnegative_lp,return_input_lp), 0)$
  99 0 $