update_supp

 script to update support

0001 % script to update support
0002 if length(gamma_xh)==M && flag == 1
0003     % swap an element...
0004     Ag = A(:,gamma_xh)'*A(:,idelta);
0005     xg = xk_1(gamma_xh);
0006     zx = -sign(pk(idelta)+delta*dk(idelta));
0007     
0008     % delx = -AtAgx\(Ag*zx);
0009     % delx = -AtAgx\(Ag*zx); % WHY?????????
0010     switch delx_mode
0011         case 'mil'
0012             delx = -iAtA*(Ag*zx);
0013         case 'qr'
0014             delx = -R\(R'\(Ag*zx));
0015     end
0016     
0017     delta_out = -xg./delx;
0018     idelta_out = find(delta_out>0);
0019     if ~isempty(idelta_out)
0020         [delta_t idelta_tmp] = min(delta_out(idelta_out));
0021         idelta_t = gamma_xh(idelta_out(idelta_tmp));
0022         
0023         x_old = xk_1;
0024         xk_1(gamma_xh) = xg+delta_t*delx;
0025         xk_1(idelta) = delta_t*zx;
0026         xk_1(idelta_t) = 0;
0027         
0028         % Remember remember the Simplex pivoting
0029         % disp(sprintf(['iter = %d, delta = %3.4g, delta_t = %3.4g, idelta = %d, idelta_t = %d, flag = %d'], iter, delta, delta_t, idelta, idelta_t, -2));
0030         
0031         outx_index = find(gamma_xh==idelta_t);
0032         gamma_xh = [gamma_xh(1:outx_index-1); gamma_xh(outx_index+1:end); gamma_xh(outx_index)];
0033         g_old = gamma_xh;
0034         
0035         gamma_xh(end) = idelta;
0036         
0037         % itr_history = [itr_history; idelta_t delta_t 2];
0038         
0039         new_order = [1:outx_index-1, outx_index+1:length(g_old), outx_index]';
0040         
0041         in_delx.new_order = new_order;
0042         in_delx.out_i = outx_index;
0043         in_delx.add_rem = 2;
0044         in_delx.Gamma = g_old;
0045         in_delx.new_x = idelta;
0046         idelta = idelta_t; % To set its constraints.
0047     else
0048         disp('Fatal error');
0049         break;
0050     end
0051 else
0052     if flag == 1
0053         % add an element
0054         g_old = gamma_xh;
0055         gamma_xh = [gamma_xh; idelta];
0056         xk_1(idelta) = 0;
0057         
0058         gamma_xc(gamma_xc==idelta) = [];
0059         
0060         in_delx.new_x = idelta;
0061         in_delx.add_rem = 1;
0062         in_delx.Gamma = g_old;
0063     else
0064         % remove an element
0065         outx_index = find(gamma_xh==idelta);
0066         len_gamma = length(gamma_xh);
0067         if len_gamma == 1
0068             break;
0069         end
0070         gamma_xh = [gamma_xh(1:outx_index-1); gamma_xh(outx_index+1:end); gamma_xh(outx_index)];
0071         g_old = gamma_xh;
0072         
0073         gamma_xh = gamma_xh(1:len_gamma-1);
0074         
0075         gamma_xc = [gamma_xc; idelta];
0076         
0077         xk_1(idelta) = 0;
0078         
0079         new_order = [1:outx_index-1, outx_index+1:length(g_old), outx_index]';
0080         
0081         in_delx.new_order = new_order;
0082         in_delx.out_i = outx_index;
0083         in_delx.add_rem = 0;
0084         in_delx.Gamma = g_old;
0085     end
0086 end
0087 
0088 % temp_gamma = zeros(N,1);
0089 % temp_gamma(gamma_xh) = gamma_xh;
0090 % gamma_xc = find([1:N]' ~= temp_gamma);