%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a slightly modified version of markov.g program that allows
% for reading parameters from the configuration file.
% Version 18 uses the Allw x nfirms coding method, and allows for starting
% at any number of firms.
% Written by: Gautam Gowrisankaran
% May 16, 1993
%
% Modified by Ron Goettler April 2011 to solve FOC within each state each 
% iteration. Modification consists of adding a function to return foc for  
% 1) entry prob  2) each firm's x  3) each firm's exit prob. The foc are 
% obtained simply by subtracting the candidate values from their optimal 
% values as determined within optimize().
%
% To run in Iterative Best Response mode (i.e., standard PM), the code calls  
% get_foc()  once, but doesn't call the Newton solver. One minor difference 
% between this iterative BR and PM is that  isentry  is updated along with 
% other policies, not beforehand as in PM. Hence, RLG added  oldentry. The 
% variable  isentry  is retained, instead of being renamed to  newentry.
% If needed, the old PM timing of updating  isentry  prior to updating 
% investments/exit could be reinstated without too much coding. But the foc 
% approach requires solving for entry simultaneously with investment/exit, 
% so cleaner to update entry similarly when not using foc.

% When both FOC and non-FOC converge, they appear to converge to the same 
% point.  But sometimes FOC converged when non-FOC doesn't, and vice-versa.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eql_ma_foc(batchnum)

disp('Entering eql_ma_foc.m ...'); 
  
% Build Par struct with general program data
Par.batchjob        = batchnum;
Par.configfile      = ['config_' num2str(Par.batchjob) '.mat'];
Par.nocheck         = -1e50; 
Par.tol             = 1e-4;  % Tolerance for convergence of value function ( 1e-13 for publishable results)
Par.foc_tol         = 1e-7;  % Tolerance for convergence of FOC solution within each state
Par.uu_competition  = 1;
Par.uu_monopoly     = 2;
Par.uu_planner      = 3;
Par.uu_quality      = 1;
Par.uu_cost         = 2;
Par.uu_capacity     = 3;
Par.uu_ran_entry    = 1;
Par.uu_det_entry    = 2;


% Load Params struct from file (if it exists) or abort
if exist(Par.configfile, 'file')
    disp('loading config file...');
    load(Par.configfile, 'Params');
else
    disp(['Could not find config file ' Par.configfile]);
    disp('Aborting!');
    return;
end

% Display model type
[prefix,eqltype,invtype] = acronym(Par, Params.eql_type, Params.ind_type);
disp(['model: ' eqltype ', investment in ' invtype]);

% Set a global for counting number of time we call the calcval function
% RP - disabled to not hamper parallelization
% global nCALCVAL; nCALCVAL=0;

% Set the default random number generator stream
RandStream.setGlobalStream(RandStream.create('mrg32k3a','Seed',58921));

%  restarting pool is time consuming  if (matlabpool('size') > 0);  matlabpool close; end; matlabpool open;   % default NumWorkers is #cores (Hyperthreading does not help matlab)
if (matlabpool('size') == 0);  matlabpool open;   end;  % default NumWorkers is #cores (Hyperthreading does not help matlab)

% Saved values to be used between repetitions of algorithm
% To enable tracing out equilibria which requires starting computation at
% previous equilibrium
ParRlg.x = 0;
ParRlg.exit = 0;
ParRlg.value = 0;
ParRlg.entry = 0;
ParRlg.v_entry = 0;
ParRlg.newvalue = 0;  
ParRlg.newx = 0;  
ParRlg.newexit = 0;  
ParRlg.isentry = 0;  

% Other Rlg params to be used
ParRlg.maxfirms_qencode        = 9;                                        % this value MUST BE SAME AS   #define maxfirms_qencode  in calcval_C.c
ParRlg.kmax                    = Params.kmax;
ParRlg.chk_foc_multi_eq        = 0;                                        % # of random restarts of within state foc search to attempt to find a different solution
ParRlg.foc                     = Params.rlg_foc;
ParRlg.no_force_exit           = 0;                                        % if 1, never bump firms off lowest rung of ladder.  Instead give free bump
ParRlg.w0_exit                 = 1;
ParRlg.damp                    = 0;                                        % weight on old values for dampening (might help if problems converging)
ParRlg.lb                      = zeros(1+2*Params.max_firms, 1);
ParRlg.ub                      = ones(1+2*Params.max_firms, 1); 
ParRlg.ub(1:Params.max_firms)  = ParRlg.ub(1:Params.max_firms)+99999999;


% Test Rlg values make sense
if ((Params.entry_high>1e99 && Params.scrap_val<-1e99) || (Params.entry_high<.0001 && Params.scrap_val<.0001))
    ParRlg.no_force_exit = 1;
end
if (ParRlg.no_force_exit)
    disp('hey: firms cannot be bumped off lowest rung since ParRlg.no_force_exit = 1 in eql_ma_foc.m');
end

if (Params.rlg_y>0 && Params.wstar>0)
    disp(' '); 
    disp('should not use both Params.rlg_y & Params.wstar');  
end
if (Params.rlg_y>0 && Params.rlg_y<Params.mc)
    disp(' '); 
    disp('must have Params.rlg_y>mc else profits = zero');
end

if (Params.rlg_outgood)
    disp('outside good is turned on');
else
    disp('outside good is turned off');
end

if (ParRlg.w0_exit == 1)
    disp('enforcing exit at w = 0 (as done in original PM code to signal room for entrants)');
else
    disp('not enforcing exit at w = 0 (as done in original PM code to signal room for entrants.  Okay to do this if no entry/exit.')
end

if (Params.rlg_scrap_hi == 0) 
    disp(' '); 
    disp('hey: must revise fixed scrap to address profits earned in exit period-->nval()for exitors depends on j in optimize! - Aborting!')
    Params.rlg_scrap_hi = 0.2;
%    return;
end

if (Params.rlg_scrap_hi-Params.scrap_val<.2)
    disp('hey: solving foc in each state may be difficult since range of scrap values is less than .2')
end

disp(['weight on old values is ParRlg.damp = ' num2str(ParRlg.damp)])
disp(['execution protocol-->' prefix 'markov.dp']);


disp(' ');  
disp('**** computing dynamic equilibrium using eql_ma_foc.m ****');   
disp(' ');


% Build binomial table
if (Params.wstar == 0)
    ParRlg.kmax = ParRlg.kmax+1;   % wstar = 0 case needs tables 1 beyond kmax @
end
binom = eye(Params.max_firms+ParRlg.kmax+1);
binom = [zeros(Params.max_firms+ParRlg.kmax+1, 1) binom];
for i = 2:Params.max_firms+ParRlg.kmax+1;
    binom(i, 2:i) = binom(i-1, 2:i)+binom(i-1, 1:i-1);
end
if (Params.wstar == 0)
    ParRlg.kmax = ParRlg.kmax-1;   % restore kmax
end

% global binomv binomcols;   binomv = reshape(binom', [], 1);   binomcols = size(binom', 1);

% RP - i keep the while here, though it only runs once.
nfirms = Params.max_firms;                                                  % hardcoded ignoring of option to use equilibrium from industry with fewer firms as starting point
while (nfirms <= Params.max_firms)

    % Start counting execution time
    d1 = clock;

    oneton = (1:nfirms)';

    ParRlg.wmax = binom(nfirms+1+ParRlg.kmax, ParRlg.kmax+2);
    disp(['number of firms: ' int2str(nfirms)  '  number of states: ' int2str(ParRlg.wmax) ]);
    disp('initializing...');

    filename = [prefix 'pr.' int2str(nfirms) 'f' '.mat'];
    load(filename);
    ParRlg.profit = profit;

    if (nfirms<ParRlg.maxfirms_qencode)
        filename = [prefix 'markov.dtable.' int2str(nfirms) 'f.' int2str(ParRlg.kmax) 'k'];
        ParRlg.two_n = 2^(nfirms-1);
        if exist([ filename,  '.mat'], 'file')
            ParRlg.multfac1 = (ParRlg.kmax+1).^(oneton-1);
            load([filename '.mat']);
            ParRlg.dtable  = dtable;
            ParRlg.etable1 = etable1;
            ParRlg.mask    = mask;
            disp('loaded dtable etable1 & mask instead of building them...')
        else
            if (nfirms>1)
                ParRlg.mask = mask_gen(nfirms, ParRlg.two_n);
                disp('mask done');
            end
            
            ParRlg.dtable = zeros(nfirms, ParRlg.wmax);
            for i = 1:ParRlg.wmax;
                ParRlg.dtable(:, i) = decode(binom,i,nfirms); 
            end
            
            ParRlg.multfac1 = (ParRlg.kmax+1).^(oneton-1);

            etable1 = zeros((ParRlg.kmax+1)^nfirms, 1);
	    kmax = ParRlg.kmax;
	    
            parfor i = 1:(kmax+1)^nfirms
                 msk = 0;
                 k = i-1;
                 if (mod(i, 100000) == 0)
                     disp(['Now at etable1 row ' int2str(i) ' of ' int2str((kmax+1)^nfirms)]);  
                 end
                 for t0 = 1:nfirms
                     j = (kmax+1)^t0;
                     msk = [((mod(k, j))*(kmax+1)/j); msk];
                     k = k-(mod(k, j));
                 end
                 d = (sortrows(msk(1:nfirms, 1), -1));
		 etable1(i) = encode2(binom,d,nfirms);
	    end
	    ParRlg.etable1 = etable1;
            disp('etable1 done');

            if (nfirms>1)
                disp('saving dtable etable1 & mask');
                filename = [prefix 'markov.dtable.' int2str(nfirms) 'f.' int2str(ParRlg.kmax) 'k'];
                dtable  = ParRlg.dtable;
                etable1 = ParRlg.etable1;
                mask    = ParRlg.mask;
                save([filename '.mat'],  'dtable', 'etable1', 'mask');
            end
        end
    else
        ParRlg.two_n = 2^(nfirms-1);
        ParRlg.etable1 = 0; 
	ParRlg.multfac1 = 0;
	ParRlg.dtable = 0; 
        ParRlg.mask = mask_gen(nfirms, ParRlg.two_n);
        disp('mask done.  Not building etable1 for qencode since nfirms>maxnfirms_qencode');
        ParRlg.dtable = zeros(nfirms, ParRlg.wmax);
        for i = 1:ParRlg.wmax
            ParRlg.dtable(:, i) = decode(binom,i,nfirms); 
        end
    end

    %initialize old values
    ParRlg.oldx = zeros(ParRlg.wmax, nfirms);
    ParRlg.oldexit = ones(ParRlg.wmax, nfirms);
    ParRlg.oldvalue = zeros(ParRlg.wmax, nfirms);
    ParRlg.oldentry = zeros(ParRlg.wmax, 1);
    ParRlg.oldv_entry = zeros(ParRlg.wmax, 1);
    
    ParRlg.foc_multi_eq = zeros(ParRlg.wmax, 1);
    
    % This if is never executed (the while look runs only once...)
    if (size(ParRlg.x, 1) == size(ParRlg.oldx, 1))
        disp('initializing with rlg_values from previous model...');
        ParRlg.oldx = ParRlg.x;
        ParRlg.oldexit = ParRlg.exit;
        ParRlg.oldvalue = ParRlg.value;
        ParRlg.oldentry = ParRlg.entry;
        ParRlg.oldv_entry = ParRlg.v_entry;
    end
    
    if (Params.delta == 0 && Params.wstar <= 0)
        do_w = (1:ParRlg.wmax)';
        do_w = do_w( ParRlg.dtable(1,:)'<ParRlg.kmax );
        ParRlg.oldx(do_w, :)       = ParRlg.oldx(do_w, :)-1;
        ParRlg.oldexit(do_w, :)    = ParRlg.oldexit(do_w, :)-1;
        ParRlg.oldvalue(do_w, :)   = ParRlg.oldvalue(do_w, :)-1;
        ParRlg.oldentry(do_w, :)   = ParRlg.oldentry(do_w, :)-1;
        do_w = (1:ParRlg.wmax)';  
        do_w = do_w( ParRlg.dtable(1,:)'==ParRlg.kmax );
        disp('only solving for values & policies over states with firm 1 at the frontier, since delta = 0');
    else
        do_w = (1:ParRlg.wmax)';
    end
    ParRlg.wmax2 = size(do_w, 1);

    ParRlg.newx = ParRlg.oldx;
    ParRlg.newexit = ParRlg.oldexit;
    ParRlg.newvalue = ParRlg.oldvalue;
    ParRlg.isentry = ParRlg.oldentry;

    %  OLD CODE
    % if Params.rlg_inv>0;
    %     if 0 && Par.batchjob<3600;
    %         aeff = Params.inv_mult*(1+Params.rlg_inv*rev((0:Params.ParRlg.kmax+1)));
    %         disp(['using linear spillover with coeff of ' num2str(Params.rlg_inv) ' , yielding aeff = ' num2str(aeff')]);
    %     else;
    %         aeff = [ones(Params.ParRlg.kmax, 1)*Params.inv_mult*(1+Params.rlg_inv); Params.inv_mult];
    %         disp(['using common spillover for all laggards, yielding aeff = ' num2str(aeff')]);
    %         if Params.delta>0 || Params.wstar>0; disp('hey: probably should not use this form of spillover(set in init.h), since leader often not at frontier');
    %         end;
    %     end;
    % else;
    %     aeff = Params.inv_mult*ones(Params.ParRlg.kmax+1, 1);
    % end;
    
    % construct aeffT table (wmax by nfirms) of investment efficiencies as function of state 
    if (Params.rlg_inv > 0)    % Params.rlg_inv > 0 --> single-step function giving common higher efficiency to all laggards
      ParRlg.aeffT = Params.inv_mult*(1+Params.rlg_inv*(repmat(max(ParRlg.dtable),nfirms,1) > ParRlg.dtable))';
    else                       % Params.rlg_inv <= 0 --> raise efficiency more for firms farther from frontier (not from kmax)
      ParRlg.aeffT = Params.inv_mult*(1-Params.rlg_inv*(repmat(max(ParRlg.dtable),nfirms,1)-ParRlg.dtable))';
    end
    
    ParRlg.minx1 = 0;
    if (Params.rlg_min_innov1>0)
      ParRlg.minx1 = Params.rlg_min_innov1 / ((1-Params.rlg_min_innov1)*Params.inv_mult);
    end
    
    % set entry state, entry_T (T for Table) based on Params.entry_at
    if (Params.entry_at > 0)
        ParRlg.entry_T = Params.entry_at + zeros(ParRlg.wmax,1);         % entry at fixed rung on quality ladder 
    else
        ParRlg.entry_T = max((ParRlg.dtable(1,:)'+Params.entry_at), 1);  % entry at higher of rung 1 or Params.entry_at below leader 
        ParRlg.entry_T = ParRlg.entry_T + (ParRlg.dtable(1,:)'==0)*(ParRlg.kmax+Params.entry_at-1);   % if empty industry,  enter at Params.entry_at below kmax
    end

    %%%
    % StartIterating:  label for goto statement in gauss was here
    %%%
    
    maxiter = 400/(1-ParRlg.damp);

    disp(['contraction...over ' int2str(ParRlg.wmax2) ' states of ' int2str(ParRlg.wmax) ' in dtable with rlg_foc = ' int2str(ParRlg.foc) ' min frontier firm innov = ' num2str(Params.rlg_min_innov1) ' maxiter = ' int2str(maxiter) ' tol = ' num2str(Par.tol) ' foc_tol = ' num2str(Par.foc_tol) '...']);
    if (ParRlg.chk_foc_multi_eq)
        disp(['in each iteration and at each state, restarting search for foc solution at ' int2str(ParRlg.chk_foc_multi_eq) ' random policies to investigate presence of multiple equilibria.']);
    end

    norm = Par.tol+99;
    irlg = 1;
    normx = norm;
    normexit = norm;
    normentry = norm;
    avgnorm = norm;
    prevnorm = norm+1;
    prevavgnorm = norm+1;
    t0 = clock;
    norm100 = norm+9;

    while (norm>Par.tol) && (irlg<maxiter || (norm-prevnorm<.01 && irlg<400) || (avgnorm-prevavgnorm<1e-7 && irlg<400))
        prevnorm = norm;
        prevavgnorm = avgnorm;

        if (mod(irlg, 100) == 0)
            norm100 = (norm+prevnorm)/2;
        else
            if (irlg>100 && irlg-100*floor(irlg/100)>31 && norm>norm100)
                disp('aborting contraction since current norm is worse than prev 100th iteration, suggesting convergence issues, ...');
                break;
            end
        end

        ParRlg = contract(Par,ParRlg,binom,Params,do_w,nfirms);
        
        abs_error = abs(ParRlg.oldvalue-ParRlg.newvalue);
        norm = max(max(abs_error));
        avgnorm = mean(mean(abs_error));
        normx = max(max(ParRlg.oldx-ParRlg.newx));
        normexit = max(max(ParRlg.oldexit-ParRlg.newexit));
        normentry = max(ParRlg.oldentry-ParRlg.isentry);

        irlg = irlg+1;
        if(mod(irlg, 10) == 0 || nfirms>9)
            disp(['iteration: ' int2str(irlg) ' avg sec ' sprintf('%8.3f', (etime(clock, t0)/(1+9*(nfirms <= 9)))) ' norm,avgnorm:' sprintf('%8.1e', [norm avgnorm]) '  norm x~exit~entry:' sprintf('%8.1e',normx) ' ~ ' sprintf('%8.1e',normexit) ' ~ ' sprintf('%8.1e',normentry) '  maxv ' num2str(max(ParRlg.newvalue(:, 1)))]);
            t0 = clock;
        end

%        if (nfirms<7 && norm>.001 && irlg>100)
%
%            [aaaaa normind1] = max(abs_error');
%            [aaaaa normind2] = max(abs_error);
%
%            if Par.batchjob == 0 && 0;  normind1 = qencode(ParRlg,binom,Params,[7; 7; 5],nfirms);  normind2 = 1;  end;
%            normcode = ParRlg.dtable(:, normind1)';
%
%            aaaaa = abs(ParRlg.newx-ParRlg.oldx);        [xind1 xind1] = max(max(aaaaa'));   xcode = ParRlg.dtable(:, xind1)';
%            bbbbb = abs(ParRlg.newexit-ParRlg.oldexit);  [eind1 eind1] = max(max(bbbbb'));   ecode = ParRlg.dtable(:, eind1)';
%            disp(['[norm avgnorm]:' num2str(padr(norm, 8, 5)) ' ' num2str(padr(avgnorm, 11, 8)) ''])
%            disp(['max norm:firm' num2str(padr(normind2, 1, 0)) 'at state' num2str(normcode) '; [old new] value:'])
%            disp([num2str(ParRlg.oldvalue(normind1, :)) '[ ]'])
%            disp([num2str(ParRlg.newvalue(normind1, :)) '[old new] x:'])
%            disp([num2str(ParRlg.oldx(normind1, :)) '[ ]'])
%            disp([num2str(ParRlg.newx(normind1, :)) '[old new] exit:'])
%            disp([num2str(ParRlg.oldexit(normind1, :)) '[ ]'])
%            disp([num2str(ParRlg.newexit(normind1, :)) 'entry:'])
%            disp(num2str([ParRlg.isentry(normind1) ParRlg.oldv_entry(normind1)]))
%            disp(['[old new x at' num2str(xcode) ' is ' num2str(ParRlg.oldx(xind1, :)) ' ]' num2str(ParRlg.newx(xind1, :))])
%            disp(['[old new exit at' num2str(ecode) ' is ' num2str(ParRlg.oldexit(eind1, :)) ' ]' num2str(ParRlg.newexit(eind1, :))])
%        end
        
        if (ParRlg.damp)
            i = max([ParRlg.damp; (.5*(norm>prevnorm))*rand(1,1)]);
            if (norm>prevnorm && ParRlg.damp == 0)
                disp(['norm>prevnorm, dampening weight on old = ' num2str(i) 'at iter = ' int2str(irlg)]);
            end
            ParRlg.oldx = ParRlg.oldx*i+(1-i)*ParRlg.newx;
            ParRlg.oldexit = ParRlg.oldexit*i+(1-i)*ParRlg.newexit;
            ParRlg.oldvalue = ParRlg.oldvalue*i+(1-i)*ParRlg.newvalue;
            ParRlg.oldentry = ParRlg.oldentry*i+(1-i)*ParRlg.isentry;
        else
            ParRlg.oldx = ParRlg.newx; 
            ParRlg.oldvalue = ParRlg.newvalue; 
            ParRlg.oldexit = ParRlg.newexit; 
            ParRlg.oldentry = ParRlg.isentry;
        end
    end
    
    disp(' ');
    disp([int2str(nfirms) ' firm(s) scenario completed;  execution time = ' num2str(etime(clock, d1)/6000) ' minutes, #iterations = ' int2str(irlg)  ]);  % ' #Calcval = ' int2str(nCALCVAL)
    disp(['tol = ' num2str(Par.tol) ' norm = ' num2str(norm) ' prevnorm~avgnorm = ' num2str([prevnorm prevavgnorm]) ' norm x~exit~entry = ' num2str([normx normexit normentry])]);

    % Display results of multiple equilibria test
    if (max(ParRlg.foc_multi_eq)>0)
        d2 = selif((1:ParRlg.wmax), ParRlg.foc_multi_eq>0);
        if (size(d2, 1)>30)
            d2 = sortrows(d2, -1);
            d2 = d2(1:30);
        end
        kk = 1.111+zeros(size(d2, 1), 1);
        disp('multiple foc solutions were obtained most frequently at the following states:#[multieq state entry investment exit]')
        disp( num2str([ParRlg.foc_multi_eq(d2) kk ParRlg.dtable(d2, :) kk ParRlg.isentry(d2) kk ParRlg.newx(d2, :) kk ParRlg.newexit(d2, :)]));
    end
    
    if (norm>Par.tol) 
        disp(' ');
%
% goto is not allowed in matlab. Disabling next two blocks which switched between FOC and non-foc when one method failed to converge
%        if (ParRlg.foc == 0 && Params.rlg_foc == 0 && (Par.batchjob == 0 || ParRlg.minx1 == 0))
%            ParRlg.foc = 1;
%            maxiter = 200; 
%            disp('non-foc failed to converge...retrying from current with foc method...')
%            goto ('startiterating');
%        end
%        if (0 && Par.batchjob>0 && ParRlg.minx1 == 0)
%            disp('non-foc & possibly foc failed to converge...retrying non-foc from current after setting Params.rlg_min_innov1 = .01...')
%            Params.rlg_min_innov1 = .01;
%            ParRlg.minx1 = Params.rlg_min_innov1/((1-Params.rlg_min_innov1)*Params.inv_mult);
%            ParRlg.foc = 0; 
%            maxiter = 200;
%            goto ('startiterating');
%        end

        if (ParRlg.foc>0) 
            disp('foc method failed to converge...');
        end
        if (ParRlg.foc == 0)
            disp('non-foc method failed to converge...');
        end
    end

    probw = ParRlg.aeffT.*ParRlg.newx;   
    probw = probw./(1+probw);

    disp(' ');     
    disp('v & policies at select states:[w v entry innovation exit]')

    w = ones(nfirms, 1);
    w_disp(ParRlg, binom, w, probw, nfirms);
    disp(' ');

    w = ParRlg.kmax*ones(nfirms, 1);
    w_disp(ParRlg, binom, w, probw, nfirms);
    for kk = nfirms:-1:2
        w(kk) = 0;
        w_disp(ParRlg, binom, w, probw, nfirms);
    end
    disp(' ');

    w = max(ParRlg.entry_T)*ones(nfirms, 1);   
    w(1) = ParRlg.kmax;
    w_disp(ParRlg, binom, w, probw, nfirms);
    for kk = nfirms:-1:2
        w(kk) = 0;
        w_disp(ParRlg, binom, w, probw, nfirms);
    end
    disp(' ');

    w = zeros(nfirms, 1); 
    w(1) = ParRlg.kmax;
    for kk = ParRlg.kmax:-1:1
        w(2) = kk;
        w_disp(ParRlg, binom, w, probw, nfirms);
    end
    disp(' ');

    w = ones(nfirms, 1);   
    w(1) = ParRlg.kmax;   
    w(nfirms) = 0;
    for kk = ParRlg.kmax:-1:1
        w(2) = kk;
        w_disp(ParRlg, binom, w, probw, nfirms);
    end
    disp(' ');

    w = ones(nfirms, 1);  
    w(1) = ParRlg.kmax;
    for kk = ParRlg.kmax:-1:1;
        w(2) = kk;
        w_disp(ParRlg, binom, w, probw, nfirms);
    end
    disp(' ');

    if (nfirms>2)
        w = zeros(nfirms, 1);   
        w(1) = ParRlg.kmax;   
        w(2) = ParRlg.kmax;
        for kk = ParRlg.kmax:-1:1
            w(3) = kk;
            w_disp(ParRlg, binom, w, probw, nfirms);
        end
    end
    disp(' ');

    if (nfirms>2)
        w = zeros(nfirms, 1);   
        w(1) = ParRlg.kmax;   
        w(2) = ParRlg.kmax-1;
        for kk = ParRlg.kmax-1:-1:1;
            w(3) = kk;
            w_disp(ParRlg, binom, w, probw, nfirms);
        end
    end
    disp(' ');

    w = (nfirms:-1:1)';   
    w = min([w'; ParRlg.kmax*ones(1, nfirms)])';
    w(1) = ParRlg.kmax;
    w_disp(ParRlg, binom, w, probw, nfirms);
    disp(' ');

    w = (sum(ParRlg.newexit'>0) == nfirms);
    if (sum(w)>1)
        irlg = selif((1:ParRlg.wmax)', w);
        disp('all firms exit at the following states: [w v x entry exit]');
        i = -9.99*ones(size(irlg, 1), 1);
        disp(num2str([ ParRlg.dtable(:, irlg)' i ParRlg.newvalue(irlg, :) i ParRlg.newx(irlg, :) i ParRlg.isentry(irlg) i ParRlg.newexit(irlg, :) ]));
    end
    disp(' ');

    if (Params.wstar>0)
        w = zeros(nfirms, 1);  
        w(1) = ParRlg.kmax-1;
        kk = ParRlg.newx(qencode(ParRlg,binom,w,nfirms):(qencode(ParRlg,binom,[ParRlg.kmax; zeros(nfirms-1, 1)],nfirms)-1), 1);
        irlg = max(kk);
        if (irlg>0)
            disp(['warning: positive investment of ' num2str(irlg) ' recorded at 2nd highest efficiency level']);
            disp('please consider increasing the maximum efficiency level (kmax) since using PM wstar to bound statespace.');
        end
    else
        w = zeros(nfirms, 1); 
        w(1) = ParRlg.kmax;
        kk = ParRlg.newx(qencode(ParRlg,binom,w,nfirms):ParRlg.wmax, 1);
        irlg = max(kk);
        disp(['max investment of ' num2str(irlg) ' recorded at highest efficiency level (which is okay since not using PM wstar to bound statespace)']);
        irlg = (kk-ParRlg.minx1)<1e-15;
        % NEED TO FIX THIS CODE  i in next line goes to size(kk) but irlg
        % is length qencode(w):wmax ...
        % if (sum(irlg)>0)
        %     i = qencode(ParRlg,binom,w,nfirms):size(kk, 1);
        %     irlg = i(irlg);
        %     disp(['minimum investment of ' num2str(ParRlg.minx1) ' by frontier firm at ' int2str(size(irlg, 1)) ' of ' int2str(size(kk, 1)) ' states with leader at frontier:[w v x entry exit]'])
        %     i = -9.99*ones(size(irlg, 1), 1);
        %     disp(num2str([ParRlg.dtable(:, irlg)' i ParRlg.newvalue(irlg, :) i ParRlg.newx(irlg, :) i ParRlg.isentry(irlg) i ParRlg.newexit(irlg, :)] ));
        % end
    end
    disp(' ');  
    
    w = sum((ParRlg.dtable(:, do_w) == 1 & ParRlg.newexit(do_w, :)'<1-1e-8));
    if (sum(w)>0)
        kk = do_w(w>0);
        disp(['firms on lowest rung w/pr(exit)<1 at ' int2str(size(kk, 1)) ' states.  They get free bump if ParRlg.w0_exit = 0,  else exit w/avg scrap when fall off ladder.  Hardcoded ParRlg.w0_exit is ' int2str(ParRlg.w0_exit)]);
        disp('first few of such states:[value invest entry exit], separated by columns of -9.99');
        kk = kk(1:min([10; size(kk, 1)]));  
        w = -9.99*ones(size(kk, 1), 1);
        disp(num2str([ParRlg.dtable(:, kk)' w ParRlg.newvalue(kk, :) w ParRlg.newx(kk, :) w ParRlg.isentry(kk) w ParRlg.newexit(kk, :) ]));
    end
    disp(' ');

    w = sum((ParRlg.dtable(:, do_w) == 1 & ParRlg.newexit(do_w, :)'<1-1e-8 & ParRlg.newx(do_w, :)'>1e-8));
    if (sum(w)>0)
        kk = do_w(w>0);
        disp(['firms on lowest rung w/invest>0 at ' int2str(size(kk, 1)) ' states.  They get free bump if ParRlg.w0_exit = 0,  else exit w/avg scrap when fall off ladder.  Hardcoded ParRlg.w0_exit is ' int2str(ParRlg.w0_exit)]);
        disp('first few of such states:[value invest entry exit], separated by columns of -9.99');
        kk = kk(1:min([10; size(kk, 1)]));  
        w = -9.99*ones(size(kk, 1), 1);
        disp(num2str([ParRlg.dtable(:, kk)' w ParRlg.newvalue(kk, :) w ParRlg.newx(kk, :) w ParRlg.isentry(kk) w ParRlg.newexit(kk, :) ]));
    end
    disp(' ');

    % Save results

    filenam2 = [prefix 'markov.' int2str(nfirms) 'ot'];

    %     w = zeros(length(ParRlg.isentry), 1)-1;
    %     
    %     i = [ParRlg.dtable' w ParRlg.newvalue];                       
    %     save([ filenam2 '.w_v' '.mat'],  'i');
    %     
    %     i = [ParRlg.dtable' w ParRlg.newx    ];                       
    %     save([ filenam2 '.w_x' '.mat'],  'i');
    %     
    %     i = [ParRlg.dtable' w probw w ParRlg.newx w ParRlg.aeffT];   
    %     save([ filenam2 '.w_innov_x_a' '.mat'],  'i');
    %     
    %     i = [ParRlg.dtable' w ParRlg.newexit ];                       
    %     save([ filenam2 '.w_exit' '.mat'],  'i');
    %     
    %     i = [ParRlg.dtable' w ParRlg.isentry w ParRlg.oldv_entry];   
    %     save([ filenam2 '.w_entry_v' '.mat'],  'i');
    %
    %     i = [ParRlg.dtable' w ParRlg.isentry w ParRlg.oldv_entry];   
    %     i = i(ParRlg.dtable(nfirms, :)' == 0, :);     
    %     save([ filenam2 '.w_entry_v_room' '.mat'],  'i');

    disp(['generating data for comparative statics in -->' filenam2]);
    newvalue    = ParRlg.newvalue;
    newx        = ParRlg.newx;
    newexit     = ParRlg.newexit;
    isentry     = ParRlg.isentry;
    oldv_entry  = ParRlg.oldv_entry;
    save([filenam2 '.mat'],  'newvalue', 'newx', 'probw', 'newexit', 'isentry', 'oldv_entry');
    % newvalue is the value function
    % newx is the investment policy function
    % probw is the probability of innovation (based on newx)
    % newexit is the exit policy function (i.e., probability of exit)
    % isentry is the entry policy function (i.e., probability of entry)
    % oldv_entry is the potential entrant's value function

    if (Params.delta == 0 && Params.wstar <= 0)
        irlg = (sum(ParRlg.newx,2)-ParRlg.minx1*sum(ParRlg.dtable == ParRlg.kmax)')<1e-15 & ParRlg.isentry == 0 & sum(ParRlg.newexit,2) == 0;
        if (sum(irlg)>0)
            disp(' ');   
            disp('absorbing states with min investment,  no entry,  no exit:');  
            disp( int2str( ParRlg.dtable(:, irlg)' ));   
            disp(' ');
        end
    end

    nfirms = nfirms+1;
end

% Save values - totally redundant, as we do nothing with them
if (norm <= Par.tol)
    ParRlg.x = ParRlg.oldx;
    ParRlg.exit = ParRlg.oldexit;
    ParRlg.value = ParRlg.oldvalue;
    ParRlg.entry = ParRlg.oldentry;
    ParRlg.v_entry = ParRlg.oldv_entry;
end

Params.eql_done = (norm <= Par.tol);
save(Par.configfile, 'Params');
end % End of function eql_ma_foc



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function contract - This procedure does one iterative step on investment 
% and the value fn.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rlgret = contract(Par,ParRlg,binom,Params,do_w,nfirms)
% First: check for which values of w_s would a firm want to enter 
%        Entry decision is made at very beginning of each period  
%        BEFORE exit is implemented (done in optimize),           
%        BUT I think isentry should account for expected exit.    
%        Is this an error in this PM code (C also) ?              
%        chkentry() calls calcval() which integrates over invest  
%        outcomes, but does not implement exit, which is done in  
%        optimize() prior to calling calcval() when getting newx  
%        Of course, when I modify code to allow for random scrap  
%        calcval() will integrate over exit as well as invest.    
  
tic;
% RP - changed wmax to wmax2, as that is the correct size of do_w in eql_ma_foc
foc = zeros(ParRlg.wmax, 2*nfirms+1);
rflag = zeros(ParRlg.wmax, 1);

% Create local matrices for parfor
tmp_foc_multi_eq    = ParRlg.foc_multi_eq;
tmp_oldv_entry      = ParRlg.oldv_entry;
tmp_isentry         = ParRlg.isentry;
tmp_newx            = ParRlg.newx;
tmp_newvalue        = ParRlg.newvalue;
tmp_newexit         = ParRlg.newexit;

parfor w1 = do_w' ;
    [ret1,ret2,ret3,ret4] = contract_w(Par,ParRlg,binom,Params,w1,nfirms);
    foc(w1, :) = ret1;
    rflag(w1, :) = ret2;
    tmp_foc_multi_eq(w1,:)  = ret3;
    tmp_oldv_entry(w1,:)    = ret4.oldv_entry(w1,:);
    tmp_isentry(w1,:)       = ret4.isentry(w1,:);
    tmp_newx(w1,:)          = ret4.newx(w1,:);
    tmp_newvalue(w1,:)      = ret4.newvalue(w1,:);
    tmp_newexit(w1,:)       = ret4.newexit(w1,:);
end

% Restore local matrices into ParRlg
ParRlg.foc_multi_eq    = tmp_foc_multi_eq;
ParRlg.oldv_entry      = tmp_oldv_entry;
ParRlg.isentry         = tmp_isentry;
ParRlg.newx            = tmp_newx;
ParRlg.newvalue        = tmp_newvalue;
ParRlg.newexit         = tmp_newexit;

if (ParRlg.foc>0)
    rmax = max(rflag(do_w));
    if (rmax>1e-5)
        disp(['mean ~ max log10(foc) = ' sprintf('%6.2f', log10(mean(rflag(do_w)))) ' ' sprintf('%6.2f', log10(rmax))  '  sec= ' num2str(toc) ]);
    end
end
rlgret = ParRlg;
end % End of function contract



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function contract_w - This procedure created so that one line of code 
% needs to be called for each state in contract(), which facilitates clean 
% multithreading.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ret1, ret2, ret3, rlgret] = contract_w(Par, ParRlg,binom, Params, w1, nfirms)

x_exit_entry = [ParRlg.oldx ParRlg.oldexit ParRlg.oldentry];
locw = ParRlg.dtable(:, w1);
foc_w1 = zeros(1, 2*nfirms+1);
rflag_w1 = 0;
foc_multi_eq_w1 = 0;
oneto2nfirms1 = (1:2*nfirms+1)';


if (ParRlg.foc <= 0)
    [tmp,ParRlg] = get_foc(Par,ParRlg,binom,Params,x_exit_entry(w1, :)',w1,oneto2nfirms1,nfirms);
    foc_w1 = tmp';
else
    if (locw(1) == 0)
        dosymm = 2*nfirms+1;
    else
        dosymm = (1:nfirms)';
        dosymm = dosymm( locw ~= ([-1; locw(1:nfirms-1)]) & locw>0 );
        dosymm = [dosymm;  (nfirms+dosymm) ];
        if (locw(nfirms) == 0)
            dosymm = [dosymm;  2*nfirms+1];
        end
    end

    do1 = dosymm( x_exit_entry(w1, dosymm)' ~= 0 & x_exit_entry(w1, dosymm)' ~= 1 );   % most 0, 1 policies persist,  so first do newton() holding them fixed.  Will check their optimality later

    if (size(do1, 1)<1)
        if (locw(1) == 0)
            do1 = 2*nfirms+1;   
        else
            do1 = 1;  
        end
    end

    if (size(do1, 1)>0)
        [tmp,ParRlg] = newton3(Par, ParRlg,binom, Params, x_exit_entry(w1, do1)', 'get_foc', w1, do1, nfirms);
        x_exit_entry(w1, do1) = tmp';
    end
    [tmp,ParRlg] = get_foc(Par,ParRlg,binom,Params,x_exit_entry(w1, :)',w1,oneto2nfirms1,nfirms);
    foc_w1 = tmp';
    rflag_w1 = max(abs(foc_w1'));

    if (size(do1, 1)<size(dosymm, 1) && rflag_w1>Par.foc_tol)
        do1 = dosymm;
        [tmp,ParRlg] = newton3(Par,ParRlg,binom,Params,x_exit_entry(w1, do1)','get_foc',w1,do1,nfirms);
        x_exit_entry(w1, do1) = tmp';
        [tmp,ParRlg] = get_foc(Par,ParRlg,binom,Params,x_exit_entry(w1, :)',w1,oneto2nfirms1,nfirms);
        foc_w1 = tmp';
        rflag_w1 = max(abs(foc_w1'));
    end

    if (rflag_w1>1e-2)
        disp(['foc: s ' int2str(w1) ' w1 ' int2str(ParRlg.dtable(:, w1)')  ' do1 ' int2str(do1') ' newx ' num2str(ParRlg.newx(w1, :),  ' %6.3f') '  newexit ' num2str(ParRlg.newexit(w1, :),  ' %5.3f') ' isentry ' num2str(ParRlg.isentry(w1),  ' %5.3f') ' log max abs foc ' num2str(log10(max(abs(foc_w1')))) ]);
        keyboard;
    end

    if (ParRlg.chk_foc_multi_eq>0)
        do1 = dosymm; 
        x_exit_entry_val = [x_exit_entry(w1, :) ParRlg.newvalue(w1, :)];
        for t0 = 0:ParRlg.chk_foc_multi_eq
            [tmp,ParRlg] = newton3(Par, ParRlg,binom, Params, ParRlg.lb(do1)+rand(size(do1, 1), 1).*(min([ParRlg.ub(do1)'; 10*ones(1, size(do1, 1))])-ParRlg.lb(do1)), 'get_foc', w1, do1, nfirms);
            x_exit_entry(w1, do1) = tmp';
            [tmp,ParRlg] = get_foc(Par,ParRlg,binom,Params,x_exit_entry(w1, :)',w1,oneto2nfirms1,nfirms);
            foc_w = tmp';
            if (max(abs(foc_w))<Par.foc_tol && max(abs(x_exit_entry(w1, do1)-x_exit_entry_val(do1)))>.001)
                foc_multi_eq_w1 = 1+foc_multi_eq_w1;
                break;
            end
        end
        ParRlg.newvalue(w1, :) = x_exit_entry_val(2*nfirms+2:3*nfirms+1);
        ParRlg.newx(w1, :) = x_exit_entry_val(1:nfirms);
        ParRlg.newexit(w1, :) = x_exit_entry_val(1+nfirms:2*nfirms);
        ParRlg.isentry(w1, :) = x_exit_entry_val(2*nfirms+1);
    end
end

ret1 = foc_w1;
ret2 = rflag_w1;
ret3 = foc_multi_eq_w1;
rlgret = ParRlg;
end % End of function contract_w



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function newton - Newton-Raphson search to find the root of the function
% objfunk. LB <= p <= UB is enforced.
% Note:  objfunk is a function name and must be passed as a string
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p, rlgret] = newton3(Par, ParRlg,binom, Params, p, objfunk, w, doo, nfirms)

lbw = ParRlg.lb; 
lbw(1:nfirms) = ParRlg.minx1*(ParRlg.dtable(:, w) == ParRlg.kmax);
np = size(p, 1); 
minstep = 1;
trymax = 200*np; 
deriv = zeros(np, np);

iter = 0;
while iter<(trymax*20)
    iter = iter + 1;
    
    [x,ParRlg] = feval(objfunk, Par, ParRlg, binom, Params, p, w, doo, nfirms);
    if (max(abs(x))<Par.foc_tol)
        break; 
    end
    
    % Get gradient
    for i = 1:np
        dp = p;
        if (ParRlg.ub(doo(i))-p(i)<0.0001)
            dp(i) = p(i)-.0001;
        else
            dp(i) = p(i)+.0001;
        end
        [dx,ParRlg] = feval(objfunk, Par, ParRlg, binom, Params, dp, w, doo, nfirms);
        deriv(:, i) = (dx-x)/(dp(i)-p(i));
    end

    % If gradient broken, restart from a different position
    if (rcond(deriv)<1e-14 || mod(iter, trymax) == 0)
        if mod(iter, trymax) == 0
            disp(['restarting newton() at rand() since not converging at iter ' int2str(iter) ' with w = ' int2str(ParRlg.dtable(:, w)') '  doo = ' int2str(doo')]);
            disp(['p = ' num2str(p') '  foc = ' num2str(x')]);
        else
            disp(['restarting newton() at rand() since singular deriv at iter ' int2str(iter) ' with w = ' int2str(ParRlg.dtable(:, w)') '  doo = ' int2str(doo')]);
            disp(['p = ' num2str(p') '  foc = ' num2str(x') ]);
            iter = trymax*(1+floor(iter/trymax));
	end
        pnew = lbw(doo)+rand(np, 1).*(min([ParRlg.ub(doo)'; 10*ones(1, np)])'-lbw(doo));
        disp(['pnew = ' num2str(pnew')]);
        minstep = max([.3;  minstep-.1]);
    else
      % pnew = inv(deriv);  
      pnew = deriv\x ;  % instead of inv(deriv)*x
        if (mod(iter, trymax)<5)
           % pnew = p-pnew*((.1+(mod(iter, trymax))/10)*rand(1,1))*x;   % start slow
            pnew = p-pnew*((.1+(mod(iter, trymax))/10)*rand(1,1));   % start slow
        else
            % pnew = p-pnew*(minstep+(1-minstep)*rand(1, 1))*x;
            pnew = p-pnew*(minstep+(1-minstep)*rand(1, 1));
        end
    end
    
    % Do newton step
    pnew = min([ParRlg.ub(doo)' ;  pnew' ])';
    pnew = max([lbw(doo)' ;  pnew' ])';
    p = pnew;
end

if (max(abs(x))>Par.foc_tol)
    disp(['newton failed with w = ' int2str(ParRlg.dtable(:, w)') '  doo = ' int2str(doo') '  returning p = ' num2str(p') '  foc = ' num2str(x')]);
end
rlgret = ParRlg;
end % End of function newton



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function get_foc - obtain foc for state w at candidate policy x_exit_ent
% x_exit_ent   are policies being changed to find root
% doo          are indices of policies being changed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ret0,rlgret] = get_foc(Par, ParRlg,binom, Params, x_exit_ent, w, doo, nfirms)

x_exit_entry_w = ([ParRlg.oldx(w, :) ParRlg.oldexit(w, :) ParRlg.oldentry(w)])';
if (size(doo, 1)>0)
    x_exit_entry_w(doo) = x_exit_ent; 
end

locw = ParRlg.dtable(:, w)';
for i = 2:nfirms
    if locw(i) == locw(i-1);
        x_exit_entry_w(i) = x_exit_entry_w(i-1);
        x_exit_entry_w(i+nfirms) = x_exit_entry_w(i-1+nfirms);
    end
end

ParRlg = chkentry_w(ParRlg, binom, Params, x_exit_entry_w, w, nfirms);
[ParRlg.newx(w, :), ParRlg.newvalue(w, :), ParRlg.newexit(w, :)] = optimize_w(Par,ParRlg,binom,Params,x_exit_entry_w,w,nfirms);

foc_w = x_exit_entry_w-([ParRlg.newx(w, :) ParRlg.newexit(w, :) ParRlg.isentry(w)])';

if (size(doo, 1) == 0)
    ret0 = 1e-99;
else
    ret0 = foc_w(doo);
end
rlgret = ParRlg;
end % End of function get_foc



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function optimize_w - This function calculates optimal investment, and 
% value fn., for a given industry structure w. Thus, a vector nfirms long 
% of each is returned.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function   [ret0, ret1, ret2] = optimize_w(Par, ParRlg,binom, Params, x_exit_entry, w, nfirms)

locw = ParRlg.dtable(:, w);
locwx = locw;
oval = ParRlg.oldvalue(w, :)';
ox = x_exit_entry(1:nfirms);
oexit = x_exit_entry(nfirms+1:2*nfirms);
nval = zeros(nfirms, 1);
nx = zeros(nfirms, 1);
nexit = zeros(nfirms, 1);
aeffw = ParRlg.aeffT(w,:)';
  
% This bit of code should not run
if (Params.rlg_scrap_hi == 0)
    disp('Hey: probably should not run fixed scrap code anymore - buggy');
    i = (min(oval) == Params.scrap_val)*(min(oval)-1)+(min(oval)>Params.scrap_val)*nfirms;
    if (i<nfirms) 
        locwx(i+1:nfirms) = zeros(nfirms-i, 1); 
    end
end

entered = x_exit_entry(2*nfirms+1);
locwe = locwx;
locwk = locwx;
locwe(nfirms) = ParRlg.entry_T(w);
locwk(nfirms) = ParRlg.kmax;
eexit = oexit;
eexit(nfirms) = 0;

doleap = (Params.rlg_leap>0 && nfirms>1); 
if doleap && nfirms>2
    doleap = (locw(nfirms-1) == 0); 
end

for j = 1:nfirms
    if (ParRlg.w0_exit && locw(j) == 0)
        v1 = Params.scrap_val;
        if (Params.rlg_scrap_hi ~= 0)
            v1 = (Params.rlg_scrap_hi+Params.scrap_val)/2; 
        end
        nval(j:nfirms) = v1*ones(nfirms-j+1, 1);
        nx(j:nfirms) = zeros(nfirms-j+1, 1);
        nexit(j:nfirms) = ones(nfirms-j+1, 1);
        break;
    end;

    v1 = 0; 
    v2 = 0;
    if (entered<1)
        [v1, v2] = calcval(ParRlg, binom, Params, j, locwx, ox, oexit, locw(j), aeffw, nfirms);
    end

    if (entered>0)
        [tempv1, tempv2] = calcval(ParRlg, binom, Params, j, locwe, ox, eexit, locw(j), aeffw, nfirms);
        if (doleap)
            [tempk1, tempk2] = calcval(ParRlg, binom, Params, j, locwk, ox, eexit, locw(j), aeffw, nfirms);
            tempv1 = (1-Params.rlg_leap)*tempv1+Params.rlg_leap*tempk1;
            tempv2 = (1-Params.rlg_leap)*tempv2+Params.rlg_leap*tempk2;
        end
        v1 = entered*tempv1+(1-entered)*v1;
        v2 = entered*tempv2+(1-entered)*v2;
    end

    if (v1 <= v2)
        r = 1.0;
    else
        r = 1.0/(Params.beta*aeffw(j)*(v1-v2));
    end
    
    r = min([max([r; 1E-10]); 1]);
    p = 1.0-sqrt(r);
    nx(j) = p/(aeffw(j)-aeffw(j)*p);

    % %% uncomment to play around with imposing minimum investment here
    % if (Par.batchjob == 0 && locw(j) == ParRlg.kmax && locw(min([nfirms; 2]))>0) %% or whatever conditions you want to trigger the min investment requirement
    %     nx(j) = max([ParRlg.minx1; nx(j)]); 
    % end

    nval(j) = Params.beta*(v1*p+v2*(1-p))-nx(j);

    if (Params.rlg_scrap_hi == 0)
        if (nval(j) <= Params.scrap_val)
            nexit(j:nfirms) = ones(nfirms-j+1, 1);
            nval(j:nfirms) = ones(nfirms-j+1, 1)*Params.scrap_val;
            nx(j) = 0;
            break;
        end
    else
        if (nval(j)<Params.rlg_scrap_hi)
            if (nval(j) <= Params.scrap_val)
                nexit(j) = 1;
                nval(j) = (Params.scrap_val+Params.rlg_scrap_hi)/2;
                nx(j) = 0;
            else
                nexit(j) = (Params.rlg_scrap_hi-nval(j))/(Params.rlg_scrap_hi-Params.scrap_val);
                nval(j) = (1-nexit(j))*nval(j)+nexit(j)*(Params.rlg_scrap_hi+nval(j))/2;
            end
        end
    end

    nval(j) = nval(j)+ParRlg.profit(w, j);

    if (Params.rlg_scrap_hi == 0)
        locwx(j) = (nval(j)>Params.scrap_val)*locw(j);
        locwe(j) = locwx(j);
    end
end

ret0 = nx';  
ret1 = nval';  
ret2 = nexit';
end % End of function optimize_w



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function chkentry_w - calculates for which value of other people's 
% omegas, would a firm want to enter, given that the market has room for 
% another firm.
%
% With FIXED scrap (i.e., orig PM), exit is implemented before looking up
% the isentry[] element, so do NOT implement exit here when getting isentry
% With RANDOM scrap, exit is integrated over in calcval()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rlgret] = chkentry_w(ParRlg, binom, Params, x_exit_entry, w, nfirms)
locw = ParRlg.dtable(:, w);
if (locw(nfirms) == 0)
    [~, v1] = calcval(ParRlg, binom, Params, nfirms, locw, x_exit_entry(1:nfirms), x_exit_entry(nfirms+1:2*nfirms), ParRlg.entry_T(w), ParRlg.aeffT(w,:)', nfirms);
    val = Params.beta*v1;

    doleap = (Params.rlg_leap>0 && nfirms>1); 
    if doleap && nfirms>2
        doleap = (locw(nfirms-1) == 0); 
    end
    
    if doleap
        [~, v1] = calcval(ParRlg, binom, Params, nfirms, locw, x_exit_entry(1:nfirms)', x_exit_entry(nfirms+1:2*nfirms)', max(locw(1),1), ParRlg.aeffT(w,:)', nfirms);
        val = (1-Params.rlg_leap)*val+Params.rlg_leap*Params.beta*v1;
    end
    ParRlg.oldv_entry(w) = val;
    ParRlg.isentry(w) = max([0; min([1; ((val-Params.entry_low)/(Params.entry_high-Params.entry_low))])]);
end
rlgret = ParRlg;
end % End of function chkentry_w



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function calcval - wrapper function for calcval_C
% This function calculates val = EEEV(.,.,.,.)p(.)p(.)p(.), where E
% represents sums, and this is the calculation of the 4-firm problem
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function   [vala, valb] = calcval(ParRlg, binom, Params, place, w, x, isex, k, aeffw, nfirms)
% Prepare location mask
if nfirms>1
    if (place == 1)
      locmask = [zeros(1, ParRlg.two_n); ParRlg.mask];
    elseif (place == nfirms)
      locmask = [ParRlg.mask; zeros(1, ParRlg.two_n)];
    else
      locmask = [ParRlg.mask(1:place-1, :); zeros(1, ParRlg.two_n); ParRlg.mask(place:nfirms-1, :)];
    end
else
    locmask = zeros(1, 1);
end
x(place) = 0;
w(place) = k;
isex(place) = 0;
p_up = aeffw.*x;   
p_up = p_up./(1+p_up);

% Prepare variables for C code
binomv = binom';  
binomv = binomv(:);  %  faster than binomv = reshape(binom', [], 1);
vala = 0;
valb = 0;
probmask_chksum = 0;

% Increase calcval counter and call the C function
% global nCALCVAL; nCALCVAL = nCALCVAL+1;
calcval_C(vala, valb, probmask_chksum, place, w, x, isex, locmask', ParRlg.kmax, nfirms, p_up, ParRlg.two_n, Params.rlg_scrap_hi, ParRlg.no_force_exit, Params.wstar, ParRlg.oldvalue', ParRlg.etable1, ParRlg.multfac1, binomv, Params.delta, Params.rlg_outgood);

% Test that returned probability sums to 1
if abs(1-probmask_chksum)>1e-10
    disp(['c code HEY: probmask_chksum not 1: 1-chksum= ' num2str(1-probmask_chksum) ' ; place k ' int2str([place k])  ' ; w ' int2str(w') ' ; isex ' num2str(isex') ' ; p_up ' num2str(p_up') ' ; vala ~ valb = ' num2str([vala valb]) ' ; locmask next block:']);
    disp(int2str(locmask));
end

end % End of function calcval



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A simple display function for results display
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function w_disp(ParRlg, binom, w, probw, nfirms)
i = qencode(ParRlg,binom,w,nfirms);
disp([ 'w:' sprintf(' %2.0f',w) '  v:' sprintf(' %10.4f', ParRlg.newvalue(i, :)) '  entry:'  sprintf(' %6.4f', ParRlg.isentry(i)) '  innov:' sprintf(' %6.4f', probw(i, :)) '  exit:' sprintf(' %6.4f', ParRlg.newexit(i, :)) ]);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decode/Encode functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ntuple = decode(binom, code, nfirms)
code = code-1;
ntuple = zeros(nfirms,1);
for i=1:nfirms
    digit = 0;
    while binom(digit+nfirms-i+2,digit+2) <= code
        digit = digit+1;
    end
    ntuple(i) = digit;
    code = code-binom(digit+nfirms-i+1,digit+1);
end
end % End of function decode


function ret2 = qencode(ParRlg, binom, ntuple, nfirms)
if (nfirms<ParRlg.maxfirms_qencode)
  ret2 = ParRlg.etable1(sum(ntuple.*ParRlg.multfac1)+1);
else
  ret2 = encode2(binom, ntuple, nfirms);
end
end % End of function qencode


function ret1 = encode2(binom, ntuple, nfirms)
binomcols = size(binom, 2);  
binomv = binom';
binomv = binomv(:);  % reshape(binom', [], 1);
% nfirms_oneton = (nfirms-(1:nfirms)');
nfirms_oneton = (nfirms-1:-1:0)';
ret1 = 1+sum(binomv((ntuple+nfirms_oneton)*binomcols+ntuple+1));
end % End of function encode2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function acronym - get textual description of types
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s,  etype,  itype] = acronym(Par, et, it)
if et == Par.uu_competition
    s = '';
    etype = 'markov perfect nash equilibrium';
elseif et == Par.uu_monopoly
    s = 'm';
    etype = 'monopoly';
elseif et == Par.uu_planner
    s = 's';
    etype = 'social planner';
end

if it == Par.uu_quality
    s = [s 'b'];
    itype = 'quality (differentiated products)';
elseif it == Par.uu_cost
    s = [s 'c'];
    itype = 'marginal cost (homogenous products)';
elseif it == Par.uu_capacity
    s = [s 'p'];
    itype = 'capacity (homogenous products)';
end

if Par.batchjob > 0 
    s = ['bj' num2str(Par.batchjob) '_' s]; 
end
end % End of function acronym



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function mask_gen - generate a mask
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mask = mask_gen(nfirms, two_n)
mask = zeros(nfirms-1, two_n);
for i = 1:two_n
    msk = 0;
    j = 2; 
    k = i-1;
    while j <= two_n
        if mod(k, j) == 0
            msk = [0;  msk];
        else
            k = k-j/2;  
            msk = [1;  msk];
        end
        j = j*2;
    end
    mask(:, i) = msk(1:nfirms-1);
end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implement gauss selif 
function y = selif(x, e)
y = x(e == 1, :);
end
