source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/extras/monolist.m @ 37

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

Added original make3d
File size: 6.6 KB
Line 
1function new_x = monolist(x,d,min)
2% MONOLIST Generate monomials
3%
4% y = MONOLIST(x,d,f)
5%
6% Returns the monomials [1 x(1) x(1)^2 ... x(1)^d(1) x(2) x(1)x(2)  etc...]
7%
8% do
9%
10% >>sdpvar x y z
11% >>sdisplay(monolist([x y z],4))
12%
13% to see the logic!
14%
15%  Input
16%     x      : Vector with SDPVAR variables
17%     d      : Integers > 0
18%
19%  See also SDPVAR, DEGREE
20
21% Author Johan Löfberg
22% $Id: monolist.m,v 1.25 2006/08/28 14:03:29 joloef Exp $
23
24% FIX : Optimal generation
25
26% Flatten
27x = reshape(x,1,length(x));
28x_orig = x;
29
30if (length(d)==1 | all(d(1)==d)) & islinear(x) & ~isa(x,'ncvar')
31    d = d(1);
32    % Powers in monomials
33    powers = monpowers(length(x),d);
34   
35    % Use fast method for x^alpha
36    if isequal(getbase(x),[zeros(length(x),1) eye(length(x))])
37        new_x = recovermonoms(powers,x);
38        return
39    end
40
41    % Monolist on dense vectors is currently extremely slow, but also
42    % needed in some applications (stability analysis using SOS) For
43    % performance issue, the code below is hard-coded for special cases
44    % FIX : Urgent, find underlying indexing...
45   
46    % Vectorize quadratic and quadrtic case
47    if d==2 & length(x)>1
48        V=x.'*[1 x];
49        ind=funkyindicies(length(x));
50        new_x = [1 V(ind(:)).'].';
51        return
52    elseif (length(x)==4 & d==6)
53       
54         ind =[    1           2           3           4           5           6           7           8,
55
56           9          10          11          12          13          14          15          16,
57
58          17          18          19          20          21          22          23          24,
59
60          25          26          27          28          29          30          31          32,
61
62          33          34          49          50          51          52          86          53,
63
64          54          55          89          56          57          91          58          92,
65
66         126          59          60          61          95          62          63          97,
67
68          64          98         132          65          66         100          67         101,
69
70
71         135          68         102         136         170         185         186         187,
72
73         188         222         256         189         190         191         225         259,
74
75
76         192         193         227         261         194         228         262         296,
77
78         330         364         195         196         197         231         265         198,
79
80         199         233         267         200         234         268         302         336,
81
82         370         201         202         236         270         203         237         271,
83
84         305         339         373         204         238         272         306         340,
85
86         374         408         442         476         510         525         526         527,
87
88         528         562         596         630         529         530         531         565,
89
90
91         599         633         532         533         567         601         635         534,
92
93
94         568         602         636         670         704         738         772         806,
95
96
97         840         535         536         537         571         605         639         538,
98
99         539         573         607         641         540         574         608         642,
100
101         676         710         744         778         812         846         541         542,
102
103         576         610         644         543         577         611         645         679,
104
105         713         747         781         815         849         544         578         612,
106
107
108         646         680         714         748         782         816         850         884,
109
110
111         918         952         986        1020        1054        1088        1122        1156,
112
113        1190          0            0          0          0             0           0            0];
114        ind = ind';
115        ind = ind(find(ind));
116        v=monolist(x,3);
117        V=v(2:end)*v.';
118        new_x = [1;V(ind(:))];
119        return
120    elseif d==4 & (1<=length(x)) & length(x)<=4 %& length(x)>1
121        v=monolist(x,2);
122        V=v(2:end)*v.';
123
124        % Cone to generate indicies
125        %p = sdpvar(n,1);
126        %v = monolist(p,2);
127        %V = v(2:end)*v';V=V(:);
128        %m = monolist(p,4)
129        %ind = [];
130        %for i = 2:length(m)
131        % ind = [ind min(find(~any(V-m(i))))];
132        %end
133
134        switch length(x)
135            case 1
136                new_x = [1; V([1 2 4 6]')];
137                return
138            case 2
139                new_x = [1;V([1 2 3 4 5 8 9 10 15 18 19 20 25 30]')];
140                return;
141            case 3
142                new_x=[1;V([1 2 3 4 5 6 7 8 9 13 14 15 24 16 17 26 18 27 36 40 41 42 51 60 43 44 53 62 45 54 63 72 81 90]')];
143                return
144            case 4
145                new_x=[1;V([    1     2     3     4     5     6     7     8     9    10    11    12    13    14    19    20    21    35    22    23    37    24    38    52    25    26    40    27    41    55    28    42 56    70    75    76    77    91   105    78    79    93   107    80    94   108   122   136  150    81    82    96   110    83    97   111   125   139   153    84    98   112   126   140  154   168   182   196   210]')];
146                return
147            otherwise
148        end
149    end
150
151
152    % Na, we have things like (c'x)^alpha
153    % precalc x^p
154    for i = 1:length(x)
155        temp = x(i);
156        precalc{i,1} = temp;
157        for j = 2:1:d
158            temp = temp*x(i);
159            precalc{i,j} = temp;
160        end
161    end
162
163    new_x = [];
164
165    for i = 1:size(powers,1) % All monomials
166        temp = 1;
167
168        for j = 1:size(powers,2) % All variables
169            if powers(i,j)>0
170                temp = temp*precalc{j,powers(i,j)};
171            end
172        end
173        new_x = [new_x temp];
174
175    end
176
177else
178
179    d = d(:)*ones(1,length(x_orig));
180
181    x = [1 x];
182
183    % Lame loop to generate all combinations
184    new_x = 1;
185    for j = 1:1:max(d)
186        temp = [];
187        for i = 1:length(x)
188            temp = [temp x(i)*new_x];
189        end
190        new_x = temp;
191        new_x = fliplr(unique(new_x));
192        new_degrees = degree(new_x,x(2:end));
193        remv = [];
194        for i = 1:length(d);
195            if new_degrees(i)>=d(i);
196                x = recover(setdiff(getvariables(x),getvariables(x_orig(i))));
197                x = [1;x(:)];
198                remv = [remv i];
199            end
200        end
201        d = d(setdiff(1:length(d),remv));
202    end
203end
204new_x = reshape(new_x(:),length(new_x),1);
205
206function ind = funkyindicies(n)
207
208M=reshape(1:n*(n+1),n,n+1);
209ind = M(:,1)';
210for i = 1:n
211    ind = [ind M(i,2:i+1)];
212end
213
214
Note: See TracBrowser for help on using the repository browser.