% Daniel J. Bates (IMA), October 19, 2007

% This file contains a quick and dirty interface between Matlab and Bertini.
% I wrote part of this script for a specific example for a project in optimal control (with A. Becutti, I. Fotiou, and M. Morari, all of ETH).
% The point of this sample is not to provide a complete interface (that will come later) but rather to indicate how to use Bertini from within Matlab.
% If you take the time to generalize this m-file (especially if you make use of the symbolic toolbox of Matlab), I would be delighted to hear about it!
% Feel free to ask questions -- my current email address is dbates1@nd.edu (although dabates@ima.umn.edu might work as well).

% NOTE:  This example uses parameter homotopies to solve systems in a parameterized family.  If you want to solve an isolated polynomial 
% system (no parameterized family) or are unsure of this means, you should begin with the other Matlab/Bertini file on my website.

% Please note that I am chopping out the guts of this file since the related optimal control problem does nothing to explain the Bertini interface.
% In particular, I am providing you with a function for executing a particular run of a parameter homotopy but no call to it.
% This script originally contained a long while loop with multiple calls to Bertini.  The idea of the control algorithm is to observe a parameter value, 
% call Bertini to solve a corresponding problem with that parameter value, and then (back in the while loop) use the solutions from Bertini to 
% come up with the next parameter value.
% If you are interested in the optimal control algorithm and the deleted code, please let me know!

% First, we know that in this case, there will always be 17 solutions of interest at the end of each run.  
% You will see later how this can be detected on the fly. 
Npaths = 17;     

% Here is the function of interest.
% We pass in the parameter values (x) and the number of interesting solutions we expect at the end of each run (Npaths).
% The function returns S, which contains the solutions from Bertini.
function [S] = sub_getsol(x,Npaths)

% Here are the random complex parameter values at which we solved the problem initially.
% The idea of a parameter homotopy is to move from a generic point in the parameter space (the one corresponding to xrand) to 
% the point in the parameter space currently of interest (the one corresponding to x).
% If you look in the other Matlab/Bertini interface file on my website, you will see that these are the numbers that we used there:
xrand = [- 5.7623 + 28.763*i; 
         - 333.88 + 51.269*i;
         7.88 - 5.269*i;
         8.31 - 31.269*i;
         3.4 + 1.269*i;
         33.88 + 51.269*i];

% Bertini I/O is largely file-based, so we open the "input" file to record our homotopy:
fid = fopen('input', 'w');  

% As in the other file, we have a number of tolerances that may be set.
% The default may be used for any of these EXCEPT USERHOMOTOPY, which in this case must be specified to 1 (since we are doing a parameter homotopy).
fprintf(fid, sprintf('CONFIG\n\nUSERHOMOTOPY: 1;\n'));
fprintf(fid, sprintf('TRACKTOLBEFOREEG: 1e-3;\n'));
fprintf(fid, sprintf('TRACKTOLDURINGEG: 1e-4;\n'));
fprintf(fid, sprintf('FINALTOL: 1e-5;\n'))
fprintf(fid, sprintf('IMAGTHRESHOLD: 1e-4;\n'));
fprintf(fid, sprintf('END;\n\n'));

% Again, as in the other file, we must declare all of the names.
% The big differences between this and the example in the other file are that:
%  1. Bertini does not homogenize in the case of parameter homotopies, so we use the type "variable" here, rather than "variable_group;"
%  2. We use parameters (not too surprising); and 
%  3. We have a path variable upon which the parameters depend.  This is assumed to go from 1 (at which we have the random complex parameter values) to 
%     0 (at which we have the target parameter values.
fprintf(fid, sprintf('INPUT\n\n'));
fprintf(fid, sprintf('variable u1, u2, lambda1, lambda2, lambda3, lambda4, mu1, mu2, mu3, mu4, mu5, mu6, iL2, iL3, vo2, vo3;\n'));
fprintf(fid, sprintf('constant iL11, vo11, ro1, uold1, iLmax1, vref1;\n'));
fprintf(fid, sprintf('constant iL12, vo12, ro2, uold2, iLmax2, vref2;\n'));
fprintf(fid, sprintf('parameter iL1, vo1, ro, uold, iLmax, vref;\n'));
fprintf(fid, sprintf('pathvariable t;\n'));
fprintf(fid, sprintf('function f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16;\n\n'));

% Here we initialize the constants relating to our starting and ending parameter values.
% Here are the starting parameter values (our random complex ones, stored in xrand).
fprintf(fid, sprintf('iL11 =  %g+(%g)*I;\n', real(xrand(1)), imag(xrand(1)) ));
fprintf(fid, sprintf('vo11 =  %g+(%g)*I;\n', real(xrand(2)), imag(xrand(2)) ));
fprintf(fid, sprintf('ro1 =  %g+(%g)*I;\n', real(xrand(3)), imag(xrand(3)) ));
fprintf(fid, sprintf('uold1 = %g+(%g)*I;\n', real(xrand(4)), imag(xrand(4)) ));
fprintf(fid, sprintf('iLmax1= %g+(%g)*I;\n', real(xrand(5)), imag(xrand(5)) ));
fprintf(fid, sprintf('vref1 = %g+(%g)*I;\n\n', real(xrand(6)), imag(xrand(6)) ));

% Here are the ending parameter values (passed in).
fprintf(fid, sprintf('iL12 =  %g;\n', x(1) ));
fprintf(fid, sprintf('vo12 =  %g;\n', x(2) ));
fprintf(fid, sprintf('ro2 =  %g;\n', x(3) ));     % Here, x is just a col. vec.
fprintf(fid, sprintf('uold2 = %g;\n', x(4) ));     % holding the last system state
fprintf(fid, sprintf('iLmax2= %g;\n', x(5) ));     % --see function passing argument.
fprintf(fid, sprintf('vref2 = %g;\n', x(6) ));

% Finally, we glue the two constants together in a homotopy, stored in the "parameters."
% WARNING:  In this sample file, we did not use the "gamma trick" as described in the book of Sommese & Wampler.
% It is important to move through the complex numbers in the parameter space in order to avoid singularities (with probability one).
% The gamma trick is a way of ensuring that paths in the parameter space will be complex and random.
% Since we are beginning at a random complex point in the parameter space in this example, we opted not to use the gamma trick.
% Be very careful about this if you use parameter homotopies, especially if going from real points to real points.
% Feel free to ask if you have questions.
fprintf(fid, sprintf('iL1 = t*iL11 + (1-t)*iL12;\n'));
fprintf(fid, sprintf('vo1 = t*vo11 + (1-t)*vo12;\n'));
fprintf(fid, sprintf('ro = t*ro1 + (1-t)*ro2;\n'));
fprintf(fid, sprintf('uold = t*uold1 + (1-t)*uold2;\n'));
fprintf(fid, sprintf('iLmax = t*iLmax1 + (1-t)*iLmax2;\n'));
fprintf(fid, sprintf('vref = t*vref1 + (1-t)*vref2;\n\n'));

% Finally, here are the functions themselves (which are actually identical to those in the other file, since the homotopy in this case 
% is wrapped into the parameters).
fprintf(fid, sprintf('f1=0.4*u1-0.2*uold-0.2*u2+lambda1+mu1*(2-2*vo1)+mu2-mu3;\n'));
fprintf(fid, sprintf('f2=-0.5*lambda1+3.52*lambda2*ro+0.48*lambda3-2.52*lambda4*ro+mu4;\n'));
fprintf(fid, sprintf('f3=20*vo2-20*vref+lambda2*(-352*ro-3.52)-lambda3+lambda4*(352*ro+2.52)-2*mu4*u2;\n'));
fprintf(fid, sprintf('f4=0.2*u2-0.2*u1+lambda3+mu4*(2-2*vo2)+mu5-mu6;\n'));
fprintf(fid, sprintf('f5=-0.5*lambda3+3.52*lambda4*ro;\n'));
fprintf(fid, sprintf('f6=20*vo3-20*vref+lambda4*(-352*ro-3.52);\n'));
fprintf(fid, sprintf('f7=0.48*iL1-vo1+u1-0.5*iL2;\n'));
fprintf(fid, sprintf('f8=-2.52*iL1*ro-352*ro*vo2+352*ro*vo1+3.52*ro*iL2+2.52*vo1-3.52*vo2;\n'));
fprintf(fid, sprintf('f9=0.48*iL2-vo2+u2-0.5*iL3;\n'));
fprintf(fid, sprintf('f10=-2.52*ro*iL2-352*ro*vo3+352*ro*vo2+3.52*ro*iL3+2.52*vo2-3.52*vo3;\n'));
fprintf(fid, sprintf('f11=mu1*(iL1-iLmax+2*u1-2*vo1*u1);\n'));
fprintf(fid, sprintf('f12=mu2*(u1-1);\n'));
fprintf(fid, sprintf('f13=-mu3*u1;\n'));
fprintf(fid, sprintf('f14=mu4*(iL2-iLmax+2*u2-2*vo2*u2);\n'));
fprintf(fid, sprintf('f15=mu5*(u2-1);\n'));
fprintf(fid, sprintf('f16=-mu6*u2;\n\nEND;\n\n'));

fclose(fid);    % so that Bertini can have access to the file

% Now we are done writing the input file.  There is also a "start" file that must be created, as described and illustrated at the end of 
% the other file.
! ./bertini input

% Finally, we grab the endpoints from the appropriate output file, as we did in the other interface file.
fid = fopen('real_solutions', 'r');
S = [];
Nvars = fscanf(fid, '%d\n', [1,1]);  % first line of file is # variables
tmp = fscanf(fid, '%d\n', [1,1]);  % second line is rubbish
for j = 1:Npaths  % for each finite solution, grab the path number and then the (complex) coordinates
  path_num = fscanf(fid, '%d\n', [1,1]);  % path number is meaningless unless = -1 (implies EOF)
  if (path_num == -1)  % indication that we have found the end of the file - this is where Npaths could be exchanged for a while loop....
    fprintf('ERROR: UNEXPECTED NUMBER OF FINITE SOLUTIONS!\n');
  end 
  point = [];
  for k = 1:Nvars  % just grab the Nvars coordinates
    point = fscanf(fid, '%e %e\n', [1,2]);
    S(k,j) = point(1) + point(2)*i;
  end
end
fclose(fid);

return