%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Copyright (c) 2004.  Takayuki Osogami.  All Rights Reserved.   %
%                                                                   %
% Permission to use, copy, modify, and distribute this source code  %
% and its documentation for any purpose, without fee, and without   %
% a written agreement, is hereby granted, provided that the above   %
% copyright notice, this paragraph and the following two paragraphs %
% appear in all copies, modifications, and distributions.           %
% Created by Takayuki Osogami, Department of Computer Science,      %
% Carnegie Mellon University.                                       %
%                                                                   %
% IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR ANY       %
% DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS   %
% SOURCE CODE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN    %
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                        %
%                                                                   %
% THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING,      %
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY     %
% AND FITNESS FOR A PARTICULAR PURPOSE.  THE SOURCE CODE AND        %
% ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS         %
% PROVIDED "AS IS."  THE AUTHOR HAS NO OBLIGATION TO PROVIDE        %
% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage: [tau,T] = matching3PH(momentsofG)
%
% INPUT, momentsofG, is a vector of size 3.
% An example of an input is:
%
%  momentsofG = [2,6,24];
%
% OUTPUT, [tau,T], is a pair of a vector and a matrix.
% When the input is [2,6,24], the output is
%
%  tau = [1,0];
%
%  T = [-1  1;
%        0 -1];
%
% Yes.  This is an Erlang-2 distribution with mean 2.
%
% Given first three moments, momentsofG,
% the function returns a PH distribution, (tau,T), whose
% first three moments is momentsofG.
% The output PH distribution has no mass probability at zero.
% That is the sum of all the elements of tau is 1.
%
% The function is not defined everywhere, but almost everywhere.
% When the function is not defined for the input, it perturbs  
% the input and returns a PH distribution that has the first three
% moments close to the input.  The error is very small for most
% practical purposes.
%
% When the output is (0,0), there was an error in the function
% or the input was illegal.
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: June 6, 2004

function [tau,T] = matching3PH(momentsofG)

%% Maximum number of phases used in the PH distribution.
%% If the number of necessary phases exceeds this limit,
%% the input is matched by the closest PH distribution.
MAX_PHASE = 50;

EPSILON = 0.001;

%% If the input moment is zero, return error.

if( (momentsofG(1) == 0) || (momentsofG(2) == 0) )

  Illegal_input = momentsofG
  tau = 0;
  T = 0;
  return;

end

%% Calculate the normalized moments of the input, G.

m2 = momentsofG(2) / momentsofG(1)^2;
m3 = momentsofG(3) / (momentsofG(1) * momentsofG(2));

%% Any PH distribution has m2 > 1 and m3 > m2.
%% If there is no PH for the input, return error.

if( (m2 <= 1 + EPSILON^2) || (m3 <= m2 + EPSILON^2) )

  Illegal_input = momentsofG
  tau = 0;
  T = 0;
  return;

end

%% If the necessary number of phases exceeds the limit, MAX_PHASE,
%% calculate the normalized moments of the "closest" PH distribution.

smallestM2 = (MAX_PHASE/(MAX_PHASE-1) + (MAX_PHASE-1)/(MAX_PHASE-2)) / 2;
if( m2 < smallestM2 )
  m2 = smallestM2;
end

smallestM3 = (MAX_PHASE+2)/(MAX_PHASE+1) * m2;
if( m3 < smallestM3 )
  m3 = smallestM3;
end

%% If there is no PH distribution with the input moments,
%% calculate the normalized moments of the "closest" PH distribution.

if( (m3 < 2*m2 - 1) && (abs(m3 - 3/2 * m2) < EPSILON^2) )

  m3 = (1+EPSILON) * m3;

elseif( m3 > 2*m2 - 1 )

 r = m3 / (2*m2-1);
 n = floor(1/(m2-1));

 if( abs(1/(m2-1) - n) < EPSILON^2 )
   m2 = 1 + 1 / (n*(1-EPSILON));
   m3 = r * (2*m2-1);
 end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is exponential %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( (abs(m2-2) < EPSILON^2) && (abs(m3-3) < EPSILON^2) )

  tau = 1;
  T = - 1 / momentsofG(1);
  return;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is 2-phase PH %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( isPH2(m2,m3) )

  [tau,T] = matching3PH2(momentsofG(1),m2,m3);
  return

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is EC (Erlang-Coxian) %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( m3 >= 2*m2 - 1 )

  [tau,T] = matching3EC(momentsofG(1),m2,m3);
  return;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is "p * EC + (1-p) * EXP" %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = floor((2*m2-m3)/(m3-m2) + EPSILON^2);

if( m3 >= ((k+1)*m2 + (k+4)) / (2*(k+2)) * m2 )

  y = (2-m2) / (4*(3/2-m3/m2));
  p = (2-m2)^2 / ((2-m2)^2+4*(2*m2-1-m3));
  m2X = 2 * y;
  m3X = 2 * m2X - 1;

  [tauX,TX] = matching3EC(1,m2X,m3X);
  [tau,T] = mixtureofPH(tauX,TX,1,-1/y,p);

  %% normalize
  mean = momentofPH(tau,T,1);
  T = T * mean / momentsofG(1);

  return;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is "EC + EXP"   %%
%% EC has mass probability at zero %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( abs(m2-2) < EPSILON^2 )

  y = (m3-2*(k+3)/(k+2)) / (3-m3);
  m2X = 2*(1+y);

else

  A = m2 * ((m3-3) - 2*(k+3)/(k+2)*(m2-2));
  B = m2 * sqrt((m3-3)^2 + 8*(k+3)/(k+2)*(m2-2)*(3/2-m3/m2));
  C = 2 * (k+3)/(k+2) * (m2-2)^2;
  y = (A+B) / C;
  m2X = (1+y) * (m2*(1+y)-2*y);

  %% This is a magic spell for Matlab, who doesn't like unnecessary precision...
%%  m2X = round(m2X * 1000) / 1000;

end
  
m3X = (k+3)/(k+2) * m2X;

[tauX,TX] = matching3EC(1,m2X,m3X);

[tau,T] = convolutionofPH(tauX,TX,1,-1/y);

%% normalize
mean = momentofPH(tau,T,1);
T = T * mean / momentsofG(1);

return;