function [n_cand,Cand_models]=gen_mv(Model,InpPat,order);

n_cand=0;
n_s_b=Model.n_s_m;
[n_pat,n_inp]=size(InpPat);

%determine current complexity
comp=0;
for i=1:n_s_b
   comp=comp+sum(Model.Sub_Model(i).a_w);
end

add_comp=n_pat-comp;

% generate all possible combinations
for i=1:n_s_b
   for j=i+1:n_s_b
      var_1=Model.Sub_Model(i).var;
      var_2=Model.Sub_Model(j).var;
      new_var=var_1;
      com_var=[];
      for k=1:length(var_2)
         if isempty((find(var_2(k)==var_1)))
            new_var=[new_var var_2(k)];
         else
            com_var=[com_var k];
         end
      end
      
      % new_var is the new multidemsional spline; see if it is already there
      found=0;
      k=1;
      while (~found)& (k<=Model.n_s_m)
         if length(Model.Sub_Model(k).var)==length(new_var)
            found=1;
            for l=1:length(new_var)
               found=found*(~isempty(find(Model.Sub_Model(k).var==new_var(l))));
            end
         end
         k=k+1;
      end
      
      if ~found
         %new candidate; see if we still can manage it
         new_comp=prod(Model.Sub_Model(i).n_fun)*prod(Model.Sub_Model(j).n_fun);
         if new_comp<=add_comp
         	n_cand=n_cand+1;
         
         	Cand_models(n_cand)=Model;

				S_M=Model.Sub_Model(i);
            [m,n]=size(S_M.Knots);
            [m1,n1]=size(Model.Sub_Model(j).Knots);
            if n>n1
               Model.Sub_Model(j).Knots=[Model.Sub_Model(j).Knots zeros(m1,n-n1)];
            else
               if n1>n
                  S_M.Knots=[S_M.Knots zeros(m,n1-n)];
               end   
            end
            
            for l=1:m1
               if ~isempty(com_var)
                 	if isempty(length(find(com_var==l)))
      					S_M.Knots=[S_M.Knots; Model.Sub_Model(j).Knots(l,:)];
      					S_M.n_Knots=[S_M.n_Knots Model.Sub_Model(j).n_Knots(l)];  
                  	S_M.n_fun=[S_M.n_fun Model.Sub_Model(j).n_fun(l)];
                  	S_M.x_min=[S_M.x_min Model.Sub_Model(j).x_min(l)];
      					S_M.x_max=[S_M.x_max Model.Sub_Model(j).x_max(l)];
      					S_M.order=[S_M.order Model.Sub_Model(j).order(l)];
                     S_M.int_knots=[S_M.int_knots Model.Sub_Model(j).int_knots(l)];
                     S_M.w=zeros(length(S_M.w)*Model.Sub_Model(j).n_fun(l),1);
                     S_M.a_w=ones(length(S_M.w),1);

                  end
               else
                  S_M.Knots=[S_M.Knots; Model.Sub_Model(j).Knots(l,:)];
      				S_M.n_Knots=[S_M.n_Knots Model.Sub_Model(j).n_Knots(l)];  
                  S_M.n_fun=[S_M.n_fun Model.Sub_Model(j).n_fun(l)];
                  S_M.x_min=[S_M.x_min Model.Sub_Model(j).x_min(l)];
      				S_M.x_max=[S_M.x_max Model.Sub_Model(j).x_max(l)];
      				S_M.order=[S_M.order Model.Sub_Model(j).order(l)];
      				S_M.int_knots=[S_M.int_knots Model.Sub_Model(j).int_knots(l)];
                  S_M.w=zeros(length(S_M.w)*Model.Sub_Model(j).n_fun(l),1);
                  S_M.a_w=ones(length(S_M.w),1);

               end   
            end   
            
            S_M.var=new_var;

				Cand_models(n_cand).n_s_m=Cand_models(n_cand).n_s_m+1;
            Cand_models(n_cand).Sub_Model(Cand_models(n_cand).n_s_m)=S_M;
         else
            disp(['This candidate -' num2str(new_var) ' - would generate ' num2str(comp+new_comp) ...
                  ' complexity. Not enough training data!']);
         end
         
      end   
      
   end
end

if n_cand==0
   Cand_models=[];
end