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

% This file contains a quick and dirty interface between Matlab and Bertini.
% I wrote 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 an automatically-generated homotopy rather than a parameter homotopy.  Another file on my website covers parameter homotopies.

function [x, uopt] = duf_closedloop_2_step_1()
clc

% This project (like many others) involved a two-stage process for solving a parameterized family of homotopies.
% In this (first) stage, we choose random complex values for the parameters and solve the system from scratch.
% The book by A. Sommese & C. Wampler provides an excellent description of the underlying ideas, although you could also look in the 
% paper (on my website) in which this example is described.
% Here are the "random" parameter values for our example (chosen at random elsewhere - this could probably be automated within Matlab).
iL1 =  -5.7623 + 28.763*i;
vo1 =  -333.88 + 51.269*i;
ro =   7.88 - 5.269*i;
uold = 8.31 -31.269*i;
iLmax= 3.4 +1.269*i;
vref = 33.88 + 51.269*i;

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

% Next we actually write to the file.
% Various configurations go at the top of the file.  Please see the Bertini user's manual for a description of the various tolerances.
% All tolerances have reasonable defaults, so this line could be entirely neglected for "vanilla" runs:
fprintf(fid, sprintf('CONFIG\n\nTRACKTOLBEFOREEG: 1e-3;\nTRACKTOLDURINGEG: 1e-4;\nFINALTOL: 1e-5;\nPRINTPATHMODULUS: 20;\nIMAGTHRESHOLD: 1e-4;\nRANDOMSEED: 2;\n\nEND;\n\n'));

% The next part of the input file records declarations of the various names to be used in the functions.  For example, in this 
% case, we have variables u1, u2, ..., vo3 (broken into two groups), several constants, and 16 functions. 
% You can choose just about any names that you want (alphanumeric, starting with a constant, of course).
fprintf(fid, sprintf('INPUT\n\nvariable_group u1, u2, lambda1, lambda2, lambda3, lambda4, mu1, mu2, mu3, mu4, mu5, mu6;\nvariable_group iL2, iL3, vo2, vo3;\nconstant iL1, vo1, ro, uold, iLmax, vref;\nfunction f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16;\n\n'));

% Next I need to initialize my constants (i.e., the random complex parameter values).  If you just want to include your 
% constants in the function definitions, that is fine, too; the next line could be omitted in that case.
fprintf(fid, sprintf('iL1 = %g+(%g)*I;\nvo1 = %g+(%g)*I;\nro = %g+(%g)*I;\nuold = %g+(%g)*I;\niLmax = %g+(%g)*I;\nvref = %g+(%g)*I;\n\n', real(iL1), imag(iL1), real(vo1), imag(vo1), real(ro), imag(ro), real(uold), imag(uold), real(iLmax), imag(iLmax), real(vref), imag(vref)));  

% Finally, we have the functions themselves:
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



% We are now done with the input file and can call Bertini.  Bertini will solve the system for us.  For details on how it does so, 
% please refer to either the user's manual or the paper on a "paradigm for software in numerical algebraic geometry" by 
% Bates, Hauenstein, Sommese, and Wampler.
% Here is how you call Bertini from within Matlab.  Of course, this assumes that the Bertini executable is in the current working directory.
% If not, either move it there or add the appropriate path before "bertini" on the next line:
! ./bertini input



% Assuming that all went well (which it should have for this example), we are ready to grab the solutions that interest us.
% Solutions are provided in several files by Bertini (finite, real, etc.).  Again, please see the user's manual for all of the various output files.
% In this case, we cared about the finite solutions, so we look in the file finite_solutions:
fid = fopen('finite_solutions', 'r');
% Next we pull in the finite solutions and store them in "S"
S = [];
% In this case, we new that the number of finite solutions would be 19, so we just hard-coded that in.
% If the number of finite solutions is not known a priori, a while loop teminating on a "path number" of -1 would suffice.
Npaths = 19;
Nvars = fscanf(fid, '%d\n', [1,1]);  % The first line of the file is the number of variables.
tmp = fscanf(fid, '%d\n', [1,1]);  % The second line is rubbish (for our purposes).
for j = 1:Npaths  % For each finite solution, we grab the path number and then the (complex) coordinates.
  path_num = fscanf(fid, '%d\n', [1,1]);  % The path number is meaningless unless it is -1 (which implies that there are no more solutions)
  if (path_num == -1)  % could be replaced with a while loop....
    fprintf('ERROR: UNEXPECTED NUMBER OF FINITE SOLUTIONS!\n');
  end 
  point = [];  % We will temporarily store our point in here.
  for k = 1:Nvars  % We just grab the Nvars coordinates.
    point = fscanf(fid, '%e %e\n', [1,2]); % The format is simple!
    S(k,j) = point(1) + point(2)*i; % The two parts are just the real and complex parts.
  end
end
fclose(fid);  % That's it - we're done reading in the finite solutions!


% This last code chunk sets up a "start" file, which contains all of our solutions from this run.
% If you just want to solve a single, isolated system (rather than a parameterized family), this is worthless.
% The point of this is that we need to know where to start our paths in a user-defined parameter homotopy.
% "start" contains the solutions to the initial system (which is exactly what we solved above).
fid = fopen('start', 'w');
fprintf(fid, '%d\n\n', Npaths);
for j = 1:Npaths
  for k = 1:Nvars
    fprintf(fid, '%e %e;\n', real(S(k,j)), imag(S(k,j)));
  end
  fprintf(fid, '\n');
end
fclose(fid);

return