%{
Copyright (c) 2015 to 2016 by Cornell University and The Regents Of
The University Of California. All Rights Reserved.

Permission to use this Procedural Yarn Fitting and Generation Tool (the "Work")
and its associated copyrights solely for educational, research and non-profit
purposes, without fee is hereby granted, provided that the user agrees as
follows:

Those desiring to incorporate the Work into commercial products or use Work and
its associated copyrights for commercial purposes should contact the Center for
Technology Licensing at Cornell University at

395 Pine Tree Road, Suite 310, Ithaca, NY 14850;
email: ctl-connect@cornell.edu;
Tel: 607-254-4698;
FAX: 607-254-5454

for a commercial license.

IN NO EVENT SHALL CORNELL UNIVERSITY ("CORNELL") OR THE UNIVERSITY OF
CALIFORNIA ("UC") BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL,
INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF
THE USE OF THE WORK AND ITS ASSOCIATED COPYRIGHTS, EVEN IF CORNELL OR UC MAY
HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

THE WORK PROVIDED HEREIN IS ON AN "AS IS" BASIS, AND NEITHER CORNELL NOR UC HAS
ANY OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
MODIFICATIONS. CORNELL AND UC MAKE NO REPRESENTATIONS AND EXTEND NO WARRANTIES
OF ANY KIND, EITHER IMPLIED OR EXPRESS, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR
THAT THE USE OF WORK AND ITS ASSOCIATED COPYRIGHTS WILL NOT INFRINGE ANY PATENT,
TRADEMARK OR OTHER RIGHTS.
%}

% Likelihood maximization

function [b1 a1] = mleR(R, N)
	bvals = 0 : 0.1 : 5;
    avals = 0 : 0.02 : 0.2;
    res = zeros(length(bvals)*length(avals), 3);
    i = 0;
	
	for b = bvals
        for a = avals
            Z = integral(@(x) pR2(x, b, a), 0, 1);
            likelihood = sum(log(pR2(R, b, a, Z)));

            i = i + 1;
            res(i, :) = [b, a, likelihood];
        end
    end
	
	[~, id] = max(res(:, 3));
    fprintf('Initial guess : b = %.2f, a = %.2f (likelihood: %.4f)\n', ...
        res(id, 1), res(id, 2), res(id, 3));

    theta = mle(R, 'pdf', @(data, b, a) pR3(data, b, a), ...
        'start', [res(id, 1), res(id, 2)], 'lowerbound', [0, 0]);
    b1 = theta(1); a1 = theta(2);
    fprintf('Fitted        : b = %.2f, a = %.2f (likelihood: %.4f)\n', ...
        b1, a1, sum(log(pR3(R, b1, a1))));
		
		
	% Generate another N sample points using fitted parameters (i.e., b1,
    % a1) and store them in P1
    %
    P1 = zeros(2, N);
    for i = 1 : N
        while true
            R = sqrt(rand());
            if rand() < pR(R, b1, a1)
                break
            end
        end
        ang = 2*pi*rand();
        P1(1, i) = R.*cos(ang);
        P1(2, i) = R.*sin(ang);
    end
	
	%figure;	scatter(P1(1, :), P1(2, :), '.');
	%title('Fitted');
end