function [mh_dist,mstr]=MHdist(x,f,g,n,infty)
% calculate MH distance between ccdfs
arguments
    x (1,:) {mustBeNumeric,mustBeReal}
    f (1,:) {mustBeNumeric,mustBeReal}%,mustBeEqualSize(f,x)}
    g (1,:) {mustBeNumeric,mustBeReal}%,mustBeEqualSize(g,x)}
    n (1,1) {mustBeNumeric,mustBeReal}=1000
    infty (1,1) logical=true
end
% checks:
if any(diff(x)<=0) | any(x<0)
    eid = 'xRange:invalid';
        msg = 'x values invalid.';
    throwAsCaller(MException(eid,msg))
end
if any([diff(f) diff(g)]>0)
    eid='function:nondecreasing'
    msg('function not decreasing.');
    throwAsCaller(MException(eid,msg));
end
if max([f(1),g(1)])>1 | min([f(end),g(end)])<0
    eid='function:notvalid'
    msg('exceeds valid ccdf range');
    throwAsCaller(MException(eid,msg))
end
% map and add end points if needed:
m=max(x);
if infty 
    t=x./(x+1); % map to [0,1] using MH
    t(end+1)=1;
    f(end+1)=0;
    g(end+1)=0;
else
    t=x/m; % scale
end
if t(1)>0
    t=[0 t];
    f=[1 f];
    g=[1 g];
end
T=linspace(0,1,n);
F=interp1(t,f,T);
G=interp1(t,g,T);
mh_dist=sum(abs(F-G))/n;
intl=[2 0.05 0.02 0.01 0.005 0.002];
i=find(mh_dist<intl,1,'last');
switch i
    case 1
        mstr='bad';
    case 2
        mstr='mediocre';
    case 3
        mstr='acceptable';
    case 4
        mstr='good';
    case 5
        mstr='excellent';
    otherwise
        mstr='perfect';
end
end



