%--------------------------------------------------------------------------------------------------%
%
% DriverOptimization.m: 
%   Main program for optimizing density coefficents or
%   density coefficients and model parameters.
%
%   Parameters:
%       none
%
%   Returns: 
%       none (saves optimization results to .mat file)
%
%   Required files: load_initial_parameters.m, load_experimental_data.m,
%       initialize_optimization_allparams.m,
%       initialize_optimization_densonly.m, costfunctionSS.m,
%       costfunctionRes.m, costfunctionDensOnly.m, subinterval_pts.m,
%       extract_densitities.m, AlphaLogNormal.m, AlphaNormal.m, pxx_he.m,
%       sma_strain.m, pxx_he_partials.m, pxx_he_Deeppartials.m,
%       resistance.m, resistance_partials.m
%                               
%                       ----------------------------
%
% Copyright (C) <2011> by <Ralph Smith (rsmith@ncsu.edu), John Crews(jhcrews@ncsu.edu)>
% <Department of Mathematics, North Carolina State University, Raleigh, NC 27695>
%
% Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
% associated documentation files (the "Software"), to deal in the Software without restriction,
% including without limitation the rights to use, copy, modify, merge, publish, distribute,
% sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all copies or
% substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
% NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
% NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT
% OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%--------------------------------------------------------------------------------------------------%


clear all;
close all;
clc

params = load_initial_parameters();
load_experimental_data;

% controls optimization
% params.optimizeAll = 1: optimizes densities and model parameters
% params.optimizeAll = 0: optimizes densities only
params.optimizeAll = 1;

if params.optimizeAll == 1
    initialize_optimization_allparams;
    app = [app '_allparams'];
    
    % optimize stress-strain parameters
    fnc_ptrSS = @(x)costfunctionSS(x,exp_t,ic,exp_stress,exp_temp,params,xScaleSS,exp_strain);
    options = optimset('Jacobian','on','MaxIter',1000,'MaxFunEvals',5000,'TolFun',1e-6,'TolX',1e-8,'Display','iter');
    [x_optSS,resnorm_optSS,residual_optSS,exitflag_optSS,output_optSS,lambda_optSS,jacobian_optSS] = lsqnonlin(fnc_ptrSS,x0,LB,UB,options);
    
    % extract previously optimized stress-strain parameters
    weights = x_optSS(1:length(params.sigmaRmean)+length(params.b1));
    params.alpha = weights(1:length(params.sigmaRmean));
    params.beta = weights(length(params.sigmaRmean)+1:length(params.sigmaRmean)+length(params.b1));
    params.alpha = params.alpha/sum(params.alpha);
    params.beta = params.beta/sum(params.beta);
    params.Ea = x_optSS(nAlpha+nBeta+1)*xScaleSS(1);
    params.Em = x_optSS(nAlpha+nBeta+2)*xScaleSS(2);
    params.sigmaL = x_optSS(nAlpha+nBeta+3)*xScaleSS(3);
    params.deltaSigmaTemp = x_optSS(nAlpha+nBeta+4)*xScaleSS(4);
    params.epsT = x_optSS(nAlpha+nBeta+5)*xScaleSS(5);
    params.tau = x_optSS(nAlpha+nBeta+6)*xScaleSS(6);
    params.V = x_optSS(nAlpha+nBeta+7)*xScaleSS(7);
    % extract densities
    params = extract_densities(params);
    
    % get phase fractions for resistance calculations
    [xp xm xa strain res] = sma_hem_constanttemp(exp_t{1},ic{1},exp_stress{1},exp_temp{1},params);
    phase_fractions{1}.xp = xp;
    phase_fractions{1}.xm = xm;
    phase_fractions{1}.xa = xa;
    [xp xm xa strain res] = sma_hem_constanttemp(exp_t{2},ic{2},exp_stress{2},exp_temp{2},params);
    phase_fractions{2}.xp = xp;
    phase_fractions{2}.xm = xm;
    phase_fractions{2}.xa = xa;
    
    % optimize resistance parameters
    fnc_ptrRes = @(x)costfunctionRes(x,exp_stress,exp_temp,params,xScaleRes,exp_res,phase_fractions);
    options = optimset('Jacobian','on','MaxIter',800,'MaxFunEvals',800,'TolFun',1e-6,'TolX',1e-8,'Display','iter');
    [x_optRes,resnorm_optRes,residual_optRes,exitflag_optRes,output_optRes,lambda_optRes,jacobian_optRes] = lsqnonlin(fnc_ptrRes,x0Res,LBRes,UBRes,options);
    
    % extract parameters
    params.rhoA0 = x_optRes(1)*xScaleRes(1);
    params.rhoM0 = x_optRes(2)*xScaleRes(2);
    params.alphaA = x_optRes(3)*xScaleRes(3);
    params.alphaM = x_optRes(4)*xScaleRes(4);
    params.nu = x_optRes(5)*xScaleRes(5);
else
    initialize_optimization_densonly;
    app = [app '_densonly'];
    
    % optimize densitiy coefficients only
    fnc_ptr = @(x)costfunctionDensOnly(x,exp_t,ic,exp_stress,exp_temp,params,exp_strain);
    options = optimset('Jacobian','on','MaxIter',1000,'MaxFunEvals',5000,'TolFun',1e-8,'TolX',1e-8,'Display','iter');
    [x_opt,resnorm_opt,residual_opt,exitflag_opt,output_opt,lambda_opt,jacobian_opt] = lsqnonlin(fnc_ptr,x0,LB,UB,options);
    
    % extract parameters
    params.alpha = x_opt(1:length(params.c1));
    params.beta = x_opt(length(params.c1)+1:length(params.c1)+length(params.b1));
    % extract densities
    params = extract_densities(params);
end

save(['OptimizationResults' app '.mat'],'params');