function m=saturate_up(m,max_y,j);

[in,out]=get_regressors(m);
range_in = max(in{j})-min(in{j});
mean_in = mean(in{j});
dyn=get(m,'dynamics');
n_in=size(in{j},2);
n_y=sum(dyn.ny(j,:));

n_u=n_in-n_y;

%Build a dataset for yk, yk-1, ... = max_y
points=generate(40,[min(in{j}(:,n_y+1:end));max(in{j}(:,n_y+1:end))]);
reg{j}=[max_y*ones(size(points,1),n_y) points];
outs=eval(m,reg,j);

%just keep points that are bigger than the limit
ids=find(outs>max_y);
sub_in=points(ids,:);
sub_out=outs(ids)-max_y;

%now build a tak_sug model with these data
n_r=2^n_u;

c_m=taksug('cancelator',n_u,1);
c_m=add_data(c_m,sub_in,sub_out);
c_m=set_static(c_m);
c_m=identify(c_m,'cluslms',struct('n_rules',n_r));

%take the symetrical of the rules
c_m.description{1}.linears=-c_m.description{1}.linears;

% do the offset
c_m.description{1}.linears(:,end,:)=c_m.description{1}.linears(:,end,:)+max_y;

% Now we have to go back to the full dimensionality
f_c=[(range_in(1)/2+max_y)*ones(n_y,n_r);c_m.description{1}.centers];

ref_pos=[max_y*ones(n_r,n_y) c_m.description{1}.centers'];
ref_memb=min(sum(membership(m,ref_pos,j)',1),0.8);

dist=n_y*(range_in(1)/2)^2;

switch m.description{j}.model_code{2},
	case 'gaussian',
		v= -log(ref_memb)/dist;
		m=add_rules(m,n_r,{'productspace' 'gaussian' 'linear'},j);
		
	case 'inversedist',
		v= ref_memb.^(1-m.description{j}.m)/dist;
		m=add_rules(m,n_r,{'productspace' 'inversedist' 'linear'},j);
end

for i=1:n_r,
	f_v(:,:,i)=[diag(v(i)*ones(n_y,1)) zeros(n_y,n_u);zeros(n_u,n_y) c_m.description{1}.ivariances(:,:,i)];
	f_l(:,:,i)=[zeros(1,n_y) c_m.description{1}.linears(:,:,i)];
end

m.description{j}.centers(:,end-n_r+1:end)=f_c;
m.description{j}.ivariances(:,:,end-n_r+1:end)=f_v;
m.description{j}.linears(:,:,end-n_r+1:end)=f_l;

% Now if we want, we can tune the linears in order to optimise the simulation
opt=foptions;
opt(1)=1;
opt(14)=100;
x=fmins('loc_err',[0.01*max_y;-0.01],opt,[],m,j,n_r);

for i=n_r-1:-1:0,
	m.description{j}.linears(:,1,end-i)=x(2);
	m.description{j}.linears(:,end,end-i)=m.description{j}.linears(:,end,end-i)+x(1);
end


% Generate the data
function r=generate(d,min_max)

n=size(min_max,2);
steps=(min_max(2,:)-min_max(1,:))/(d-1);
num=d^n;
siz=1;
r=zeros(num,n);
i=1;
while num>1,
	num=num/d;
	range=(min_max(1,i):steps(i):min_max(2,i))';
	generator=range(:,ones(siz,1))';
	r(:,i)=repmat(generator(:),[num 1]);
	siz=siz*d;
	i=i+1;
end

