function c=init_controller(m,pert_u);
% Works for 1 output only
% Determine the center of the space

n_out=get(m,'n_out');
n_in=get(m,'n_in');

in=get(m,'data_in');
out=get(m,'data_out');
y_r=mean(m.famimo.data.bds_out);
u_r=mean(m.famimo.data.bds_in);
u_off=(m.famimo.data.bds_in-u_r([1 1],:))*0.6;
dyn=get(m,'dynamics');

mask=build_regressors(ones(max(max(max(dyn.ny+(dyn.ny>0))),max(max(dyn.nu+dyn.nd))),1)*[1:n_in+n_out],dyn.ny,dyn.nu,dyn.nd,1);
u_r = fmins('taksug_out',u_r,[],[],y_r,mask,m);

% This is a variant of fit_linears design for the case where there is no offset
% Now if you have many output you are in the shit because you do not need to have
% the same partition (the same rules positions) for each output

n_y_reg=dyn.ny;
n_u_reg=sum(dyn.nu);
tm = add_data(m,in - u_r(ones(size(in,1),1),:),out - y_r(ones(size(in,1),1),:));

j=1;

[xx,yy]=get_regressors(tm);

x=xx{j};
y=yy{j};
	
n_rules=m.description{j}.n_rules;

members=membership(tm,x,j)';

sum_memb=sum(members,1)+realmin;

alphas=(members./sum_memb(ones(1,n_rules),:))';

i1=1:n_rules;
i2=i1(ones(1,n_u_reg+n_y_reg),:);
i3=(1:n_u_reg+n_y_reg)';
i4=i3(:,ones(1,n_rules));

theA=alphas(:,i2).*x(:,i4);

theB=y;
theL=theA\theB;
l = reshape(theL,1,n_u_reg+n_y_reg,n_rules);

% Build the controller
c_n_in = 2 + size(pert_u,2);
c_n_out = n_in - size(pert_u,2);

cont_u=setdiff(1:n_in,pert_u);
x_ids = [1:n_y_reg, n_y_reg+pert_u];
u_ids = n_y_reg+cont_u;

c = taksug([m.famimo.name '_controller'], c_n_in, c_n_out);
c = add_dynamics(c,zeros(c_n_out),[n_y_reg * ones(c_n_out,2) ones(c_n_out, c_n_in - 2)],[zeros(c_n_out,2) zeros(c_n_out, c_n_in - 2)]);
c = add_rules(c,n_rules);
set(c,{1:c_n_out 'centers'},[m.description{j}.centers(1:n_y_reg,:); m.description{j}.centers(x_ids,:)]);

for i=1:n_rules,
	
	A_i = [l(:,x_ids,i);eye(size(x_ids,2)-1) zeros(size(x_ids,2)-1,1)];
	B_i = [l(:,u_ids,i); zeros(size(x_ids,2)-1,c_n_out)];
	
	K_i = dlqr(A_i,B_i,eye(c_n_in),0.2*eye(c_n_out));
	
	for j_c=1:c_n_out,
		v = projclus(m.description{j}.ivariances(:,:,i),(x_ids)');
		c.description{j_c}.ivariances(:,:,i) = [ones(2)*1e-3 zeros(2,size(v,2)); zeros(size(v,1),2) v];
		c.description{j_c}.linears(:,:,i) = [-K_i(j_c,1:2) K_i(j_c,:) u_r(cont_u(j_c))-K_i(j_c,:)*[zeros(n_y_reg,1);u_r(pert_u)]];
	end
end

	



