feed_irreg.m

software to find the threshold of an LDPC code in a Gaussian channel

function [y,dec_fail] = feed_irreg(rho,lambda,f_in,K,delta,lpad,iter)

% FEED = Fast Enough Evolution of Densities

y = zeros(1,iter);

K_pow2 = 2^(ceil(log2(K)));
x0 = 0:delta:(K*delta);
v_varm = exp((-1/2).*x0);
p_infc = 0;

f_in_pad = zeros(1,K_pow2);
f_in_pad(1:(K+1)) = f_in((K+1):(2*K+1)).*v_varm;
f_in_pad((2*K_pow2-K+1):(2*K_pow2)) = f_in_pad((K+1):-1:2);
F_f_in = real(fft(f_in_pad));

fp = f_in((K+1):(2*K+1));
fn = f_in((K+1):-1:1);
fp(1) = fp(1)/2;
fn(1) = fn(1)/2;

A = zeros(K+2*lpad+1,K);
x = delta:delta:(K*delta);

for zi = (-1*lpad*delta):delta:((K+lpad)*delta)
	z_index = round(zi/delta)+lpad+1;
	A(z_index,:) = tanh(x/2).^exp(zi);
end

v_varm = exp((-1/2).*x0);

% check node

nnls_opts = optimset('TolX',1e-14); % this makes lsqnonneg happy

ic = 1;
keep_going = 1;

while((ic <= iter) & (keep_going))

p_infc_last = p_infc;

[A_r,map_r,K_r] = active_rgn(fp,A,delta,lpad);

G1 = 1 - test_tanh2(delta,K_r,fp(1:(K_r+1))+fn(1:(K_r+1)),lpad);
H1 = 1 - test_tanh2(delta,K_r,fp(1:(K_r+1))-fn(1:(K_r+1)),lpad);
H1 = H1 - H1(length(H1)) + G1(length(G1));
% this will guarantee that the "infinity component" in H is zero

G = (G1+H1)/2;
H = (G1-H1)/2;



%keyboard

Gp = zeros(1,length(G));
Gn = zeros(1,length(G));

for dc = 2:length(rho)

	if (rho(dc) > 0)

		for c = 0:(dc-1)

		        % c represents the number of negative messages

			if (mod(c,2) == 0)
				% the result is positive
				Gp = Gp + rho(dc)*nchoosek(dc-1,c)*((H.^c).*(G.^(dc-1-c)));
			else
				Gn = Gn + rho(dc)*nchoosek(dc-1,c)*((H.^c).*(G.^(dc-1-c)));
			end

		end

	end

end

alphap = zeros(1,K);
alphan = zeros(1,K);

p_infc = Gp(length(Gp)); % this is the component of the probability mass "at infinity"
Gp = Gp - p_infc; % remove it from the Fourier transform ...

% decimation -- to keep lsqnonneg happy

if (K_r < 50)
	alphapr = lsqnonneg(A_r,Gp',[],nnls_opts);
	alphap(1:K_r) = alphapr;
elseif (K_r < 75)
        alphapr = lsqnonneg(A_r(:,1:2:K_r),Gp',[],nnls_opts);
        alphap(1:2:K_r) = alphapr;
else
	alphapr = lsqnonneg(A_r(:,1:3:K_r),Gp',[],nnls_opts);
	alphap(1:3:K_r) = alphapr;
end

%alphapr = lsqnonneg(A_r,Gp',[],nnls_opts);
%alphanr = lsqnonneg(A_r,Gn',[],nnls_opts);


% construct negatives and normalize
f_x = zeros(1,2*K+1);
f_x((K+2):(2*K+1)) = alphap;
f_x(K:-1:1) = alphap.*exp(-1*x);
f_x(K+1) = 1-sum(f_x)-p_infc;
p_zero = f_x(K+1);
if (f_x(K+1) < 0)
        f_x = f_x./(sum(f_x) - f_x(K+1));
        f_x(K+1) = 0;
end

% p_infc

% variable nodes

f_var = zeros(1,2*K_pow2);
f_var(1:(K+1)) = f_x((K+1):(2*K+1)).*v_varm;
f_var((2*K_pow2-K+1):(2*K_pow2)) = f_var((K+1):-1:2);

F_f_out = zeros(1,2*K_pow2);
F_f_var = real(fft(f_var));

p_infv = 0;

for dv = 2:length(lambda)

	if (lambda(dv) > 0)
		F_f_out = F_f_out + lambda(dv)*(F_f_in.*(F_f_var.^(dv-1)));
		p_infv = p_infv + lambda(dv)*(1-(1-p_infc)^(dv-1));
	end

end

iF_f_var = real(ifft(F_f_out));



fp = iF_f_var(1:(K+1));
fn = fp .* v_varm; % by the symmetry property
fp = fp ./ v_varm;
fp(1) = fp(1)/2;
fn(1) = fn(1)/2;
s = sum(fp) + sum(fn) + p_infv;
%fp = fp./s;
%fn = fn./s;

[ic sum(fn) s p_infc p_infv p_zero]
%sum(fn)

y(ic) = sum(fn);


% stopping conditions
%

%if ((ic > 5) & (abs(y(ic)-y(ic-1)) < 1e-5) & (abs(p_infc_last - p_infc) < 1e-5))
%	keep_going = 0; % force termination
%else
if (sum(fn) < 1e-5)
	keep_going = 0;
else
	keep_going = 1;
end

%keyboard

ic = ic + 1;

end

if (y(ic-1) > 1e-5)
	dec_fail=1;
else
	dec_fail=0;
end