function approxBarMu=zeroCumulantsClosure(net,mdyn)
% Auxiliary m-script used by momentClosure.m to compute the moment closure
% functions obtained by setting high-order cumulants equal to zero.
%
% Copyright (C) 2006  Joao Hespanha

% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% Modification history:
%
% Created May 12, 2007
%
% Modified June 27, 2007
% Get moment dynamics as a single structure

verbose=0;

nMoments=size(mdyn.Mu.ndx,1);
nBarMu=size(mdyn.barMu.ndx,1);

nSpecies=size(mdyn.barMu.ndx,2);

nSpecies=length(net.species);
nReactions=length(net.reaction);

%% create symbolic species vector

symSpecies=[];
for thisSpecies=1:nSpecies
  symSpecies=[symSpecies,sym(net.species(thisSpecies).id)];
end

% Compute polynomial expansion of the moment generating function, up
% to the order of the largest moment in BarMu
maxorder=max(sum(mdyn.barMu.ndx,2));
syms lx
g=taylor(exp(sqrt(-1)*lx),lx,maxorder+1);
lx=0;
for i=1:nSpecies
   lx=lx+sym(sprintf('l%d',i))*symSpecies(i);
 end
%lx
g=sym(subs(g,'lx',lx,0));
[coef,mon,idx]=monomialindexes(g,symSpecies);

% take expected value
g=0;
for thisMon=1:size(idx,1)
  if any(idx(thisMon,:))
    g=g+coef(thisMon)*expname(net,idx(thisMon,:));
  else
    g=g+coef(thisMon); 
  end    
end  

% find cumulants that must be set to zero:
% . moments that appear in BarMu
% . moments of non-maximal order that appear in moment generating function
% but not in mdyn.Mu.ndx

k=find(sum(idx,2)<maxorder);
zerocum=[mdyn.barMu.ndx;idx(k,:)];
for thisMon=size(zerocum,1):-1:1
  % monomial in state?    
  k=find(all(mdyn.Mu.ndx==ones(nMoments,1)*zerocum(thisMon,:),2));
  if ~any(zerocum(thisMon,:)) || ~isempty(k)
    zerocum(thisMon,:)=[];
  end    
end
%zerocum

% compute cummulants that must be set to zero
equation='';
unknown='';
for thisCum=1:size(zerocum,1)
  cg=log(g);
  for i=1:nSpecies
    cg=diff(cg,sprintf('l%d',i),zerocum(thisCum,i));
  end
%  cg
%  [sprintf('l%d=0,',1:nSpecies-1),sprintf('l%d=0',nSpecies)]
  cum=maple('subs',[sprintf('l%d=0,',1:nSpecies-1),sprintf('l%d=0',nSpecies)],cg);
  equation=[equation,',',char(cum)];
  unknown=[unknown,',',char(expname(net, zerocum(thisCum,:)))];
end  
if length(equation)>1
  equation(1)=[];
  unknown(1)=[];
end  

% find moments
sol=solve(equation,unknown);

% extract moment to output

approxBarMu=sym(zeros(size(mdyn.barMu.ndx,1),1));

for thisMon=1:size(mdyn.barMu.ndx,1)
  if any(mdyn.barMu.ndx(thisMon,:))
    if verbose  
      fprintf('cumulantsClosure: closing %s ~ ',char(expname(net,mdyn.barMu.ndx(thisMon,:))))
    end
    
    if isstruct(sol)
      approxBarMu(thisMon)=expand(getfield(sol,char(expname(net,mdyn.barMu.ndx(thisMon,:)))));
    else % single unknown?
      approxBarMu(thisMon)=expand(sol);
    end  
    
    if verbose  
      fprintf('%s\n',char(approxBarMu(thisMon)))
    end    
  else
    approxBarMu(thisMon)=1;
  end    
end

%mdyn.barMu.ndx
%approxBarMu