function [f1_f, f2_f, f3_f,OD_f, theta_f, phi_f] = nod(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_noddi(ydata);
[ydata1, sigma]=norm_meas_gen(ydata,b);
 %% Estimation using GA
 %GA Module used from http://www.acse.dept.shef.ac.uk/cgi-bin/gatbx-download
NIND = 48;           
MAXGEN = 60;         % maximum Number of ietrations
GGAP = 0.9;          % Generation gap, how many new individuals are created
NVAR = 4;            % Number of variables (OD theta phi v_ic)
PRECI = 15;          % Precision of binary representation

% Build field descriptor
 FieldD =[15 15  15  15  ;    0.01  0.01  0.01 0.1;...
     0.98    pi    pi  1;     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(@nod_eval_nlin,par(i,:),D_intra,D_iso,G_mag,Del,del,g,b,sigma,ydata1);
end 
ObjV=Objv';
   Best(gen+1) = min(ObjV);
   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(@nod_eval_nlin,par1(i,:),D_intra,D_iso,G_mag,Del,del,g,b,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,:);                % OD theta phi v_ic
%% Linear parameters fit using cvx
% toolbox available at http://cvxr.com/cvx/
[A1, A2, A3] = NOD_A(x1,D_intra,D_iso,G_mag,Del,del,g,b);
A=[A1 A2 A3];
 cvx_begin quiet
     variable f1e 
     variable f2e 
     variable f3e 
     fe=[f1e f2e f3e]';
         minimize(norm((ydata1)- A*fe)-sigma);   
     1.0e-3 <= f1e <= 0.9
     1.0e-3 <= f2e <= 0.9
     0 <= f3e <= 0.99
     f1e + f2e + f3e  == 1
  cvx_end
    v_ic=x1(4); v_iso=f3e;
    f1ga=(v_ic-(v_ic*v_iso));
    f2ga=(1-v_iso)*(1-v_ic);
    x2=[((f1e+f1ga)/2)  ((f2e+f2ga)/2)  f3e x1(1) x1(2) x1(3)];
%% Trust-region-reflective algorithm for non-linear data fitting
 opts = optimset('Display','off');
    lb=[1.0e-3 1.0e-3 0 0.01 0.01 0.01 ];
    ub=[0.9 0.9 0.99 0.98 pi pi ];
    [xf,resnorm,~,exitflag,output] = lsqcurvefit(@noddi_eval,x2,xdata,ydata1,lb,ub,opts);
f1_f=xf(1); f2_f= xf(2); f3_f=xf(3); OD_f=xf(4);theta_f=xf(5);  phi_f=xf(6);

