function [f1_f, f2_f, f3_f,f4_f,theta_f, phi_f, R_f] = Aax_exvivo(xdata,Gamma,D_intra,D_iso,ydata)
%% Data Shaping
g=[xdata(:,1) xdata(:,2) xdata(:,3)];
G_mag= (xdata(:,4)./1000000); %(Tesla/micro meter)
Del=xdata(:,5);          %(in seconds)
del=xdata(:,6);          %(in seconds)
td=(Del-(del/3));        %(in seconds)
b=(td).*((Gamma.*del.*G_mag).^2); %(sec/mircometer^2)
ydata1=norm_meas_Aax(ydata);
sigma=0.05;  % for SNR 20 of ex-vivo data
%% Estimation of non-linear parameters using GA 
%GA module used from http://www.acse.dept.shef.ac.uk/cgi-bin/gatbx-download
NIND = 48;         
MAXGEN = 80;        % maximum Number of iterations
GGAP = 0.9;         % Generation gap, how many new individuals are created
NVAR = 4;           % Number of variables (i.e., theta, phi, R, est_f)
PRECI = 15;         % Precision of binary representation

% Build field descriptor
  FieldD =[ 15    15    15  15;...
    0.01  0.01   0.1   0.1; pi     pi    11    0.9;...
    1     1     1  1;  0   0     0     0;  1   1     1     1;  1   1     1     1];
% Initialise population
   Chrom = crtbp(NIND, NVAR*PRECI);

% Reset counters
   Best = NaN*ones(MAXGEN,1);	% best in current population
   gen = 0;			            % generational counter

% Evaluate initial population
par=bs2rv(Chrom,FieldD);
for i=1:NIND
Objv(i)=feval(@Aax_exvivo_nlin,par(i,:),g,Gamma,G_mag,del,Del,b,D_intra,D_iso,sigma,ydata1); 
end 
ObjV=Objv';
   Best(gen+1) = min(ObjV);
% Generational loop
   while gen < MAXGEN,

    % Assign fitness-value to entire population
       FitnV = ranking(ObjV);

    % Select individuals for breeding
       SelCh = select('sus', Chrom, FitnV, GGAP);

    % Recombine selected individuals (crossover)
       SelCh = recombin('xovsp',SelCh,0.7);

    % Perform mutation on offspring
       SelCh = mut(SelCh);

    % Evaluate offspring, call objective function
    par1=bs2rv(SelCh,FieldD);
    for i=1:max(size(par1))
    ObjVsel(i)=feval(@Aax_exvivo_nlin,par1(i,:),g,Gamma,G_mag,del,Del,b,D_intra,D_iso,sigma,ydata1); 
    end    
    ObjVSel = ObjVsel';  
    % Reinsert offspring into current population
       [Chrom ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel);

    % Increment generational counter
       gen = gen+1;

    % Update display and record current best individual
       Best(gen+1) = min(ObjV);
   end 
[ObjVSel,ind]=sort(ObjVSel);
par1=par1(ind,:);
x1=par1(1,:);
%% Linear parameters fit using cvx
% toolbox available at http://cvxr.com/cvx/
[A1, A2, A3, A4] = Aax_exvivo_A(x1,g,Gamma,G_mag,del,Del,b,D_intra,D_iso);
Aa=[A1 A2 A3 A4];
 cvx_begin quiet
     variable f1e 
     variable f2e 
     variable f3e 
     variable f4e 
     fe=[f1e f2e f3e f4e]';
         minimize(norm((ydata1)- (Aa*fe)-sigma));   % for offset gaussian
        %  minimize(norm((ydata1)- (Aa*fe)))   
     0 <= f1e <= 1
     0 <= f2e <= 1
     0 <= f3e <= 1
     0 <= f4e <= 1
     f1e + f2e + f3e + f4e == 1
     f1e <= x1(4)
     cvx_end
  x2=[f1e  f2e  f3e x1(1) x1(2) x1(3) f4e];
%% Trust-region-reflective algorithm for non-linear data fitting
 opts = optimset('Display','off');
  lb=[0.01 0.01  0.01 0.01 0.01 0.1 0.01];
  ub=[0.9  0.9  0.9   pi pi 11  0.9];
  [xf,resnorm,~,exitflag,output] = lsqcurvefit(@Aax_exvivo_eval,x2,xdata,ydata1,lb,ub,opts);
f1_f=xf(1); f2_f= xf(2); f3_f=xf(3); f4_f=xf(7); theta_f=xf(4); phi_f=xf(5); R_f=xf(6);

