map_mmpp2.m

function MAP=map_mmpp2(MEAN,SCV,SKEW,ACF1)
% MAP=map_mmpp2(MEAN,SCV,SKEW,ACF1) - Fit a MMPP(2) as a MAP
%
%  Input:
%  MEAN: mean inter-arrival time of the process
%  SCV: squared coefficient of variation of inter-arrival times
%  SKEW: skewness of inter-arrival times (-1 => automatic minimization,
%  applies only to SCV>1)
%  ACF1: lag-1 autocorrelation coefficient (-1 => maximum feasible
%  autocorrelation)
%
%  Output:
%  MAP: a MAP in the form of {D0,D1}
%
%  Examples:
%  - MAP=map_mmpp2(1,2,-1,0.2) MMPP(2) process with minimal skewness
%  - MAP=map_mmpp2(1,2,-1,-1) MMPP(2) process with minimal skewness and
%  maximal autocorrelation
%
 
SCV_REQ=SCV; % keep the request for diagnostics, SCV is recomputed below
FEASTOL=10^(-map_feastol);
 
E1=MEAN;
E2=(1+SCV)*E1^2;
E3=-(2*E1^3-3*E1*E2-SKEW*(E2-E1^2)^(3/2));
 
% The closed form below solves the moment-matching equations, but nothing in
% it constrains the solution to be a MAP: outside the MMPP(2) feasible set it
% returns negative rates, i.e. a D1 with negative entries and a D0 with a
% positive diagonal. Rowsums stay zero, so the usual generator check does not
% catch it. Reject the request instead of returning a non-MAP.
if SCV < 1-FEASTOL
    error('map_mmpp2: SCV=%g is infeasible, the inter-arrival times of an MMPP(2) are over-dispersed (SCV>=1).', SCV_REQ);
end
if abs(SCV-1) <= FEASTOL
    error('map_mmpp2: SCV=1 is the Poisson boundary, where the MMPP(2) fit is degenerate: the decay rate G2=ACF1/(1-1/SCV)/0.5 divides by zero and every rate comes back NaN. Use map_exponential(%g) for a Poisson process.', MEAN);
end
 
% ACF1=RHO0MAX is attained in the limit of a decay rate G2->1; no MMPP(2)
% exceeds it, and none is negatively autocorrelated
RHO0MAX=0.5*(1-1/SCV);
if ACF1~=-1
    if ACF1 < -FEASTOL
        error('map_mmpp2: ACF1=%g is infeasible, an MMPP(2) cannot be negatively autocorrelated. Pass ACF1=-1 to request the maximum feasible autocorrelation.', ACF1);
    end
    if ACF1 > RHO0MAX+FEASTOL
        error('map_mmpp2: ACF1=%g exceeds the maximum lag-1 autocorrelation %g feasible at SCV=%g. Pass ACF1=-1 to request it.', ACF1, RHO0MAX, SCV_REQ);
    end
end
 
if ACF1==-1
    G2=1-10*10^(-map_feastol); % autocorrelation decay rate
else
    G2=ACF1/(1-1/SCV)/0.5; % autocorrelation decay rate
end
if SKEW==-1 && SCV>1 % determine MAP with nearly minimum third moment
    E3=(3/2+0.001)*E2^2/E1;
end
 
% the SKEW==-1 branch above sits just above E3MIN, which is the infimum of the
% third moment over the class
E3MIN=(3/2)*E2^2/E1;
if E3 < E3MIN-FEASTOL
    error('map_mmpp2: SKEW=%g gives E3=%g, below the minimum third moment %g feasible at SCV=%g. Pass SKEW=-1 for the minimum-skewness fit.', SKEW, E3, E3MIN, SCV_REQ);
end
 
SCV=(E2-E1^2)/E1^2;
if G2<1e-6
    mu00=2*(6*E1^3*SCV-E3)/E1/(6*E1^3*SCV+3*E1^3*SCV^2+3*E1^3-2*E3);
    mu11=0;
    q01= 9*E1^5*(SCV-1)*(SCV^2-2*SCV+1)/(6*E1^3*SCV-E3)/(6*E1^3*SCV+3*E1^3*SCV^2+3*E1^3-2*E3);
    q10=-3*(SCV-1)*E1^2/(6*E1^3*SCV-E3);
else
    mu00=G2*(-4*E3*G2+4*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3*G2-18*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2-18*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV^2-12*E1^3*G2^2-12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV+12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV+12*E1^3*G2*SCV^2-9*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV+3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)+12*E1^3*G2^2*SCV+9*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2+12*E1^3*G2+12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2-3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^3)/(12*E1^3*G2^3*SCV+3*E1^3*SCV^3*G2-12*E1^3*G2^3+18*E1^3*G2^2*SCV^2-3*E1^3*G2+27*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV^2-9*E1^3*G2*SCV^2+18*E1^3*G2^2-12*E1^3*G2^2*SCV+9*E1^3*G2*SCV-12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^3*SCV-9*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^3*G2-24*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV^2-(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3*SCV^2+4*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3*G2^2+12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^3-(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3+2*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3*SCV+9*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2+24*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV-27*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV+6*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV-12*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2-24*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2+6*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^3-4*E3*G2^2)/E1;
    mu11=(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3);
    q01=-3*E1^2*(-6*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV+12*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV-6*G2*SCV*E1^2-3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2+(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3+3*E1^2*G2+6*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2-9*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2*G2+3*E1^2*G2*SCV^2-E3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV-6*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV+6*E1^2*G2^2*SCV+3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2-G2*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3-3*E1^2*G2^2+3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2*G2^2-3*E1^2*SCV^2*G2^2+G2*SCV*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3)/(-45*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV^2+18*G2^2*E1^5*SCV+18*E1^5*G2^3-27*E1^5*G2^2*SCV^2+6*E1^2*G2^2*E3-27*E1^5*G2^2-18*E1^5*G2^3*SCV-18*E1^5*G2*SCV+18*E1^5*G2*SCV^2+3*E1^2*G2*E3-3*E1^2*G2*E3*SCV+(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/E1/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3^2+3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV*E3-36*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV+36*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2+36*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2+45*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV-12*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV*E3-3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*E3+9*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV^3+36*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*SCV^2-6*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^2/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^2*E3+18*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^3*SCV-18*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2^3-9*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))*E1^5/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2);
    q10=3*(-3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^3-3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV^2+6*E1^3*SCV^2+3*E1^3*G2*SCV^2+3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV^2+6*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2*SCV-E3*SCV-6*E1^3*SCV+(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3*SCV-6*E1^3*G2*SCV-3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*SCV-(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*E3+3*E1^3*G2-3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)*G2+3*E1^3*(-3*E1^3*G2+3*E1^3*G2*SCV-6*E1^3*SCV+E3+(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2)^(1/2))/(-3*E1^3*SCV^2-6*E1^3*SCV-3*E1^3+2*E3)+E3)*E1^2*(-1+G2)/(E3^2-12*E1^3*SCV*E3+6*E1^3*G2*E3-6*G2*SCV*E1^3*E3+18*G2*SCV^3*E1^6-18*E1^6*G2*SCV^2+9*E1^6*G2^2+36*E1^6*SCV^2+18*E1^6*G2*SCV-18*E1^6*SCV*G2^2+9*E1^6*SCV^2*G2^2-18*E1^6*G2);
end
% Catch-all: the checks above cover the known infeasible directions, but the
% authoritative test is the solution itself. An MMPP(2) has non-negative rates
% by construction, so anything else is not a MAP and must not be returned.
RATES=[mu00,mu11,q01,q10];
if ~isreal(RATES) || any(isnan(RATES)) || any(RATES < -FEASTOL)
    error('map_mmpp2: (MEAN=%g,SCV=%g,SKEW=%g,ACF1=%g) is not MMPP(2)-feasible: the fit gives [mu00 mu11 q01 q10]=%s, which is not a MAP.', ...
        MEAN, SCV_REQ, SKEW, ACF1, mat2str(RATES,6));
end
% a request on the feasibility boundary lands on a zero rate up to roundoff:
% clear that sign flip, but leave small POSITIVE rates alone -- ACF1=-1 asks
% for G2=1-1e-7, whose near-uncoupled chain has legitimate rates around 1e-9
RATES(RATES<0)=0;
mu00=RATES(1); mu11=RATES(2); q01=RATES(3); q10=RATES(4);
 
D0=[ -mu00-q01,       q01;
    q10, -mu11-q10];
D1=[ mu00,   0;
    0, mu11];
MAP={D0,D1};
end