symMinimization.m

function [G, angle, Type, kappa] = symMinimization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath) %(Q_hom,B_eff)


syms v1 v2 q1 q2 q3 q12 b1 b2 b3


% -------- Options ----------
% % print_Input = false;  %effective quantities
% print_Input = true;
% % print_statPoint = false;
% print_statPoint = true;
% print_Minimizer = true;
% % make_FunctionPlot = false;
% make_FunctionPlot = true;

print_Uniqueness = false;         
check_StationaryPoints = false;   % could be added to Input-Parameters..
compare_with_Classification = false;  %maybe added later --- 


%fprintf('Functions to be minimized:')
f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2-2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
                                        + q12*(v1^2*v2^2-b2*v1^2-b1*v2^2+b1*b2);
f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2)...
                                        + q12*(v1^2*v2^2+b2*v1^2+b1*v2^2+b1*b2);


% ---- Fix parameters


% Epsilon used:
epsilon = 1e-8;


if ~exist('InputPath','var')
     % third parameter does not exist, so default it to something
      absPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs";
 end


% 1. Import effective quantities from CellSolver-Code:
%read as sparse Matrix...
try %absolutePath
    Qmat = spconvert(load(absPath + '' + "/QMatrix.txt"));
    Bmat = spconvert(load(absPath + '' + "/BMatrix.txt"));
%     fprintf('Use absolute Path')
catch ME  % use relativePath
    Qmat = spconvert(load('../outputs/QMatrix.txt'));
    Bmat = spconvert(load('../outputs/BMatrix.txt'));
%     fprintf('Use relative Path')
end
%convert to full matrix... 
Qmat = full(Qmat);
Bmat = full(Bmat);




% --- TODO CHECK: assert if Q is orthotropic ???  check ifq13=q31=q23=q32= 0 ?


if print_Input
    fprintf('effective quadratic form:')
    Qmat
    fprintf('effective prestrain')
    Bmat
    
    % check if Q is (close to..) symmetric
    % könnte Anti-symmetric part berechnen und schauen dass dieser klein?
    % Test: issymmetric(Qmat) does not work for float matrices?
    % symmetric part 0.5*(Qmat+Qmat')
    % anti-symmetric part 0.5*(Qmat-Qmat')
    if norm(0.5*(Qmat-Qmat'),'fro') < 1e-8
        fprintf('Qmat (close to) symmetric \n')
        norm(0.5*(Qmat-Qmat'),'fro') % TEST
    else
        fprintf('Qmat not symmetric \n')
    end
    
    
    % Check if B_eff is diagonal this is equivalent to b3 == 0
    if abs(Bmat(3)) < 1e-8
        fprintf('B_eff is diagonal (b3 == 0) \n')
    else
        fprintf('B_eff is  NOT diagonal (b3 != 0)  \n')
    end
end




% CAST VALUES TO SYM FIRST? This is done anyway..
% % Substitute effective quantitites      
f_plus = subs(f_plus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
                                                     Bmat(1), Bmat(2), Bmat(3)});

f_minus = subs(f_minus,{q1, q2, q3, q12, b1, b2, b3}, {Qmat(1,1), Qmat(2,2), Qmat(3,3), Qmat(1,2), ...
                                                     Bmat(1), Bmat(2), Bmat(3)});


% Compute the Gradients 
df_plusx = diff(f_plus,v1);
df_plusy = diff(f_plus,v2);
df_minusx = diff(f_minus,v1);
df_minusy = diff(f_minus,v2);

% Setup Equations Grad(f) = 0 
eq1 = df_plusx == 0;
eq2 = df_plusy == 0;
eqns_plus = [eq1, eq2];
eq3 = df_minusx == 0;
eq4 = df_minusy == 0;
eqns_minus = [eq3, eq4];

% -------  Symbolically Solve Equations: 
% More robust (works even for values b_3 ~ 1e-08 ): 
S_plus = solve(eqns_plus,v1,v2,'MaxDegree' , 5);  
S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5);

A_plus = S_plus.v1;
B_plus = S_plus.v2;
A_minus = S_minus.v1;
B_minus = S_minus.v2;

if check_StationaryPoints
    %---------- TEST if Grad(f) = 0 ---------------------
    fprintf('Testing equation grad(f) = 0  with stationary points')
    for i = 1:size(A_plus,1)
        fprintf('Testing %d.point (f_plus): ',i )
        [ double(subs(subs(df_plusx,v1,A_plus(i)),v2,B_plus(i))), double(subs(subs(df_plusy,v1,A_plus(i)),v2,B_plus(i))) ]
    end
    for i = 1:size(A_minus,1)
        fprintf('Testing %d.point (f_minus): ',i )
        [double(subs(subs(df_minusx,v1,A_minus(i)),v2,B_minus(i))), double(subs(subs(df_minusy,v1,A_minus(i)),v2,B_minus(i)))]
    end
    % ------------------------------------
end

% --- Extract only Real-Solutions
% fprintf('real stationary points of f_plus:')
tmp1 = A_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
tmp2 = B_plus(imag(double(A_plus))==0 & imag(double(B_plus)) == 0);
A_plus = tmp1;
B_plus = tmp2;
SP_plus = [A_plus,B_plus];

% fprintf('real stationary points of f_minus:')
tmp1 = A_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
tmp2 = B_minus(imag(double(A_minus))==0 & imag(double(B_minus)) == 0);
A_minus = tmp1;
B_minus = tmp2;
SP_minus = [A_minus,B_minus];

% TODO one should use f_plus.subs(A_plus..) to compute function value symbolically?
% in the end only the stationaryPoints are used.. should be ok to compare function values numerically

% Determine global Minimizer from stationary points:
% fprintf('function values at stationary points (f_plus):')
T_plus = arrayfun(@(v1,v2) double(f_plus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_plus,B_plus,'UniformOutput', false);
T_plus = cell2mat(T_plus);

%Test: use Substitution
% subs(f_plus,{v1, v2}, {A_plus,B_plus})

% fprintf('function values at stationary points (f_minus):')
T_minus = arrayfun(@(v1,v2) double(f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3)),A_minus,B_minus,'UniformOutput', false);
T_minus = cell2mat(T_minus);

%Test: use Substitution
% T_minus = subs(f_minus,{v1, v2}, {A_minus,B_minus})
% double(T_minus)

if print_statPoint
    fprintf('real stationary points of f_plus: ')
%     SP_Plus     %alternative: output as symbolic (can be unwieldy)
    double(SP_plus)
    fprintf('real stationary points of f_minus:')
%     SP_Minus     %alternative: output as symbolic (can be unwieldy)
    double(SP_minus)
    fprintf('function values at stationary points (f_plus):')
    T_plus
    fprintf('function values at stationary points (f_minus):')
    T_minus
end

% --- Find Stationary Point(s) with minimum Value 
[Min_plus,MinIdx_plus] = min(T_plus, [], 'all', 'linear');    %find one min...
[Min_minus,MinIdx_minus] = min(T_minus, [], 'all', 'linear');
% [Min_minus,MinIdx_minus] = min(T_minus)                     % works with symbolic too

% Compare Minimizers of f_plus & f_minus
[globalMinimizerValue,GlobalIdx] = min([Min_plus,Min_minus]); 

if GlobalIdx == 1     %Min_plus   % i.e. Global Minimizer is given by f_plus
    GlobalMinimizer = SP_plus(MinIdx_plus,:);
    sign = 1.0;
elseif GlobalIdx == 2  %Min_minus % i.e. Global Minimizer is given by f_minus
    GlobalMinimizer = SP_minus(MinIdx_minus,:);
    sign = -1.0; 
end

% ------ Check if there are more SP with the same value...
 MinIndices_minus = find(T_minus(:) == globalMinimizerValue);    % Find indices of All Minima
 MinIndices_plus  = find(T_plus(:) == globalMinimizerValue);    % Find indices of All Minima                   
 numMinSP_minus = size(MinIndices_minus,1);   % One of these is always >= 2 due to the structure of the roots..
 numMinSP_plus  = size(MinIndices_plus,1);
    
%     AllMinSP_minus = SP_minus(MinIndices_minus,:)
%     AllMinSP_minus = double(SP_minus(MinIndices_minus,:))
%     AllMin = T_minus(MinIndices)  %bereits klar dass diese selben funktionswert haben..
      
    Minimizer = sign*(GlobalMinimizer'*GlobalMinimizer);  % global minimizing Matrix G*
    MinimizerCount = 1;
    % different Stationary Points might correspond to the same minimizing
    % Matrix G*... check this: 
    
    % Compare only with other StationaryPoints/Minimizers
    % remove Index of Minimizer
    if GlobalIdx == 1     
        MinIndices_plus = MinIndices_plus(MinIndices_plus~=MinIdx_plus);
    elseif GlobalIdx == 2  
        MinIndices_minus = MinIndices_minus(MinIndices_minus~=MinIdx_minus);
    end
    
    MinIndices = cat(1,MinIndices_plus,MinIndices_minus);  %[Minimizers-Indices f_plus, Minimizer-Indices f_minus]
    
    for i = 1:(numMinSP_minus+numMinSP_plus-1)  % -1: dont count Minimizer itself.. 
        idx = MinIndices(i);   
        
        if i > numMinSP_plus
            SP = SP_minus(idx,:);
        else
            SP = SP_plus(idx,:); 
        end
        
%         SP_value = T_minus(idx) % not needed? 
        Matrix = sign*(SP'*SP);

        if norm(double(Matrix-Minimizer),'fro') < 1e-8   %check is this sufficient here?
%             fprintf('both StationaryPoints correspond to the same(Matrix-)Minimizer')
        else
%             fprintf('StationaryPoint corresponds to a different (Matrix-)Minimizer')
            MinimizerCount = MinimizerCount + 1;
        end
    end
% ----------------------------------------------------------------------------------------------------------------



% Output Uniqueness of Minimizers:
if print_Uniqueness
    if MinimizerCount == 1
        fprintf('Unique Minimzier')
    elseif MinimizerCount == 2
        fprintf('Two Minimziers')
    else
        fprintf('1-Parameter family of Minimziers')
    end
end





% --- determine the angle of the Minimizer
% a1 = Minimizer(1,1)
% a2 = Minimizer(2,2)
a1 = double(Minimizer(1,1));
a2 = double(Minimizer(2,2));

% compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
e = [sqrt((a1/(a1+a2))), sqrt((a2/(a1+a2)))];    % always positive under sqrt here .. basically takes absolute value here
angle = atan2(e(2), e(1)); 

% compute curvature kappa
kappa = (a1 + a2);

% % CHeck off diagonal entries:
% sqrt(a1*a2);
% double(Minimizer);

G = double(Minimizer);



% --- "Classification" / Determine the TYPE of Minimizer by using 
% the number of solutions (Uniqueness?)
% the angle (axial- or non-axial Minimizer)

% (Alternative compute det[GlobalMinimizer' e1']  where e1 = [1 0] ?)

% Check Uniqueness -- Options: unique/twoMinimizers/1-ParameterFamily
if MinimizerCount == 1
%     fprintf('Unique Minimzier')
    % Check if Minimizer is axial or non-axial:
    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
     Type = 3;
    else % non-axial Minimizer   
     Type = 1;
    end
elseif MinimizerCount == 2
%     fprintf('Two Minimziers')
    % Check if Minimizer is axial or non-axial:
    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
        Type = 3;
    else % non-axial Minimizer   
        fprintf('ERROR: Two non-axial Minimizers cannot happen!')
    end
else
%     fprintf('1-Parameter family of Minimziers')
    % Check if Minimizer is axial or non-axial:
    if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
%         fprintf('ERROR: axial Minimizers cannot happen for 1-Parameter Family!')
    else % non-axial Minimizer   
        Type = 2;
    end
end
% ------------------------------------------------------------------------------------------------------




if print_Output
    fprintf(' --------- Output symMinimization --------')
    fprintf('Global Minimizer v: (%d,%d) \n', GlobalMinimizer(1),GlobalMinimizer(2) )
    fprintf('Global Minimizer Value f(v): %d \n', sym(globalMinimizerValue) ) %cast to symbolic
    %     fprintf('Global Minimizer Value : %d', globalMinimizerValue )
    
    fprintf('Global Minimizer G: \n' )
    G 
    fprintf("Angle = %d \n", angle)
    fprintf("Curvature = %d \n", kappa)
    fprintf("Type = %i \n", Type)
    fprintf(' --------- -------------------- --------')
end



if make_FunctionPlot  
    fsurf(@(x,y) f_plus(x,y,q1,q2,q3,q12,b1,b2,b3)) % Plot functions
    hold on
    plot3(double(A_plus),double(B_plus),T_plus,'g*')
    %Plot GlobalMinimizer:
    hold on
    plot3(double(GlobalMinimizer(1)),double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
    % view(90,0)
    % view(2)

    figure
    fsurf(@(x,y) f_minus(x,y,q1,q2,q3,q12,b1,b2,b3))
    hold on
    plot3(double(A_minus),double(B_minus),T_minus,'g*')
    hold on
    plot3(double(GlobalMinimizer(1)), double(GlobalMinimizer(2)),globalMinimizerValue, 'o', 'Color','c')
end


return 


% Write symbolic solution to txt-File in Latex format
% fileID = fopen('txt.txt','w');
% fprintf(fileID,'%s' , latex(S_plus.v1));
% fclose(fileID);