Repository: POptUS/IBCDFO (language: python)

Issue:
Adding gams files



Discussion / hints:

Fix patch:
diff --git a/goombah/m/goombah.m b/goombah/m/goombah.m
index ac1954a2..fb060c68 100644
--- a/goombah/m/goombah.m
+++ b/goombah/m/goombah.m
@@ -28,6 +28,7 @@
     % Deduce p from evaluating Ffun at x0
     try
         F0 = Ffun(x0);
+        F0 = F0(:)';
     catch
         warning('Problem using Ffun. Exiting');
         X = [];
diff --git a/goombah/m/subproblems/minimize_sum_quad_models_squared.gms b/goombah/m/subproblems/minimize_sum_quad_models_squared.gms
new file mode 100644
index 00000000..7f1c67a1
--- /dev/null
+++ b/goombah/m/subproblems/minimize_sum_quad_models_squared.gms
@@ -0,0 +1,150 @@
+* Solve the problem:
+*
+*       minimize sum M(x)**2
+*         s.t.   Low <= (x - x0) <= Upp
+*
+* where M is a defined by P quadratics mapping from R^n (0.5*(x - x0)^T * H_p * (x - x0) + g_p * (x - x0) + c_p)
+
+SETS
+  N   'elements of x'
+  P   'set of quadratic models'
+;
+
+Alias(N, M)
+
+PARAMETERS
+H(N, M, P)  'Model quadratic terms'
+g(N, P)     'Model linear terms'
+b(P)        'Model constant terms'
+solver      'Flag for solver'
+x0(N)       'Trust region center'
+x1(N)       'Candidate starting point for TRSP'
+Low(N)      'Lower bound on step'
+Upp(N)      'Upper bound on step'
+dsc(N)      'automatic scaling factor for d := x - x0'
+gamma(N)    'manual scaling factor for d := x - x0'
+;
+
+$if NOT exist quad_model_data.gdx $abort File quad_model_data.gdx does not exist: run the calling script from Matlab to create it.
+$gdxin quad_model_data.gdx
+$load N P H g b x0 x1 solver Low Upp
+$gdxin
+
+VARIABLES
+  tau     Objective value
+  x(N)    decision
+  d(N)    'd := (x - x0) / dsc, i.e. d * dsc := (x - x0)'
+  m_F(P)  value of each quadratic
+
+option decimals=8;
+
+* Declare model equations
+EQUATIONS
+obj                     Objective
+each_model
+* bounds_LB
+* bounds_UB
+;
+
+* Define model equations
+obj..               tau =e= sum(P, power(m_F(P), 2));
+
+* each_model(P)..     m_F(P) =e= 0.5*sum((N, M), (x(N) - x0(N))*H(N, M, P)*(x(M) - x0(M))) + sum(N, g(N, P)*(x(N) - x0(N))) + b(P);
+each_model(P)..     m_F(P) =e= 0.5*sum{(N, M), (d(N)*gamma(N))*H(N, M, P)*(d(M)*gamma(N))} + sum{N, g(N, P)*d(N)*gamma(N)} + b(P);
+
+* bounds_LB(N)..      x(N) - x0(N) =g= Low(N);
+* bounds_UB(N)..      x(N) - x0(N) =l= Upp(N);
+
+model TRSP / ALL /;
+
+dsc(N) = 1;
+gamma(N) = max[abs(low(N)),abs(upp(N))];
+
+$ifthen set AUTOSCALE
+  dsc(N) = gamma(N);
+  gamma(N) = 1;
+  TRSP.scaleopt = 1;
+$endif
+
+d.lo(N) = Low(N) / gamma(N);
+d.up(N) = Upp(N) / gamma(N);
+
+x.l(N) = x0(N);
+d.l(N) = 0;
+d.scale(N) = dsc(N);
+m_F.l(P) = b(P);
+tau.l = smin(P, power(m_F.l(P), 2));
+
+option limrow = 1000;
+option limcol = 1000;
+option Optca = 0;
+option Optcr = 0;
+option reslim = 30;
+
+* display x0;
+* display scale;
+* display Q;
+* display z;
+* display b;
+* display H;
+* display g;
+* display c;
+* display 'lastly, ', tau.l;
+* display 'lastly, ', m_F.l;
+
+TRSP.optfile = 1;
+TRSP.trylinear = 1;
+
+$onecho > minos.opt
+Feasibility tolerance  0
+Optimality  tolerance  0
+$offecho
+
+* option DNLP=LINDOGLOBAL;
+* option DNLP=CONOPT;
+* option DNLP=couenne;
+* option DNLP=mosek;
+* option DNLP=minos;
+* option DNLP=examiner;
+* option DNLP=snopt;
+
+if (solver = 1,
+  option DNLP=conopt;
+);
+if (solver = 2,
+$onecho > minos.opt
+$offecho
+  option DNLP=minos;
+);
+if (solver = 3,
+$onecho > snopt.opt
+$offecho
+  option DNLP=snopt;
+);
+if (solver = 4,
+$onecho > baron.opt
+threads 100
+$offecho
+  option DNLP=baron;
+);
+
+SOLVE TRSP MINIMIZING tau USING DNLP;
+
+scalars modelStat, solveStat;
+modelStat = TRSP.modelstat;
+solveStat = TRSP.solvestat;
+
+parameters
+  Nerr(N)
+  xd(N)    'x - x0'
+  ;
+xd(N) = d.L(N)*gamma(N);
+Nerr(N) = max(0, Low(N) - xd(N));
+* abort$[sum{N, Nerr(N)}] 'x-x0 < Low', Nerr;
+Nerr(N) = max(0, xd(N)-Upp(N));
+* abort$[sum{N, Nerr(N)}] 'x-x0 > Upp', Nerr;
+
+m_F.L(P) = 0.5*sum{(N, M), xd(N)*H(N, M, P)*xd(M)} + sum{N, g(N, P)*xd(N)} + b(P);
+
+x.L(N) = x0(N) + xd(N);
+execute_unload 'solution', modelStat, solveStat, x, tau, m_F;
diff --git a/manifold_sampling/m/check_inputs_and_initialize.m b/manifold_sampling/m/check_inputs_and_initialize.m
index c5fcef25..a6dbd9ff 100644
--- a/manifold_sampling/m/check_inputs_and_initialize.m
+++ b/manifold_sampling/m/check_inputs_and_initialize.m
@@ -4,7 +4,7 @@
 
     n = length(x0);
     p = length(F0);
-    delta = 0.3;
+    delta = 0.1;
     printf = 1;
 
     h_activity_tol = min(1e-8, delta);
@@ -14,10 +14,11 @@
     fq_pars.Par(2) = max(10, sqrt(n)); % [dbl] delta multiplier for all interp. points
     fq_pars.Par(3) = 1e-3;  % [dbl] Pivot threshold for validity (1e-5)
     fq_pars.Par(4) = .001;  % [dbl] Pivot threshold for additional points (.001)
-    fq_pars.npmax = (n + 1) * (n + 2) / 2;     % [int] number of points in model building
+    % fq_pars.npmax = (n + 1) * (n + 2) / 2;     % [int] number of points in model building
+    fq_pars.npmax = 2 * n + 1;     % [int] number of points in model building
 
     % Internal parameters/tolerances for manifold sampling
-    tol.maxdelta = 1;
+    tol.maxdelta = 1e8;
     tol.mindelta = 1e-13;
     tol.gtol = 1e-13;
     tol.norm_g_change =  1e-3;  % Tolerance for the change of norm(g_k)
diff --git a/manifold_sampling/m/h_examples/sum_squared.m b/manifold_sampling/m/h_examples/sum_squared.m
new file mode 100644
index 00000000..077110f7
--- /dev/null
+++ b/manifold_sampling/m/h_examples/sum_squared.m
@@ -0,0 +1,32 @@
+function [h, grads, Hash] = sum_squared(z, H0)
+% Evaluates the sum-sqaured objective
+%   min sum_j { z_j^2 }
+%
+% Inputs:
+%  z:              [1 x p]   point where we are evaluating h
+%  H0: (optional)  [1 x l cell of strings]  set of hashes where to evaluate z
+
+% Outputs:
+%  h: [dbl]                       function value
+%  grads: [p x l]                 gradients of each of the l manifolds active at z
+%  Hash: [1 x l cell of strings]  set of hashes for each of the l manifolds active at z (in the same order as the elements of grads)
+
+z = z(:);
+n = length(z);
+
+if nargin == 1
+
+    h = sum(z.^2);
+
+    grads = zeros(n, 1);
+    Hash = cell(1, 1);
+    Hash{1} = int2str(1);
+    grads(:) = 2 * z;
+
+elseif nargin == 2
+    h = sum(z.^2);
+    grads = zeros(length(z), 1);
+    grads(:) = 2 * z;
+else
+    error('Too many inputs to function');
+end


Test patch:
diff --git a/goombah/m/tests/compare_goombah_and_pounders.m b/goombah/m/tests/compare_goombah_and_pounders.m
new file mode 100644
index 00000000..820d6776
--- /dev/null
+++ b/goombah/m/tests/compare_goombah_and_pounders.m
@@ -0,0 +1,98 @@
+% This wrapper tests various algorithms against the Benchmark functions from the
+% More and Wild SIOPT paper "Benchmarking derivative-free optimization algorithms"
+function [] = compare_goombah_and_pounders()
+addpath('../../../../BenDFO/data/');
+addpath('../../../../BenDFO/m/');
+addpath('../../../minq/m/minq5');
+addpath('../../../goombah/m');
+addpath('../../../pounders/m/general_h_funs');
+addpath('../../../pounders/m');
+addpath('../../../goombah/m/subproblems/');
+addpath('../../../manifold_sampling/m');
+addpath('../../../manifold_sampling/m/h_examples/');
+
+load dfo.dat;
+
+spsolver = 2; % TRSP Solver
+nfmax = 100;
+gtol = 1e-13;
+factor = 10;
+
+for row = 1:length(dfo)
+    nprob = dfo(row, 1);
+    n = dfo(row, 2);
+    m = dfo(row, 3);
+    factor_power = dfo(row, 4);
+
+    BenDFO.nprob = nprob;
+    BenDFO.m = m;
+    BenDFO.n = n;
+
+    objective = @(x)calfun_wrapper(x, BenDFO, 'smooth');
+
+    X0 = dfoxs(n, nprob, factor^factor_power)';
+    npmax = 2 * n + 1;  % Maximum number of interpolation points [2*n+1]
+    L = -Inf(1, n);
+    U = Inf(1, n);
+    nfs = 0;
+    F0 = [];
+    xkin = 1;
+    delta = 0.1;
+    printf = 0;
+
+    GAMS_options.file = '../../../goombah/m/subproblems/minimize_sum_quad_models_squared.gms';
+    GAMS_options.solvers = 1:4;
+    subprob_switch = 'linprog';
+
+    for hfun_cases = 1:1
+        Results = cell(2, 3, 53);
+        if hfun_cases == 1
+            hfun = @(F)sum(F.^2);
+            combinemodels = @leastsquares;
+        elseif hfun_cases == 2
+            error("Not implemented in GOOMBAH yet");
+            alpha = 0; % If changed here, also needs to be adjusted in squared_diff_from_mean.m
+            hfun = @(F)sum((F - 1 / length(F) * sum(F)).^2) - alpha * (1 / length(F) * sum(F))^2;
+            combinemodels = @squared_diff_from_mean;
+        elseif hfun_cases == 3
+            error("Not implemented in GOOMBAH yet");
+            if m ~= 3 % Emittance is only defined for the case when m == 3
+                continue
+            end
+            hfun = @emittance_h;
+            combinemodels = @emittance_combine;
+        end
+        disp([row, hfun_cases]);
+
+        filename = ['./benchmark_results/poundersM_and_GOOMBAH_nfmax=' int2str(nfmax) '_gtol=' num2str(gtol) '_prob=' int2str(row) '_spsolver=' num2str(spsolver) '_hfun=' func2str(combinemodels) '.mat'];
+
+        for method = [1, 2]
+            if method == 1
+                [X, F, flag, xk_best] = pounders(objective, X0, n, npmax, nfmax, gtol, delta, nfs, m, F0, xkin, L, U, printf, spsolver, hfun, combinemodels);
+                Results{method, hfun_cases, row}.alg = 'POUNDERs';
+            elseif method == 2
+                [X, F, h, xkin] = goombah(@sum_squared, objective, nfmax, X0, L, U, GAMS_options, subprob_switch);
+                Results{method, hfun_cases, row}.alg = 'GOOMBAH';
+            end
+
+            evals = size(F, 1);
+            h = zeros(evals, 1);
+            for i = 1:evals
+                h(i) = hfun(F(i, :));
+            end
+        end
+
+        Results{method, hfun_cases, row}.problem = ['problem ' num2str(row) ' from More/Wild'];
+        Results{method, hfun_cases, row}.Fvec = F;
+        Results{method, hfun_cases, row}.H = h;
+        Results{method, hfun_cases, row}.X = X;
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %     save('-mat7-binary', filename, 'Results') % Octave save
+        save(filename, 'Results');
+    end
+end
+end
+
+function [fvec] = calfun_wrapper(x, struct, probtype)
+[~, fvec, ~] = calfun(x, struct, probtype);
+end