function [modelName,xx,ff,rc,totClockTime] = dynareFindSteadyState(varargin)
%function [xx,ff,rc]=dynareFindSteadyState(optional args)
%compute steady state generate initval assignments file(<modName>InitValFile)
%inputs(all args optional, order irrelevant)
%a logical value(true,false) determines the displayed 
%    amount of algorithm iteration info, 
%    default false, true to display copious amounts of iteration info
%a pair of doubles determines global search iteration limits
%    [maxTimeAllowed itmax] default [10 100]
%    max cpu secs for global search phase
%    max number of iterations for global search phase
%anRSpec -- "global" search region specification  
%    type help regionSpec for explanation 
%(note, dynareFindSteadyState ignores any  initval's in the .mod file.)
%outputs
%xx -- approximation to steady state
%ff -- result of evaluating model eqns at approx steady state(should be zeroes)
%rc -- return code 1 for success 0 for failure
%file(<modName>InitValFile)  -- initval assignments to add to
%  <modName>.mod file
%  prints output of running dynare's check(ignores any steady state where
%  check crashes, reports check results when available)
%$Id: dynareFindSteadyState.m,v 1.14 2007/04/27 18:27:20 m1gsa00 Exp $
global M_ oo_ 
modelName='not set';
xx='not set';
ff='not set';
rc='not set';
totClockTime='not set';

try
  theFssInput=fssInput(varargin{:});
  disp(sprintf('%s',char(theFssInput)));
catch
  error('dynareFindSteadyState:algortihm input specification failed')
end
showEm=getShowEm(theFssInput);
itmax=getItmax(theFssInput);
maxTimeAllowed=getMaxTimeAllowed(theFssInput);
anRSpec=getAnRSpec(theFssInput);
modelName=getModelName(theFssInput);
if(~isempty(modelName))
  eval(['dynare ' modelName])
  global M_ oo_ 
  if(~strcmp(M_.fname,modelName))
    error(['dynareFindSteady: request for >>' modelName '<< ineffective ' ...
           'active model still >>' M_.fname])
  end
else
  try
    modelName=M_.fname;
  catch
    error('dynareFindSteadyState: no active model, but no model name specified')
  end
end


rc=0;
timeBeginPgm=clock;
%   if exist([M_.fname '_steadystate.m'])%this is a noop need to delete
%    else
warning('OFF','all')
[xval,rc,manyf,manyx,manyc,manyFeas]=findSteadyState(@justRes,...
                                                  genBounds(anRSpec),...
                                                  showEm,itmax,maxTimeAllowed,...
                                                  [oo_.exo_steady_state;oo_.exo_det_steady_state]);
warning('ON','all')
%  end 
theTol=1e-8;theMaxIter=1000;

%display('validate starting point for newton')

%fVal=feval(@forNewtonSys, xval,[]);
%display('sorted prenewton')

%reportBadEqns(showEm,fVal);


%chkEncounterdBadVals(showEm,fVal);
%display('done validation')
%display('applying newton method')
%warning off MATLAB:nearlySingularMatrix
%warning off MATLAB:divideByZero
%[xx,ff,numit,rc]=newtonsys(@forNewtonSys,xval, theTol, theMaxIter,showEm);
%warning on MATLAB:nearlySingularMatrix
%warning on MATLAB:divideByZero
chkRes=0;
numOthers=size(manyf,2);
disp(sprintf('%i init vals to try',numOthers));
timeStart=clock;
for(ii=1:numOthers)
  if(manyc(ii)>0|manyFeas(ii)>0)
    %disp(sprintf('ii=%i not feasible, skipping',ii));
  else
    xval=manyx(:,ii);
chkRes=chkAGuess(ii,xval,theTol,theMaxIter,showEm,timeStart,numOthers);
rc=chkRes;
end
if(chkRes==1)
return;
end
end
timeEndPgm=clock;
totClockTime=etime(timeEndPgm,timeBeginPgm);
disp('failed to compute steady state')
disp(theFssInput);
unsuccessfulInitValFile([M_.fname,'InitValFile']);

function theRes=forNewtonSys(xx,jFlag)
global M_ oo_ 
if(isempty(jFlag))
  theRes=feval([M_.fname '_static'], xx, [oo_.exo_steady_state;oo_.exo_det_steady_state],M_.params);
else
  [ig,theRes]=feval([M_.fname '_static'], xx, [oo_.exo_steady_state;oo_.exo_det_steady_state],M_.params);
end

function [res,NASFlag]=justRes(yy,varargin)
global M_ oo_ 

try
  notRes=feval([M_.fname '_static'], yy, [oo_.exo_steady_state;oo_.exo_det_steady_state],M_.params);
  res= norm(notRes);
  if(isfinite(res))
    NASFlag=0;
  else
    %reportBadEqns(notRes);
    NASFlag=1;
  end
catch
  res=inf;NASFlag=1;
  if(exist('notRes'))
    %reportBadEqns(notRes);
  else
  end
end


function encounterGoodVals(fVal)
realSmallVal=1e-10
realSmall=find((abs(fVal)<=realSmallVal))
if(length(realSmall>0))
  realSmall
end



function chkEncounterdBadVals(showEm,fVal)
if(showEm)
  realBig=find((fVal==inf)+(fVal==-inf));
  realBad=find(isnan(fVal));
  if(length(realBig)>0)
    realBig
    fVal
    %error('infinity for some equations')
    disp('infinity for some equations')
  end
  if(length(realBad)>0)
    realBad
    disp('NaN for some equations: quitting newtonsys');
    %error('NaN for some equations: quitting newtonsys');
  end
end

function reportBadEqns(showEm,fVal)
if(showEm)
  realBig=find((fVal==inf)+(fVal==-inf));
  realBad=find(isnan(fVal));
  if(length(realBig)>0)
    'infinity for some equations'
    realBig
  end
  if(length(realBad)>0)
    'NaN for some equations'
    realBad
  end
  [toFindBig,idx]=sort(abs(fVal),1,'descend');
  aTenth=max(1,min(10,round(length(fVal)/10)));
  disp('largest discrepancies and their locations')
  disp(toFindBig(1:aTenth))
  disp(idx(1:aTenth))
end

function createInitValFile(fileName,varNames,initVals)
numVals=size(varNames,1);
if(and(and(isa(varNames,'char'),isa(initVals,'double')),isa(fileName,'char')))
  if(size(initVals,1)==numVals)
    [fid,msg]=fopen(fileName,'w');
    if(isempty(msg))
      fprintf(fid,'initval;\n');
      for(ii=1:numVals)
        fprintf(fid,'%s=%20.10f;\n',varNames(ii,:),initVals(ii));
      end
      fprintf(fid,'end;\n');
      fclose(fid);
    else
      fclose(fid);
      error(['createInitValFile: cannot open file >>%s<<', msg],fileName)
    end
  else
    fclose(fid);
    error('createInitValFile:number of values to assign must equal number of names')
  end
else
  fclose(fid);
  error('createInitValFile:expects char,char,double')
end



function unsuccessfulInitValFile(fileName)
[fid,msg]=fopen(fileName,'w');
if(isempty(msg))
  fprintf(fid,['steady state computation unsuccessful: ' msg]);
else
  error('unsuccessfulInitValFile:cannot open file >>%s<<',fileName);
end



function chkRes=chkAGuess(ii,xval,theTol,theMaxIter,showEm,timeStart,numOthers)
global M_ oo_ 
chkRes=0;
    %disp(sprintf('ii=%d of %d trying newton',ii,numOthers));
    warning off MATLAB:nearlySingularMatrix
    warning off MATLAB:singularMatrix
    warning off MATLAB:divideByZero
    warning off all
    [xx,ff,numit,rc]=newtonsys(@forNewtonSys,xval,...
                               theTol, theMaxIter,showEm);
    timeNow=clock;

    %      if(and(showEm,0==mod(ii,floor(numOthers/100))))%show about every 1 in four
    if(and(showEm,1))%show about every 1 in four
      disp(sprintf('time so far=%15.2f,remaining timeat present rate=%15.2f',...
                   etime(timeNow,timeStart),((numOthers/ii)-1)*(etime(timeNow,timeStart))));
    end
    try
      if(rc==1)
        oo_.steady_state=xx;
        %'do check'
        %       'not doing check'
        try
          chkRes=dynareAMACheck;
        catch
          rc=0;
          error(sprintf('dynareFindSteadyState: !!!!!!!!!!!!!!!!!check pgm crashes  for convergent steady state, ii=%d', ii))
        end
        disp(sprintf('ii=%d found start that converges,norm of f(x)=%d',ii,norm(ff)));
        if(chkRes==1)
          disp(sprintf(['also passes unique stable solution test']))
          disp(sprintf(['creating InitValFile']))
          createInitValFile([M_.fname,'InitValFile'],M_.endo_names,xx);
          timeEndPgm=clock;
          totClockTime=etime(timeEndPgm,timeBeginPgm);
          return
        else
        end
      else
      end
    catch
    end
    %disp('steady state not successfully computed with current start guess');
    warning on MATLAB:nearlySingularMatrix
    warning off MATLAB:singularMatrix
    warning on MATLAB:divideByZero
    warning on all