function [q,dqdparam,elqp] = gtype(DIR,JOB,imap,x0,p,tfinal,ilist,klist,iplot,solution_only,NextGen,R0_only)
%
%  Inputs
%    x0 = initial conditions  
%    p  = parameters
%    psi = adjoint data
%    tfinal = final time
%    iplot  = plot level
%    ilist  = variables for which sensitivity information is required
%    klist  = parameters for which sensitivity information is required
%    iplot  = structure containing plotting information
%    solution_only = Option for computing only the solutions without
%                   sensitivities and elasticities, for faster computation
%    
% Variables
%    t runs from 1 to tfinal and is the time step
%    i runs from 1 to xdim and is the row index
%    j runs from 1 to xdim and is the column index
%    k runs from 1 to kdim and is the parameter index
%
%
% Treats the system as a general nonlinear iteration
%       x(t+1,theta) = g(x(t,theta),theta)
%

global NCOLUMNS
NCOLUMNS=[1  2*ones(1,3)  3*ones(1,2)  4*ones(1,6)  5*ones(1,8)  6*ones(1,10)  7*ones(1,12)  8*ones(1,14)  9*ones(1,16)];

xdim=length(x0);
kdim=length(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                  %
%     Compute solutions (and sensitivities)                                                       %
%                                                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if imap == 1  
    [nt,t,x,dxdp] = solve_map(x0,p,tfinal,xdim,kdim,solution_only);
else 
    [nt,t,x,dxdp] = solve_ode(x0,p,tfinal,xdim,kdim,solution_only);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                  %
%     Print solutions                                                             %
%                                                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Prints the solution at tfinal (Always printed) 
fprintf('\n');
for i=1:xdim
    fprintf('x(%3i) at tfinal = %13.6e \n', i, x(i,nt))      
end
fprintf('\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                  %
%     Plot solutions                                                                               %
%                                                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

output_directory=[DIR,'/',JOB];
[s,mess,messid]=mkdir(output_directory);

% Plot solutions -- always done by default 
plot_solutions(output_directory,imap,t,x,ilist)

if(solution_only == 0)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     Print and plot sensitivities                                                                 %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Prints sensitivity of solution if i_plot_sensitivity = 1
    if (iplot.sensitivities == 1 && (iplot.dxdp.true == 1 || iplot.dxdp.var == 1 || iplot.dxdp.param == 1))
        for i=1:xdim
            for k=1:kdim
                fprintf('Sensitivity of x(%3i) wrt parameter %3i at tfinal = %13.6e \n', i, k, dxdp(i,k,nt))
            end
        end
        fprintf('\n');
    end
    
    % Prints sensitivity of solution if i_plot_sensitivity = 1
    if (iplot.sensitivities == 1 && iplot.dxdz)
        for i=1:xdim
            for k=1:xdim
                fprintf('Sensitivity of x(%3i) wrt initial condition %3i at tfinal = %13.6e \n', i, k, dxdp(i,kdim+k,nt))
            end
        end
        fprintf('\n');
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     Plot solutions and sensitivities                                                             %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Plot sensitivities wrt variables and parameters in lists
    if (iplot.sensitivities == 1)
        plot_sensitivities(output_directory,imap,t,x,dxdp,xdim,kdim,ilist,klist,iplot)
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     Quantity of Interest                                                                         %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % Calculate quantity of interest
    [q,dqdparam]=solve_qoi(t,x,[p;x0],dxdp,xdim,kdim);
    
    % Print quantity of interest sensitivities wrt parameters at tfinal
    if (iplot.sensitivities == 1 && iplot.dqdp == 1)
        fprintf('qoi at tfinal = %13.6e \n', q(nt))
        for k=1:kdim
            fprintf('Sensitivity of qoi wrt parameter %3i at tfinal = %13.6e \n', k, dqdparam(k,nt))
        end
        fprintf('\n');
    end
     
    % Print quantity of interest sensitivities wrt initial conditions at tfinal
    if (iplot.sensitivities == 1 && iplot.dqdz == 1)
        fprintf('qoi at tfinal = %13.6e \n', q(nt))
        for k=1:xdim
            fprintf('Sensitivity of qoi wrt initial condition %3i at tfinal = %13.6e \n', k, dqdparam(kdim+k,nt))
        end
        fprintf('\n');
    end
    
    
    % Plot sensitivities wrt qoi
    if (iplot.sensitivities == 1)
        plot_qoi_sensitivities(output_directory,imap,t,q,dqdparam,xdim,kdim,ilist,klist,iplot)
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     R0                                                                                           %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if NextGen(1) ~= 0
        % Compute R0 and its sensitivities
        DFx0 = x(:,1);     % create a disease-free initial condition
        for i=NextGen
            DFx0(i)= 0;
        end
        [r0,dR0dp,R0warnings] = solve_r0(DFx0,p,dxdp(:,1:kdim,end),kdim,NextGen);   % Not quite correct
        
        % See if there is a condition in the theorem that does not hold!
        if R0warnings.w0 == 1
            fprintf('R0 sensitivities are not available:\nMuPad is unable to calculate the analytical R0 from the Next Generation Matrix.\nPerhaps there are too many infected classes.\nFor efficiency, set R0_only = 1.\n');
            msgbox('WARNING: R0 sensitivities are not available.  MuPad is unable to calculate the analytical R0 from the Next Generation Matrix.  Perhaps there are too many infected classes.  For efficiency, set R0_only = 1.');
            % force R0_only = 0 for the remainder of run
            R0_only = 1;
        end
        if R0warnings.w1 == 1
            fprintf('R0 is NOT VALID:\nThe transition matrix is not asymptotically stable\n');
            msgbox('WARNING: R0 is NOT VALID.  The transition matrix is not asymptotically stable.');
        end
        if R0warnings.w2F == 1
            fprintf('R0 is NOT VALID:\nThe fecundity matrix F is not nonnegative\n');
            msgbox('WARNING: R0 is NOT VALID.  The fecundity matrix F is not nonnegative.');
        end
        if R0warnings.w2V == 1
            fprintf('R0 is NOT VALID:\nThe transition matrix T is not nonnegative\n');
            msgbox('WARNING: R0 is NOT VALID.  The transition matrix T is not nonnegative.');
        end
        if R0warnings.w3 == 1
            fprintf('R0 is NOT VALID:\nThe equilibrium is not asymptotically stable in the absence of disease\n');
            msgbox('WARNING: R0 is NOT VALID.  The equilibrium is not asymptotically stable in the absence of disease.');
        end
        if r0<0
            fprintf('R0 is NOT VALID:\nR0 must be a nonnegative number\n');
            msgbox('WARNING: R0 is NOT VALID.  R0 must be a nonnegative number.');
        end
        if R0warnings.w4 ~= 0
            if R0warnings.w4 == 100 
                fprintf('R0 is NOT VALID:\nThe transition matrix V singular.\nPerhaps SENSAI does not correctly identify the placement of terms in F and V.\n');
                msgbox('WARNING: R0 is NOT VALID.  The transition matrix V singular.  Perhaps SENSAI does not correctly identify the placement of terms in F and V.');
            else        
                fprintf('R0 is NOT VALID:\nEquation %d is identically 0 and makes the transition matrix V singular\n', R0warnings.w4);        
                msgbox(['WARNING: R0 is NOT VALID.  Equation ' num2str(R0warnings.w4),' is identically 0 and makes the transition matrix V singular.']);
            end
        end
        if R0warnings.w5 ~= 0 
            fprintf('R0 is NOT VALID:\nThe disease-free subspace is not invariant.  Disease can enter state %d even in the absence of infection\n',R0warnings.w5);
            msgbox(['WARNING: R0 is NOT VALID.  The disease-free subspace is not invariant.  Disease can enter state ',num2str(R0warnings.w5), ' even in the absence of infection.']);
        end
        
        % Print R0 and its sensitivities
        fprintf('\nR0 = %5.4d \n', r0)
        
        if iplot.sensitivities == 1 && R0_only == 0
            for k=1:kdim
                fprintf('Sensitivity of R0 wrt parameter %3i = %13.6e \n', k, dR0dp(k))
            end
            fprintf('\n');
        end
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     User-defined parameter                                                                       %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %   Calculate sensitivities of variable and qoi wrt "new" parameter
    [cparam,dxdc,dqdc]=solve_new_param_sensitivities(t,x,[p;x0],dxdp,dqdparam,xdim,kdim);
    if NextGen(1) ~= 0
        [dR0dc] = solve_new_param_sensitivities_R0(x0,p,dR0dp,xdim,kdim);
    end
    
    % Print/plot sensitivity of variables and qoi wrt "new" parameter
    if (iplot.sensitivities == 1 && iplot.cp == 1)
        fprintf('\nSensitivity wrt user defined parameter \n')
        for i=1:xdim
            fprintf('Sensitivity of x(%3i) wrt new parameter  at tfinal = %13.6e \n', i, dxdc(i,nt))
        end
        fprintf('Sensitivity of qoi wrt new parameter at tfinal = %13.6e \n\n', dqdc(nt))
        if NextGen(1) ~= 0 && R0_only == 0
            fprintf('Sensitivity of R0 wrt new parameter at t=0 is %13.6e \n\n', dR0dc)
        end
        
        % Plot sensitivity of variables and qoi wrt "new" parameter
        plot_new_param_sensitivities(output_directory,imap,t,dxdc,dqdc,xdim,kdim,ilist);

    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                                                                                                  %
    %     Elasticities                                                                                 %
    %                                                                                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % Calculate elasticities
    [elxp,elqp]=solve_elasticities(t,x,p,x0,dxdp,q,dqdparam,xdim,kdim,iplot);
    [elxc,elqc]=solve_new_param_elasticities(t,x,q,cparam,dxdc,dqdc,xdim,kdim,iplot);
    if NextGen(1) ~= 0
        [elR0p] = solve_R0_elasticities(p,r0,dR0dp,kdim,iplot);
        [elR0c] = solve_new_param_elasticities_R0(r0,cparam,dR0dc,iplot);
    end
    
    % Print elasticity to parameters
    if(iplot.elasticities == 1 && (iplot.dxdp.true == 1 || iplot.dxdp.var == 1 || iplot.dxdp.param ==1 ))
        fprintf('\n');
        for i=1:xdim
            for k = 1:kdim
                fprintf('Elasticity of x(%3i) to parameter %3i at tfinal = %13.6e \n', i, k, elxp(i,k,nt))
            end
        end
        fprintf('\n');
    end
    
    % Print elasticity to initial conditions
    if(iplot.elasticities == 1 && iplot.dxdz == 1)
        for i=1:xdim
            for k = 1:xdim
                fprintf('Elasticity of x(%3i) to initial condition %3i at tfinal = %13.6e \n', i, k, elxp(i,kdim+k,nt))
            end
        end
        fprintf('\n');
    end
    
    % Print QoI elasticity to parameters
    if(iplot.elasticities == 1 && iplot.dqdp == 1)
        for k=1:kdim
            fprintf('Elasticity of qoi to parameter %3i at tfinal = %13.6e \n', k, elqp(k,nt))
        end
        fprintf('\n');
    end
    
    % Print QoI elasticity to initial conditions
    if(iplot.elasticities == 1 && iplot.dqdz == 1)
        for k=1:xdim
            fprintf('Elasticity of qoi to initial condition %3i at tfinal = %13.6e \n', k, elqp(k+kdim,nt))
        end
        fprintf('\n');
    end
    
    % Print R0 elasticity to parameters
    if NextGen(1) ~= 0  && R0_only == 0
        if iplot.elasticities == 1
            for k=1:kdim
                fprintf('Elasticity of R0 to parameter %3i at t=0 is %13.6e\n', k, elR0p(k))
            end
            fprintf('\n');
        end
    end
    
    % Print elasticity of variables and QoI wrt user-defined parameter
    if (iplot.elasticities == 1 && iplot.cp == 1)
        fprintf('\nElasticity wrt user defined parameter \n')
        for i=1:xdim
            fprintf('Elasticity of x(%3i) wrt new parameter  at tfinal = %13.6e \n', i, elxc(i,nt))
        end
        fprintf('Elasticity of qoi wrt new parameter at tfinal = %13.6e \n\n', elqc(nt))
    end
    
    % Print elasticity of R0 wrt user-defined parameter
    if (iplot.elasticities == 1 && iplot.cp == 1 && NextGen(1)~=0 && R0_only == 0)
        fprintf('Elasticity of R0 wrt new parameter at t=0 is %13.6e \n\n', elR0c)
    end
    
    % Plot elasticities with respect to parameters
    if (iplot.elasticities == 1)
        plot_elasticities(output_directory,imap,t,x,elxp,xdim,kdim,ilist,klist,iplot)
    end
    
    % Plot elasticities with respect to QoI
    if (iplot.elasticities == 1)
        plot_qoi_elasticities(output_directory,imap,t,q,elqp,xdim,kdim,ilist,klist,iplot)
    end
    
    % Plot elasticities with respect to user-defined parameter
    if(iplot.elasticities == 1 && iplot.cp == 1)
        plot_new_param_elasticities(output_directory,imap,t,elxc,elqc,xdim,kdim,ilist);
    end
    
    % Plot R0 information in a bar plot?
    if NextGen(1) ~= 0
        plot_R0(output_directory,r0,dR0dp,elR0p,dR0dc,elR0c,klist,iplot,R0_only)
    end
    
else
    
    %Display a warning if solution_only selected but NextGen indices are specified
    if solution_only == 1 && NextGen(1) ~= 0
        fprintf('Warning:  You have specified solution only but also to compute R0.  To save time, set NextGen = 0.\n');
        msgbox('Warning:  You have specified solution only but also to compute R0.  To save time, set NextGen = 0.');
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                  %
%     Combinational Warning Options                                                               %
%                                                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot_warnings(iplot,solution_only)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                  %
%     Save input data and solutions                                                                %
%                                                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename=[output_directory,'/inputs.txt'];  
fid = fopen(filename, 'wt');
fprintf(fid,'\nInitial solution \n');
for i = 1:xdim  
  fprintf(fid, 'x0(%3i,1) = %13.6e \n', i, x0(i,1));
end 

fprintf(fid,'\nParameters\n');
for k = 1:kdim  
  fprintf(fid, 'p(%3i,1) = %13.6e \n', k, p(k,1));
end 

fprintf(fid,'\nFinal time \n');
fprintf(fid, 'tfinal = %13.6e \n', tfinal);
fclose(fid);

copyfile('gvec.m',output_directory);
if(solution_only == 0)
    copyfile('dgvec_dxvec.m',output_directory);
    copyfile('dgvec_dparam.m',output_directory);
    copyfile('qoi.m',output_directory);
    copyfile('cp.m',output_directory);
end

filename=[output_directory,'/output.mat'];   
if(solution_only == 0)
    if(NextGen(1) ~= 0)
        save(filename, 't', 'x', 'p', 'dxdp', 'q', 'dqdparam', 'elxp', 'elqp', 'xdim', 'kdim', 'tfinal','cparam','dxdc','dqdc','elxc','elqc','r0','dR0dp','dR0dc','elR0p','elR0c')
    else
        save(filename, 't', 'x', 'p', 'dxdp', 'q', 'dqdparam', 'elxp', 'elqp', 'xdim', 'kdim', 'tfinal','cparam','dxdc','dqdc','elxc','elqc')
    end
else
    save(filename, 't', 'x', 'p', 'xdim', 'kdim', 'tfinal')
end

% % Special version for pine12
% if ~isempty(strfind(DIR,'pine12b'))
%     pine12_raw_output(output_directory)
% end
 
end






