source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/modules/global/cutsdp.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 13.1 KB
Line 
1function output = cutsdp(p)
2% CUTSDP
3%
4% See also SOLVESDP, BNB, BINVAR, INTVAR, BINARY, INTEGER, LMI
5
6% Author Johan Löfberg
7% $Id: cutsdp.m,v 1.5 2006/05/23 21:55:59 joloef Exp $
8
9% *************************************************************************
10%% INITIALIZE DIAGNOSTICS IN YALMIP
11% *************************************************************************
12bnbsolvertime = clock;
13showprogress('Cutting plane solver started',p.options.showprogress);
14
15% *************************************************************************
16%% If we want duals, we may not extract bounds. However, bounds must be
17% extracted in discrete problems.
18% *************************************************************************
19if p.options.cutsdp.recoverdual
20    warning('Dual recovery not implemented yet in CUTSDP')
21end
22
23% *************************************************************************
24%% Define infinite bounds
25% *************************************************************************
26if isempty(p.ub)
27    p.ub = repmat(inf,length(p.c),1);
28end
29if isempty(p.lb)
30    p.lb = repmat(-inf,length(p.c),1);
31end
32
33% *************************************************************************
34%% ADD CONSTRAINTS 0<x<1 FOR BINARY
35% *************************************************************************
36if ~isempty(p.binary_variables)
37    p.ub(p.binary_variables) =  min(p.ub(p.binary_variables),1);
38    p.lb(p.binary_variables) =  max(p.lb(p.binary_variables),0);
39end
40
41% *************************************************************************
42%% Extract better bounds from model
43% *************************************************************************
44if ~isempty(p.F_struc)
45    [lb,ub,used_rows] = findulb(p.F_struc,p.K);
46    if ~isempty(used_rows)
47        lower_defined = find(~isinf(lb));
48        if ~isempty(lower_defined)
49            p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
50        end
51        upper_defined = find(~isinf(ub));
52        if ~isempty(upper_defined)
53            p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
54        end
55        p.F_struc(p.K.f+used_rows,:)=[];
56        p.K.l = p.K.l - length(used_rows);
57    end
58end
59
60% *************************************************************************
61%% ADD CONSTRAINTS 0<x<1 FOR BINARY
62% *************************************************************************
63if ~isempty(p.binary_variables)
64    p.ub(p.binary_variables) =  min(p.ub(p.binary_variables),1);
65    p.lb(p.binary_variables) =  max(p.lb(p.binary_variables),0);
66end
67
68p.ub = min(p.ub,p.options.cutsdp.variablebound');
69p.lb = max(p.lb,-p.options.cutsdp.variablebound');
70
71% *************************************************************************
72%% PRE-SOLVE (nothing fancy coded)
73% *************************************************************************
74if isempty(find(isinf([p.ub;p.lb]))) & p.K.l>0
75    [p.lb,p.ub] = tightenbounds(-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end),p.F_struc(1+p.K.f:p.K.f+p.K.l,1),p.lb,p.ub,p.integer_variables);
76end
77
78% *************************************************************************
79%% PERTURBATION OF LINEAR COST
80% *************************************************************************
81p.corig = p.c;
82if nnz(p.Q)~=0
83    g = randn('seed');
84    randn('state',1253); %For my testing, I keep this the same...
85    % This perturbation has to be better. Crucial for many real LP problems
86    p.c = (p.c).*(1+randn(length(p.c),1)*1e-4);
87    randn('seed',g);
88end
89
90% *************************************************************************
91%% We don't need this
92% *************************************************************************
93p.options.savesolverinput  = 0;
94p.options.savesolveroutput = 0;
95
96% *************************************************************************
97%% Display logics
98% 0 : Silent
99% 1 : Display cut progress
100% 2 : Display node solver prints
101% *************************************************************************
102switch max(min(p.options.verbose,3),0)
103    case 0
104        p.options.cutsdp.verbose = 0;
105    case 1
106        p.options.cutsdp.verbose = 1;
107        p.options.verbose = 0;
108    case 2
109        p.options.cutsdp.verbose = 2;
110        p.options.verbose = 0;
111    case 3
112        p.options.cutsdp.verbose = 2;
113        p.options.verbose = 1;
114    otherwise
115        p.options.cutsdp.verbose = 0;
116        p.options.verbose = 0;
117end
118
119% *************************************************************************
120%% START CUTTING
121% *************************************************************************
122[x_min,solved_nodes,lower,feasible,D_struc] = cutting(p);
123%% --
124
125% *************************************************************************
126%% CREATE SOLUTION
127% *************************************************************************
128output.problem = 0;
129if ~feasible
130    output.problem = 1;
131end
132if solved_nodes == p.options.cutsdp.maxiter
133    output.problem = 3;
134end
135output.solved_nodes = solved_nodes;
136output.Primal       = x_min;
137output.Dual = D_struc;
138output.Slack = [];
139output.solverinput  = 0;
140output.solveroutput =[];
141output.solvertime   = etime(clock,bnbsolvertime);
142%% --
143
144function [x,solved_nodes,lower,feasible,D_struc] = cutting(p)
145
146% *************************************************************************
147%% Sanity check
148% *************************************************************************
149if any(p.lb>p.ub)
150    x = zeros(length(p.c),1);
151    solved_nodes = 0;
152    lower = inf;
153    feasible = 0;
154    D_struc = [];
155    return
156end
157
158% *************************************************************************
159%% Create function handle to solver
160% *************************************************************************
161cutsolver = p.solver.cutsolver.call;
162
163% *************************************************************************
164%% Create copy of model without
165%  the SDP part
166% *************************************************************************
167p_lp = p;
168p_lp.F_struc = p_lp.F_struc(1:p.K.l+p.K.f,:);
169p_lp.K.s = 0;
170
171% *************************************************************************
172%% DISPLAY HEADER
173% *************************************************************************
174if p.options.cutsdp.verbose
175    disp('* Starting YALMIP cutting plane for MISDP based on MILP');
176    disp(['* Lower solver   : ' p.solver.cutsolver.tag]);
177    disp(['* Max iterations : ' num2str(p.options.cutsdp.maxiter)]);
178end
179
180if p.options.bnb.verbose;            disp(' Node       Infeasibility.     Lower    LP cuts');end;
181
182%% Initialize diagnostic
183infeasibility = -inf;
184solved_nodes = 0;
185feasible = 1;
186lower = -inf;
187saveduals = 1;
188
189% *************************************************************************
190%% Add diagonal cuts to begin with
191% *************************************************************************
192savedCuts = [];
193savedIndicies = [];
194if p.K.s(1)>0
195    top = p.K.f+p.K.l+1;
196    for i = 1:length(p.K.s)
197        n = p.K.s(i);
198        newF=[];
199        for m = 1:p.K.s(i)
200            d = eyev(p.K.s(i),m);
201            index = (1+(m-1)*(p.K.s(i)+1));
202            newF = [newF;p.F_struc(top+index-1,:)];
203        end
204        % Clean
205        newF(abs(newF)<1e-12) = 0;
206        keep=find(any(newF(:,2:end),2));
207        newF = newF(keep,:);
208
209        p_lp.F_struc = [p_lp.F_struc;newF];
210        p_lp.K.l = p_lp.K.l + size(newF,1);
211        top = top+n^2;
212    end
213end
214
215goon = 1;
216rr = p_lp.integer_variables;
217rb = p_lp.binary_variables;
218only_solvelp = 0;
219pplot = 0;
220
221% *************************************************************************
222% Crude lower bound
223% FIX for quadratic case
224% *************************************************************************
225lower = 0;
226if nnz(p.Q) == 0
227    for i = 1:length(p.c)
228        if p.c(i)>0
229            if isinf(p.lb(i))
230                lower = -inf;
231                break
232            else
233                lower = lower + p.c(i)*p.lb(i);
234            end
235        elseif p.c(i)<0
236            if isinf(p.ub(i))
237                lower = -inf;
238                break
239            else
240                lower = lower + p.c(i)*p.ub(i);
241            end
242        end
243    end
244end
245%lower = sum(sign(p.c).*(p.lb));
246if isinf(lower) | nnz(p.Q)~=0
247    lower = -1e6;
248end
249
250% *************************************************************************
251% Experimental stuff for variable fixing
252% *************************************************************************
253if p.options.cutsdp.nodefix & (p.K.s(1)>0)
254    top=1+p.K.f+p.K.l+sum(p.K.q);
255    for i=1:length(p.K.s)
256        n=p.K.s(i);
257        for j=1:size(p.F_struc,2)-1;
258            X=full(reshape(p.F_struc(top:top+n^2-1,j+1),p.K.s(i),p.K.s(i)));
259            X=(X+X')/2;
260            v=real(eig(X+sqrt(eps)*eye(length(X))));
261            if all(v>=0)
262                sdpmonotinicity(i,j)=-1;
263            elseif all(v<=0)
264                sdpmonotinicity(i,j)=1;
265            else
266                sdpmonotinicity(i,j)=nan;
267            end
268        end
269        top=top+n^2;
270    end
271else
272    sdpmonotinicity=[];
273end
274
275while goon
276
277    % Add lower bound
278    if ~isinf(lower)
279        p_lp.F_struc = [p_lp.F_struc;-lower p_lp.c'];
280        p_lp.K.l = p_lp.K.l + 1;
281    end
282
283    if p.options.cutsdp.nodetight
284        % Extract LP part Ax<=b
285        A = -p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),2:end);
286        b = p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),1);
287        c = p_lp.c;
288        % Tighten bounds and find redundant constraints
289        [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
290        A(redundant,:) = [];
291        b(redundant) = [];
292        p_lp.F_struc(p_lp.K.f+redundant,:) = [];
293        p_lp.K.l = p_lp.K.l-length(redundant);
294    end
295
296    if p.options.cutsdp.nodefix
297        % Try to find variables to fix w.l.o.g
298        [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity);
299        p_lp.lb(fix_up)   = p_lp.ub(fix_up);
300        p_lp.ub(fix_down) = p_lp.lb(fix_down);
301        while ~(isempty(fix_up) & isempty(fix_down))
302            [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
303            A(redundant,:) = [];
304            b(redundant) = [];
305            p_lp.F_struc(p_lp.K.f+redundant,:) = [];
306            p_lp.K.l = p_lp.K.l-length(redundant);
307            fix_up = [];
308            fix_down = [];
309            % Try to find variables to fix w.l.o.g
310            [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity);
311            p_lp.lb(fix_up)   = p_lp.ub(fix_up);
312            p_lp.ub(fix_down) = p_lp.lb(fix_down);
313        end
314    end
315
316
317    output = feval(cutsolver,p_lp);
318
319    % Remove lower bound (avoid accumulating them)
320    if ~isinf(lower)
321        p_lp.K.l = p_lp.K.l - 1;
322        p_lp.F_struc = p_lp.F_struc(1:end-1,:);
323    end
324
325    if output.problem == 1 | output.problem == 12
326        % LP relaxation was infeasible, hence problem is infeasible
327        feasible = 0;
328        lower = inf;
329        goon = 0;
330        x = zeros(length(p.c),1);
331        lower = inf;
332    else
333        % Relaxed solution
334        x = output.Primal;
335        lower = p.f+p.c'*x+x'*p.Q*x;
336
337        infeasibility = 0;
338        if p.K.s(1)>0
339            % Add cuts
340            top = p.K.f+p.K.l+1;
341            for i = 1:1:length(p.K.s)
342                n = p.K.s(i);
343                X = p.F_struc(top:top+n^2-1,:)*[1;x];
344                X = reshape(X,n,n);
345                Y = randn(n,n);
346                newcuts = 1;
347                newF = zeros(n,size(p.F_struc,2));
348                [d,v] = eig(X);
349                infeasibility = min(infeasibility,min(diag(v)));
350                dummy=[];
351                newF = [];
352                if infeasibility<0
353                    [ii,jj] = sort(diag(v));
354                    for m = jj(1:min(length(jj),p.options.cutsdp.cutlimit))'%find(diag(v<0))%1:1%length(v)
355                        if v(m,m)<0
356                            bA =  d(:,m)'*(kron(d(:,m),speye(n))'*p.F_struc(top:top+n^2-1,:));
357                            b = bA(:,1);
358                            A = -bA(:,2:end);
359                            newF = [newF;b -A];
360                            %f = b-floor(b);
361                            %fj = A - floor(A);
362                            %A = floor(A) + max((fj-f),0)./(1-f);
363                            %b = floor(b);
364                            %newF = [newF;b -A];
365                            newcuts = newcuts + 1;
366                        end
367                    end
368                end
369                newF(abs(newF)<1e-12) = 0;
370                keep=find(any(newF(:,2:end),2));
371                newF = newF(keep,:);
372                if size(newF,1)>0
373                    p_lp.F_struc = [p_lp.F_struc;newF];
374                    p_lp.K.l = p_lp.K.l + size(newF,1);
375                    [i,j] = sort(p_lp.F_struc*[1;x]);
376                end
377                top = top+n^2;
378            end
379        end
380      %    res = find(p_lp.F_struc*[1;x] > mean(abs(p_lp.F_struc*[1;x])));
381      %    p_lp.F_struc(res,:) = [];
382      %    p_lp.K.l = p_lp.K.l-length(res);
383        goon = infeasibility <= p.options.cutsdp.feastol;
384        goon = goon & feasible;
385        goon = goon & (solved_nodes < p.options.cutsdp.maxiter-1);
386    end
387
388    solved_nodes = solved_nodes + 1;
389    if p.options.cutsdp.verbose;fprintf(' %4.0f :      %12.3E      %12.3E      %2.0f\n',solved_nodes,infeasibility,lower,p_lp.K.l-p.K.l);end
390end
391
392D_struc = [];
Note: See TracBrowser for help on using the repository browser.