source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/@sdpvar/mtimes.m @ 37

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

Added original make3d
File size: 18.7 KB
Line 
1function Z = mtimes(X,Y)
2%MTIMES (overloaded)
3
4% Author Johan Löfberg
5% $Id: mtimes.m,v 1.57 2006/10/19 12:42:36 joloef Exp $
6
7% Check classes
8X_is_spdvar = isa(X,'sdpvar');
9Y_is_spdvar = isa(Y,'sdpvar');
10
11% Convert block objects
12if ~X_is_spdvar
13    if isa(X,'blkvar')
14        X = sdpvar(X);
15        X_is_spdvar = isa(X,'sdpvar');
16    end
17end
18
19if ~Y_is_spdvar
20    if isa(Y,'blkvar')
21        Y = sdpvar(Y);
22        Y_is_spdvar = isa(Y,'sdpvar');
23    end
24end
25
26% Lame special cases, make sure to reurn
27% empty matrices in the sense that the
28% used MATLAB version
29if isempty(X)
30    YY = full(reshape(Y.basis(:,1),Y.dim(1),Y.dim(2)));
31    Z = X*YY;
32    return
33elseif isempty(Y)
34    XX = full(reshape(X.basis(:,1),X.dim(1),X.dim(2)));
35    Z = XX*Y;
36    return
37end
38
39
40% Optimized calls in different order?
41if X_is_spdvar & Y_is_spdvar
42    manytimesfew = length(X.lmi_variables) > 5*length(Y.lmi_variables);
43    if manytimesfew
44        Z = (Y'*X')'; % Optimized for this order (few variables * many variables)
45        return
46    end
47end
48
49% Different code for
50% 1 : SDPVAR * DOUBLE
51% 2 : DOUBLE * SDPVAR
52% 3 : SDPVAR * SDPVAR
53switch 2*X_is_spdvar+Y_is_spdvar
54    case 3
55        try
56            x_isscalar =  (X.dim(1)*X.dim(2)==1);
57            y_isscalar =  (Y.dim(1)*Y.dim(2)==1);
58
59            % Optimized unique
60            all_lmi_variables = uniquestripped([X.lmi_variables Y.lmi_variables]);
61
62            % Create clean SDPVAR object
63            Z = X;Z.dim(1) = 1;Z.dim(2) = 1;Z.lmi_variables = all_lmi_variables;Z.basis = [];
64
65            % Awkward code due to bug in ML6.5
66            Xbase = reshape(X.basis(:,1),X.dim(1),X.dim(2));
67            Ybase = reshape(Y.basis(:,1),Y.dim(1),Y.dim(2));
68            if x_isscalar
69                Xbase = sparse(full(Xbase));
70            end
71            if y_isscalar
72                Ybase = sparse(full(Ybase));
73            end
74
75            index_X = double(ismembc(all_lmi_variables,X.lmi_variables));
76            index_Y = double(ismembc(all_lmi_variables,Y.lmi_variables));
77            iX=find(index_X);
78            iY=find(index_Y);
79            index_X(iX)=1:length(iX);index_X=index_X(:);
80            index_Y(iY)=1:length(iY);index_Y=index_Y(:);
81
82            ny = Y.dim(1);
83            my = Y.dim(2);
84            nx = X.dim(1);
85            mx = X.dim(2);
86
87            % Pre-allocate sufficiently long
88            Z.lmi_variables = [Z.lmi_variables zeros(1,length(X.lmi_variables)*length(Y.lmi_variables))];
89
90            % Pre-calc identity (used a lot
91            speyemy = sparse(1:my,1:my,1,my,my);
92
93            % Linear terms
94            inner_vector_product = (X.dim(1)==1 & Y.dim(2)==1 & (X.dim(2) == Y.dim(1)));
95            if inner_vector_product
96                base1=Xbase*Y.basis;base1=base1(2:end);
97                base2=Ybase.'*X.basis;base2=base2(2:end);
98                [i1,j1,k1]=find(base1);
99                [i2,j2,k2]=find(base2);
100                base1 = sparse(i1,iY(j1),k1,1,length(all_lmi_variables));
101                base2 = sparse(i2,iX(j2),k2,1,length(all_lmi_variables));
102                Z.basis = [Xbase*Ybase base1+base2];
103            else
104                base0 = Xbase*Ybase;
105                if x_isscalar
106                    base1 = Xbase*Y.basis(:,2:end);
107                    base2 = Ybase(:)*X.basis(:,2:end);
108                elseif y_isscalar
109                    base1 = Xbase(:)*Y.basis;base1=base1(:,2:end);
110                    base2 = X.basis*Ybase;base2=base2(:,2:end);
111                else
112                    base1 = kron(speyemy,Xbase)*Y.basis(:,2:end);
113                    base2 = kron(Ybase.',speye(nx))*X.basis(:,2:end);
114                end
115                [i1,j1,k1]=find(base1);
116                [i2,j2,k2]=find(base2);
117                base1 = sparse(i1,iY(j1),k1,size(base0(:),1),length(all_lmi_variables));
118                base2 = sparse(i2,iX(j2),k2,size(base0(:),1),length(all_lmi_variables));
119                Z.basis = [base0(:) base1+base2];
120            end
121
122            % Loop start for nonlinear terms
123            i = length(all_lmi_variables)+1;
124
125            [mt,oldvariabletype,mt_hash,hash] = yalmip('monomtable');
126
127            % Check if problem is bilinear. We can exploit this later
128            % to improve performance significantly...
129            bilinearproduct = 0;
130            candofastlocation  = 0;
131            if all(oldvariabletype(X.lmi_variables)==0) & all(oldvariabletype(Y.lmi_variables)==0)
132               % if isempty(intersect(X.lmi_variables,Y.lmi_variables))
133                if ~any(ismembc(X.lmi_variables,Y.lmi_variables))
134                    bilinearproduct = 1;
135                    try
136                        dummy = ismembc2(1,1); % Not available in all versions (needed in ismember)
137                        candofastlocation = 1;
138                    catch
139                    end
140                end
141            end
142
143            oldmt = mt;
144            local_mt = mt(all_lmi_variables,:);
145            used_variables = any(local_mt,1);
146            local_mt = local_mt(:,used_variables);
147
148            possibleOld = find(any(mt(:,used_variables),2));
149            if all(oldvariabletype <=3)
150                % All monomials have non-negative integer powers
151                % no chance of x^2*x^-1, hence all products
152                % are nonlinear
153                possibleOld = possibleOld(find(oldvariabletype(possibleOld)));
154                if size(possibleOld,1)==0
155                    possibleOld = [];
156                end
157            end
158
159            if bilinearproduct & ~isempty(possibleOld)
160                if length(X.lmi_variables)<=length(Y.lmi_variables)
161                    temp = mt(:,X.lmi_variables);
162                    temp = temp(possibleOld,:);
163                    possibleOld=possibleOld(find(any(temp,2)));
164                else
165                    temp = mt(:,Y.lmi_variables);
166                    temp = temp(possibleOld,:);
167                    possibleOld=possibleOld(find(any(temp,2)));
168                end
169            end
170
171            theyvars = find(index_Y);
172            thexvars = find(index_X);
173
174            possibleOldHash = mt_hash(possibleOld);
175            oldhash = hash;
176            hash = hash(used_variables);
177            new_mt_hash = [];
178            new_mt_hash_transpose = [];
179            new_mt = [];
180            changed_mt = 0;
181            local_mt = local_mt';
182            nvar = size(mt,1);
183
184            for ix = thexvars(:)'
185
186                mt_x = local_mt(:,ix);
187
188                testthese = theyvars(:)';
189
190                % Compute [vec(Xbasis*Ybasis1) vec(Xbasis*Ybasis2) ...]
191                % in one shot using vectorization and Kronecker tricks
192                % Optimized and treat special case scalar*matrix etc
193                if x_isscalar
194                    Xibase = X.basis(:,1+index_X(ix));
195                    allprodbase = Xibase * Y.basis(:,1+index_Y(testthese));
196                elseif y_isscalar
197                    Xibase = X.basis(:,1+index_X(ix));
198                    allprodbase =  Xibase * Y.basis(:,1+index_Y(testthese));
199                elseif inner_vector_product
200                    Xibase = X.basis(:,1+index_X(ix)).';
201                    allprodbase = Xibase*Y.basis(:,1+index_Y(testthese));
202                else
203                    Xibase = reshape(X.basis(:,1+index_X(ix)),nx,mx);                 
204                    temp = kron(speyemy,Xibase);
205                    allprodbase = temp * Y.basis(:,1+index_Y(testthese));
206                end
207
208                % Keep non-zero matrices
209                nonzeromatrices = find(sum(abs(allprodbase),1)>1e-12);
210                testthese = testthese(nonzeromatrices);
211                allprodbase = allprodbase(:,nonzeromatrices);
212
213                % Some data for vectorization
214                nyvars = length(testthese);
215                if prod(size(mt_x))==1 % Bug in Solaris and Linux, ML 6.X
216                    allmt_xplusy = local_mt(:,testthese) + sparse(repmat(full(mt_x),1,nyvars));
217                else
218                    allmt_xplusy = local_mt(:,testthese) + repmat(mt_x,1,nyvars);
219                end
220                allhash = allmt_xplusy'*hash;
221                thesewhereactuallyused = zeros(1,nyvars);
222                copytofrom  = ones(1,nyvars);
223                acounter =  0;
224
225
226                % Special case : x*inv(x) and similiar...
227                sum_to_constant = abs(allhash)<eps;
228                add_these = find(sum_to_constant);
229                if ~isempty(add_these)
230                    prodbase = allprodbase(:,add_these);
231                    Z.basis(:,1) = Z.basis(:,1) + sum(prodbase,2);
232                    copytofrom(add_these) = 0;
233                end
234                indicies = find(~sum_to_constant);
235                indicies = indicies(:)';
236
237                allbefore_in_old = 1;
238                if bilinearproduct & candofastlocation
239                    [dummy,allbefore_in_old] = ismember(allhash,possibleOldHash);
240                end
241
242
243                if bilinearproduct & candofastlocation & (nnz(allbefore_in_old)==0)
244                    % All nonlinear variables are new, so we can create them at once
245                    changed_mt=1;
246                    thesewhereactuallyused = thesewhereactuallyused+1;
247                    Z.lmi_variables(i:(i+length(indicies)-1)) = (nvar+1):(nvar+length(indicies));
248                    nvar = nvar + length(indicies);
249                    i = i + length(indicies);
250                else
251                    isemptynew_mt_hash = isempty(new_mt_hash);                   
252                   
253                    for acounter = indicies
254
255                        current_hash = allhash(acounter);
256
257                        % Ok, braze your self for some horrible special case
258                        % treatment etc...
259                        if isemptynew_mt_hash | bilinearproduct % only search among old monomials
260                            if bilinearproduct & candofastlocation
261                                before = allbefore_in_old(acounter);
262                                if before==0
263                                    before = [];
264                                else
265                                    before = possibleOld(before);
266                                end
267                            else
268                                before = possibleOld(findhash(possibleOldHash,current_hash));                               
269%                                before = possibleOld(findhash(possibleOldHash,current_hash,length(possibleOldHash)));                               
270                            end
271                        else
272                            before = findhash(new_mt_hash,current_hash); % first among new monomials
273%                            before = findhash(new_mt_hash_transpose,current_hash,topp); % first among new monomials
274                            if before % Try new if we failed finding any match
275                                before=before+size(mt_hash,1);
276                            else
277                                before = possibleOld(findhash(possibleOldHash,current_hash));                               
278%                                before = possibleOld(findhash(possibleOldHash,current_hash,length(possibleOldHash)));                               
279                            end
280                        end
281                        if before
282                            Z.lmi_variables(i) = before;
283                        else
284                            changed_mt=1;
285                            isemptynew_mt_hash=0;
286                            thesewhereactuallyused(acounter) = 1;
287                            if ~bilinearproduct % We don't need to update until later since it not is used inside the loop
288                                new_mt_hash = [new_mt_hash;current_hash];
289 %                               new_mt_hash_transpose(topp) = current_hash;topp = topp + 1;
290                            end
291                            nvar = nvar + 1;
292                            Z.lmi_variables(i) = nvar;                           
293                        end
294                        i = i+1;
295                    end
296                end
297                if all(copytofrom)
298                    Z.basis = [Z.basis allprodbase];
299                else
300                    Z.basis = [Z.basis allprodbase(:,find(copytofrom))];
301                end
302                if all(thesewhereactuallyused)
303                    new_mt = [new_mt allmt_xplusy];
304                else
305                    new_mt = [new_mt allmt_xplusy(:,find(thesewhereactuallyused))];
306                end
307                if bilinearproduct
308                    new_mt_hash = [new_mt_hash;allhash(find(thesewhereactuallyused))];
309                end
310            end
311%             
312%             try
313%            new_mt_hash_transpose = new_mt_hash_transpose(1:topp-1);
314%            norm(new_mt_hash_transpose-new_mt_hash);
315%             catch
316%                 1
317%             end
318           
319            if ~isempty(new_mt)
320                [i1,j1,k1] = find(mt);
321                [ii1,jj1,kk1] = find(new_mt');
322                uv = find(used_variables);uv=uv(:);
323                mt = sparse([i1(:);ii1(:)+size(mt,1)],[j1(:);uv(jj1(:))],[k1(:);kk1(:)],size(mt,1)+size(new_mt,2),size(mt,2));
324            end
325
326            % We pre-allocated a sufficiently long, now pick the ones we
327            % actually filled with values
328            Z.lmi_variables = Z.lmi_variables(1:i-1);
329
330            % Fucked up order (lmi_variables should be sorted)
331            if any(diff(Z.lmi_variables)<0)
332                [i,j]=sort(Z.lmi_variables);
333                Z.basis = [Z.basis(:,1) Z.basis(:,j+1)];
334                Z.lmi_variables = Z.lmi_variables(j);
335            end
336
337            [un_Z_vars2] = uniquestripped(Z.lmi_variables);
338            if length(un_Z_vars2) < length(Z.lmi_variables)
339                [un_Z_vars,hh,jj] = unique(Z.lmi_variables);
340                if length(Z.lmi_variables) ~=length(un_Z_vars)
341                    Z.basis = Z.basis*sparse([1 1+jj],[1 1+(1:length(jj))],ones(1,1+length(jj)))';
342                    Z.lmi_variables = un_Z_vars;
343                end
344            end
345
346            if changed_mt%~isequal(mt,oldmt)
347                newmt = mt(size(oldmt,1)+1:end,:);
348                nonlinear = ~(sum(newmt,2)==1 & sum(newmt~=0,2)==1);
349                newvariabletype = spalloc(size(newmt,1),1,nnz(nonlinear))';
350                nonlinearvariables = find(nonlinear);
351                newvariabletype = sparse(nonlinearvariables,ones(length(nonlinearvariables),1),3,size(newmt,1),1)';
352                if ~isempty(nonlinear)
353                    %mt = internal_sdpvarstate.monomtable;
354                    %newvariabletype(nonlinear) = 3;
355                    quadratic = sum(newmt,2)==2;
356                    newvariabletype(quadratic) = 2;
357                    bilinear = max(newmt,[],2)<=1;
358                    newvariabletype(bilinear & quadratic) = 1;
359                    sigmonial = any(0>newmt,2) | any(newmt-fix(newmt),2);
360                    newvariabletype(sigmonial) = 4;
361                end
362                yalmip('setmonomtable',mt,[oldvariabletype newvariabletype],[mt_hash;new_mt_hash],oldhash);
363            end
364
365            if ~(x_isscalar | y_isscalar)
366                Z.dim(1) = X.dim(1);
367                Z.dim(2) = Y.dim(2);
368            else
369                Z.dim(1) = max(X.dim(1),Y.dim(1));
370                Z.dim(2) = max(X.dim(2),Y.dim(2));
371            end
372        catch
373            error(lasterr)
374        end
375        % Reset info about conic terms
376        Z.conicinfo = [0 0];
377        Z = clean(Z);
378
379    case 2
380
381        n_X = X.dim(1);
382        m_X = X.dim(2);
383        [n_Y,m_Y] = size(Y);
384
385        x_isscalar =  (n_X*m_X==1);
386        y_isscalar =  (n_Y*m_Y==1);
387
388        if ~x_isscalar
389            if ((m_X~= n_Y & ~y_isscalar))
390                error('Inner matrix dimensions must agree.')
391            end
392        end
393
394        n = n_X;
395        m = m_Y;
396        Z = X;
397
398        if x_isscalar
399            if y_isscalar
400                if Y==0
401                    Z = 0;
402                    return
403                else
404                    Z.basis = Z.basis*Y;
405                    % Reset info about conic terms
406                    Z.conicinfo = [0 0];
407                    return                               
408                end                   
409            else
410                Z.dim(1) = n_Y;
411                Z.dim(2) = m_Y;
412                Z.basis = kron(Z.basis,Y(:));
413                Z.conicinfo = [0 0];
414                Z = clean(Z);
415                return
416            end
417        elseif y_isscalar
418            Z.dim(1) = n_X;
419            Z.dim(2) = m_X;
420            Z.basis = Z.basis*Y;
421            Z.conicinfo = [0 0];
422            Z = clean(Z);
423            return
424        end
425
426        Z.dim(1) = n;
427        Z.dim(2) = m;
428        Z.basis = kron(Y.',speye(n_X))*X.basis;
429        Z.conicinfo = [0 0];
430        Z = clean(Z);
431
432    case 1
433
434        n_Y = Y.dim(1);
435        m_Y = Y.dim(2);
436        [n_X,m_X] = size(X);
437
438        x_isscalar =  (n_X*m_X==1);
439        y_isscalar =  (n_Y*m_Y==1);
440
441        if ~x_isscalar
442            if ((m_X~= n_Y & ~y_isscalar))
443                error('Inner matrix dimensions must agree.')
444            end
445        end
446
447        n = n_X;
448        m = m_Y;
449        Z = Y;
450
451        % Special cases
452        if x_isscalar
453            if y_isscalar
454                if X==0
455                    Z = 0;
456                    return
457                else
458                    Z.basis = Z.basis*X;
459                    Z.conicinfo = [0 0];
460                    return
461                end
462            else
463                Z.dim(1) = n_Y;
464                Z.dim(2) = m_Y;
465                try
466                    Z.basis = sparse(X)*Y.basis;
467                catch
468                    % This works better when low on memory in some cases
469                    [i,j,k] = find(Y.basis);
470                    Z.basis = sparse(i,j,X*k,size(Y.basis,1),size(Y.basis,2));
471                end
472                Z.conicinfo = [0 0];
473                Z = clean(Z);
474                return
475            end
476        elseif y_isscalar
477            Z.dim(1) = n_X;
478            Z.dim(2) = m_X;
479            Z.basis = X(:)*Y.basis;
480            Z = clean(Z);
481            return
482        end
483
484        if m_Y==1
485            Z.basis = X*Y.basis;
486        else
487            try               
488                speyemy = speye(m_Y);
489                kronX = kron(speyemy,X);
490                Z.basis = kronX*Y.basis;
491            catch
492                for i = 1:size(Y.basis,2);
493                    dummy = X*reshape(Y.basis(:,i),Y.dim(1),Y.dim(2));
494                    Z.basis(:,i) = dummy(:);
495                end
496            end
497        end
498        Z.dim(1) = n;
499        Z.dim(2) = m;
500        Z.conicinfo = [0 0];
501        Z = clean(Z);
502
503    otherwise
504        error('Logical error in mtimes. Report bug')
505end
506
507
508
509function Z=clean(X)
510temp = any(X.basis,1);
511temp = temp(2:end);
512index = find(temp);
513if ~isempty(index)
514    Z = X;
515    if length(index)~=length(Z.lmi_variables)
516        Z.basis = Z.basis(:,[1 1+index]);
517        Z.lmi_variables = Z.lmi_variables(index);
518    end
519else
520    Z = full(reshape(X.basis(:,1),X.dim(1),X.dim(2)));
521end
Note: See TracBrowser for help on using the repository browser.