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

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

Added original make3d
File size: 7.3 KB
Line 
1function y = subsasgn(X,I,Y)
2%SUBASGN (overloaded)
3
4% Author Johan Löfberg
5% $Id: subsasgn.m,v 1.5 2006/07/26 20:17:58 joloef Exp $   
6
7try
8    if strcmp('()',I.type)
9        X_is_spdvar = isa(X,'sdpvar');
10        Y_is_spdvar = isa(Y,'sdpvar');
11        if any(I.subs{1} <=0)
12            error('Index into matrix is negative or zero.');
13        end
14        switch 2*X_is_spdvar+Y_is_spdvar
15            case 1
16                % This code does not work properly
17                % Only work if b is undefined!!?!!
18                % generally ugly code...
19                y = Y;
20                [n_y,m_y] = size(Y);
21                y_lmi_variables = y.lmi_variables;
22                try
23                    X0 = sparse(subsasgn(full(X),I,full(reshape(Y.basis(:,1),n_y,m_y))));
24                    [n_x,m_x] = size(X0);
25                    y.basis = reshape(X0,n_x*m_x,1);
26                    X = full(X)*0;
27                    for i = 1:length(y_lmi_variables)
28                        X0 = full(sparse(subsasgn(X,I,full(reshape(Y.basis(:,i+1),n_y,m_y)))));
29                        y.basis(:,i+1) = reshape(X0,n_x*m_x,1);
30                    end
31                    y.dim(1) = n_x;
32                    y.dim(2) = m_x;
33                    % Reset info about conic terms
34                    y.conicinfo = [0 0];
35                catch
36                    error(lasterr)
37                end
38            case 2
39                if ~isempty(Y)
40                  Y = sparse(Y);
41                end
42                y = X;
43               
44                % Special code for speed
45                % elements in vector replaced with constants
46                if min(X.dim(1),X.dim(2))==1 & (length(I.subs)==1)
47                     y = X;
48                     y.basis(I.subs{1},1) = Y;
49                     y.basis(I.subs{1},2:end) = 0;                     
50                     y = clean(y);
51                     % Reset info about conic terms
52                     if isa(y,'sdpvar')
53                         y.conicinfo = [0 0];
54                     end
55                     return;
56                end
57                   
58               
59                x_lmi_variables = X.lmi_variables;
60                lmi_variables = [];
61               
62                % y.basis = [];
63                n = y.dim(1);
64                m = y.dim(2);
65                subX = sparse(subsasgn(full(reshape(X.basis(:,1),n,m)),I,Y));
66                y.basis = subX(:);
67               
68                j = 1;
69                Z = 0*Y;
70                for i = 1:length(x_lmi_variables)
71                    subX = sparse(subsasgn(full(reshape(X.basis(:,i+1),n,m)),I,Z));
72                    if (norm(subX,inf)>0)
73                        y.basis(:,j+1) = subX(:);
74                        lmi_variables = [lmi_variables x_lmi_variables(i)];
75                        j = j+1;
76                    end
77                end 
78                y.dim(1) = size(subX,1);
79                y.dim(2) = size(subX,2);
80                if isempty(lmi_variables) % Convert back to double!!
81                    y=full(reshape(y.basis(:,1),y.dim(1),y.dim(2)));
82                    return
83                else %Nope, still a sdpvar
84                    y.lmi_variables = lmi_variables;
85                     % Reset info about conic terms
86                    y.conicinfo = [0 0];
87                end
88               
89            case 3
90                z = X;
91               
92                x_lmi_variables = X.lmi_variables;
93                y_lmi_variables = Y.lmi_variables;
94               
95                               
96                % In a first run, we fix the constant term and null terms in the X basis
97                lmi_variables = [];
98                nx = X.dim(1);
99                mx = X.dim(2);
100                ny = Y.dim(1);
101                my = Y.dim(2);
102               
103                if (mx==1) & (my == 1) & isempty(setdiff(y_lmi_variables,x_lmi_variables)) & (max(I.subs{1}) < nx);
104                    % Fast specialized code for Didier
105                     y = specialcode(X,Y,I);
106                     return
107                end
108                %      subX = sparse(subsasgn(full(reshape(X.basis(:,1),nx,mx)),I,reshape(Y.basis(:,1),ny,my)));
109                subX = subsasgn(reshape(X.basis(:,1),nx,mx),I,reshape(Y.basis(:,1),ny,my));
110               
111                z.basis = subX(:);
112                j = 1;
113                yz = 0*reshape(Y.basis(:,1),ny,my);
114                for i = 1:length(x_lmi_variables)
115                    %        subX = sparse(subsasgn(full(reshape(X.basis(:,i+1),nx,mx)),I,yz));
116                    subX = subsasgn(reshape(X.basis(:,i+1),nx,mx),I,yz);
117                   
118                    if (norm(subX,inf)>0)
119                        z.basis(:,j+1) = subX(:);
120                        lmi_variables = [lmi_variables x_lmi_variables(i)];
121                        j = j+1;
122                    end
123                end
124                z.lmi_variables=lmi_variables;
125                all_lmi_variables = union(lmi_variables,y_lmi_variables);
126                in_z = ismembc(all_lmi_variables,lmi_variables);
127                in_y = ismembc(all_lmi_variables,y_lmi_variables);
128                z_ind = 2;
129                y_ind = 2;
130                basis=z.basis(:,1);
131                nz = size(subX,1);
132                mz = size(subX,2);
133                for i = 1:length(all_lmi_variables)
134                    switch 2*in_y(i)+in_z(i)
135                        case 1
136                            basis(:,i+1) = z.basis(:,z_ind);z_ind = z_ind+1;
137                        case 2
138                            temp = sparse(subsasgn(full(0*reshape(X.basis(:,1),nx,mx)),I,full(reshape(Y.basis(:,y_ind),ny,my))));
139                            basis(:,i+1) = temp(:);
140                            y_ind = y_ind+1;
141                        case 3
142                            Z1 = z.basis(:,z_ind);
143                            Z4 = Y.basis(:,y_ind);
144                            Z3 = reshape(Z4,ny,my);
145                            Z2 = sparse(subsasgn(0*reshape(full(X.basis(:,1)),nx,mx),I,Z3));
146                            temp = reshape(Z1,nz,mz)+Z2;                           
147%                            temp = reshape(z.basis(:,z_ind),nz,mz)+sparse(subsasgn(0*reshape(full(X.basis(:,1)),nx,mx),I,reshape(Y.basis(:,y_ind),ny,my)));
148                            basis(:,i+1) = temp(:);
149                            z_ind = z_ind+1;
150                            y_ind = y_ind+1;
151                        otherwise
152                    end
153                end;
154                z.dim(1) = nz;
155                z.dim(2) = mz;
156                z.basis = basis;
157                z.lmi_variables = all_lmi_variables;
158                y = z;                 
159                % Reset info about conic terms
160                y.conicinfo = [0 0];                 
161            otherwise
162        end
163    else
164        error('Reference type not supported');
165    end
166   
167catch
168    error(lasterr)
169end
170
171
172function y = specialcode(X,Y,I)
173
174y = X;
175X_basis = X.basis;
176Y_basis = Y.basis;
177ind = I.subs{1};ind = ind(:);
178yvar_in_xvar = zeros(length(Y.lmi_variables),1);
179for i = 1:length(Y.lmi_variables);
180    yvar_in_xvar(i) = find(X.lmi_variables==Y.lmi_variables(i));
181end
182y.basis(ind,:) = 0;
183mapper = [1 1+yvar_in_xvar(:)'];mapper = mapper(:);
184[i,j,k] = find(y.basis);
185[ib,jb,kb] = find(Y_basis);
186i = [i(:);ind(ib(:))];
187j = [j(:);mapper(jb(:))];
188k = [k(:);kb(:)];
189y.basis = sparse(i,j,k,size(y.basis,1),size(y.basis,2));
190y = clean(y);
191
192
Note: See TracBrowser for help on using the repository browser.