function [calmatrixv, btl0]=run_oxygen_cal_ml(sbecoef0,sbecoefd,gistat_cal,iswitch,wswitch,dirjpg,direps,dirmerge,cname,p_final,csuf,p_maxd,gdeprange,p_equation,p_plot);
%
% PLEASE UPDATE COEFFICIENTS FROM CALIBRATION SHEET in o2_params BEFORE HAND !!
%
% For oxygen fitting.  First run match_sgn to get the
% files needed.  gistat_cal is a vector of station numbers to fit.
%
% This fits up to 5 parameters: according to routine oxyg_eq and equation chosen (p_equation)
%
% If iswitch is set to:
% 0 the program fits the stations as a group
% 1 the program fits the stations as a group, but
% with a slope that is a linear function of station number.
% 2 the program first fits the stations as a group,
% then goes back and fits a slope and bias to individual
% stations, keepin other factors constant at the group values
%
% If wswitch is set to:
% 0 the fit is a regular L2 (least squares) norm.
% 1 the fit is L2, but it is weighted program weights  by
% the nominal oxygen bottle spacing, thus fitting the deep
% portion of the water column better.
% 2 the fit is an L1 norm (not least squares, sometimes
% called a robust fit, it tends to ignore outliers.  It
% might not work too well in fitting high and low values.
% 3 the fit is an L4 norm, which one may actually want to try
% to fit high and low values...
%
% deal with kludgy switches, anyting other than 1, 2 or 3 for
% wswitch is set to the default of zero...
% from G. Johnson, June 2008
% modifs A. Ganachaud, 2014

unix('rm TMP/*') % remove all temporary files

if wswitch~=1&wswitch~=2&wswitch~=3
    wswitch=0;
end % if

close all

gistat_cal=gistat_cal(:);

% load data from match_sgn

eval(['load ' dirmerge 'oxydata_' cname 'p'])
clear z % use btl. variable now


% first selects stations being fit: this will condition the output
% coefficent for a given station range (gistat_cal)
gis=find((ismember(btl.stat,gistat_cal)==1));

btl=keep_data(btl,gis);
% first and last stations; station linear function for iswitch=1 option
btl.stat_beg=min(btl.stat);
btl.stat_last=max(btl.stat);
dstas=btl.stat_last-btl.stat_beg+1;
%btl.dnstat contains linear "station" variable for all bottles (btl. vectors)
btl.dnstat=(btl.stat-btl.stat_beg+1)/dstas;

%now remove wrong bottles
if nargin>=11 & exist('p_maxd') %further filter by maximum depth to take into account
    %AG August 2015: added .*btl.wnklo2 as we may have kept all bottles including those without O2 samples
    gis=find(~isnan(btl.dsali.*btl.dtemp68.*btl.wnklo2ml) &  (btl.dpres<p_maxd) );
else
    gis=find(~isnan(btl.dsali.*btl.dtemp68.*btl.wnklo2ml) );
end

% ADD FLAG TO BTL
btl.wnklo2_flag=3*ones(length(btl.stat),1); %put all flags to 3 (bad measurements)
% A priori flag 2 (good measurement) for those bottles that will be
% processed here
btl.wnklo2_flag(gis)=2;
% Define btl0 as the original bottle data (btl will be shrunk as bad
% bottles are eliminated through iterations)
btl.gibtl=(1:length(btl.stat))'; % define the original sample index (to find our way in btl0 variable when we shrink btl vectors)
btl0=btl; %to have final bottle data with flags for outliers

%remove all data from btl but those of the selected stations from btl structure
btl=keep_data(btl,gis);
ndata=length(gis);

% save the file for use by oxygen_cal_ml.m during minimization
save TMP/temp_oxydata_ml btl wswitch gistat_cal sbecoefd p_equation

% initialize new coefficients and outlier matrix

sbecoef=sbecoef0;

if iswitch~=1&iswitch~=2
    iswitch=0;
end

% now deal with the various grouping options...

switch iswitch
    case 0 %simple fit
        p_slope=0;
        
        % iteratively loop throwing out 2.8 standard deviation outliers...
        noutliers=0;
        gioutlier=0;
        while length(gioutlier)>0
            
            sbecoef=fminsearch(@oxygen_cal_ml,sbecoef);
            
            btl.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl);
            
            % compute standard deviation
            
            btl.doxres=(btl.wnklo2um-btl.oxcal);
            stdres=std(btl.doxres);
            
            %figure(10);pl1=plot(btl.stat(:,1),btl.doxres,'ro');hold on;title('minimization progress...');ppause;set(pl1,'color','b')
            % throw out the outliers and flag them in btl0
            
            gioutlier=find(abs(btl.doxres)> 2.8*stdres);
            noutliers=noutliers+length(gioutlier);
            gis=find(abs(btl.doxres)<=2.8*stdres);
            %flag outliers for btl0
            btl0.wnklo2_flag(btl.gibtl(gioutlier))=3;
            
            %remove outliers from btl for iterations
            btl=keep_data(btl,gis);
            
            % save the file for use by oxygen_cal_ml.m
            save TMP/temp_oxydata_ml btl wswitch gistat_cal sbecoefd p_equation
            
        end % while
        
        % save the results in a file temp_oxydata
        percentused=100*(ndata-noutliers)/ndata;
        
        save TMP/temp_oxydata_ml btl btl0 sbecoef sbecoefd wswitch gistat_cal
        
        calmatrixv.sbecoefs=ones(size(gistat_cal))*sbecoef;
        calmatrixv.ndata=ones(size(gistat_cal))*ndata;
        calmatrixv.percentused=ones(size(gistat_cal))*percentused;
        calmatrixv.stdres=ones(size(gistat_cal))*stdres;
        %
        % re calculate oxygen from original data
        %
        btl0.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl0);
        
        %title for graphics
        ttlstr='Group Fit';
        btl0=diags_plots(btl0,ttlstr, gdeprange, p_final, csuf, iswitch, wswitch, direps, dirjpg);
        
    case 1 %iswitch
        %
        % Linear station dependance
        %
        p_slope=1;
        % start with a first guess of no station dependence...
        sbecoef(6)=0; %this is the linear dependence coefficient
        
        btl.dnstat=(btl.stat-btl.stat_beg+1)/dstas;
        %dnstat_outputmatrix is linear variable only for each station
        dnstat_outputmatrix=(gistat_cal-btl.stat_beg+1)/dstas;
        % ** this is for gistat_cal that will contain calib coeff for each
        % **station** not profile
        
        %need to save temp file again with dnstas variable
        save TMP/temp_oxydata_ml btl wswitch gistat_cal sbecoefd p_equation
        
        % iteratively loop throwing out 2.8 standard deviation outliers...
        noutliers=0;
        gioutlier=0;
        while length(gioutlier)>0
            
            sbecoef=fminsearch(@oxygen_cal_staslope_ml,sbecoef,  optimset('MaxFunEvals',400*length(sbecoef)));%'MaxIter',400*length(sbecoef)
            btl.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl);
            % compute standard deviation
            
            btl.doxres=(btl.wnklo2um-btl.oxcal);
            stdres=std(btl.doxres);
            %figure(10);pl1=plot(btl.stat(:,1),btl.doxres,'ro');hold on;title('minimization progress...');ppause;set(pl1,'color','b')
            
            % throw out the outliers and flag them in btl0
            gioutlier=find(abs(btl.doxres)> 2.8*stdres);
            noutliers=noutliers+length(gioutlier);
            gis=find(abs(btl.doxres)<=2.8*stdres);
            %flag outliers for btl0
            btl0.wnklo2_flag(btl.gibtl(gioutlier))=3;
            
            %remove outliers from btl for iterations
            btl=keep_data(btl,gis);
            
            % save the file for use by oxygen_cal_staslope_ml.m
            save TMP/temp_oxydata_ml btl wswitch gistat_cal sbecoefd p_equation
            
        end % while
        
        % save the results in a file temp_oxydata
        percentused=100*(ndata-noutliers)/ndata;
        
        save TMP/temp_oxydata_ml btl btl0 sbecoef sbecoefd wswitch gistat_cal
        
        calmatrixv.sbecoefs=[sbecoef(1)+sbecoef(6)*dnstat_outputmatrix,ones(size(gistat_cal))*sbecoef(2:5)];
        calmatrixv.ndata=ones(size(gistat_cal))*ndata;
        calmatrixv.percentused=ones(size(gistat_cal))*percentused;
        calmatrixv.stdres=ones(size(gistat_cal))*stdres;
        
        %
        % re calculate oxygen from original data
        %
        btl0.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl0);
        
        %title for graphics
        ttlstr='Group Fit with station dependent slope';
        btl0=diags_plots(btl0,ttlstr, gdeprange, p_final, csuf, iswitch, wswitch, direps, dirjpg);
        
    case 2 %iswitch
        % option iswitch=2: first fits the stations as a group,
        % then goes back and fits a slope and bias to individual
        % stations, keepin other factors constant at the group values
        
        % first fits the stations as a group
        p_slope=0;
        
        % iteratively loop throwing out 2.8 standard deviation outliers...
        noutliers=0;
        gioutlier=0;
        while length(gioutlier)>0
            
            sbecoef=fminsearch(@oxygen_cal_ml,sbecoef);
            
            btl.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl);
            
            % compute standard deviation
            
            btl.doxres=(btl.wnklo2um-btl.oxcal);
            stdres=std(btl.doxres);
            
            %figure(10);pl1=plot(btl.stat(:,1),btl.doxres,'ro');hold on;title('minimization progress...');ppause;set(pl1,'color','b')
            % throw out the outliers and flag them in btl0
            
            gioutlier=find(abs(btl.doxres)> 2.8*stdres);
            noutliers=noutliers+length(gioutlier);
            gis=find(abs(btl.doxres)<=2.8*stdres);
            %flag outliers for btl0
            btl0.wnklo2_flag(btl.gibtl(gioutlier))=3;
            
            %remove outliers from btl for iterations
            btl=keep_data(btl,gis);
            
            % save the file for use by oxygen_cal_ml.m
            save TMP/temp_oxydata_ml btl wswitch gistat_cal sbecoefd p_equation
            
        end % while
        
        % save the results in a file temp_oxydata
        percentused=100*(ndata-noutliers)/ndata;
        
        save TMP/temp_oxydata_ml btl btl0 sbecoef sbecoefd wswitch gistat_cal
        
        calmatrixv.sbecoefs=ones(size(gistat_cal))*sbecoef;
        calmatrixv.ndata=ones(size(gistat_cal))*ndata;
        calmatrixv.percentused=ones(size(gistat_cal))*percentused;
        calmatrixv.stdres=ones(size(gistat_cal))*stdres;
        
        %first diagnostic
        %
        % re calculate oxygen from original data
        %
        btl0.oxcal=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, btl0);
        %title for graphics
        ttlstr='Group Fit before station-by station slope/bias adj (iswitch =2)';
        btl0=diags_plots(btl0,ttlstr, gdeprange, 0);
        disp('END OF FIRST STEP: NOW PROCEED TO SLOPE AND BIAS OPTIMIZATION FOR EACH STATION INDIVIDUALLY')
        ppause
        
        % now go through station by station and adjust the slope and bias
        save TMP/oxystaadjust_ml btl sbecoef sbecoefd wswitch gistat_cal
        
        gslope=[];
        gbias=[];
        gndata=[];
        gpercentused=[];
        gsstdres=[];
        
        % loop through stations
        
        for igistat=1:length(gistat_cal);
            
            istat=gistat_cal(igistat);
            
            load TMP/oxystaadjust_ml
            slp_bias=sbecoef(1:2); %sb contains slope and bias
            gis=find(btl.stat==istat);
            btl=keep_data(btl,gis);
            save TMP/temp_oxystaadjust_ml btl sbecoef gistat_cal sbecoefd wswitch p_equation
            
            % accumulate number of data points for each station
            gndata=[gndata;length(gis)];
            
            % do the minimization
            
            % iterate on 2.8 standard deviations.
            gioutlier=0;
            while length(gioutlier)>0
                %somehow at each iteration we start with slp_bias from original
                %optimization for all station: we do not re-start with
                %slp_bias from latest step
                
                slp_bias=fminsearch(@slopeandbias_ml,slp_bias);
                
                sbecoefsp=sbecoef;
                sbecoefsp(1)=slp_bias(1);
                sbecoefsp(2)=slp_bias(2);
                btl.oxcal=oxyg_eq(p_equation, p_slope, sbecoefsp, sbecoefd, btl);
                
                % compute standard deviation
                btl.doxres=(btl.wnklo2um-btl.oxcal);
                stdres=std(btl.doxres);
                
                %flag outliers
                gioutlier=find(abs(btl.doxres)> 2.8*stdres);
                noutliers=noutliers+length(gioutlier);
                gis=find(abs(btl.doxres)<=2.8*stdres);
                %flag outliers for btl0
                btl0.wnklo2_flag(btl.gibtl(gioutlier))=3;
                
                %remove outliers from btl for iterations
                btl=keep_data(btl,gis);
                
                save TMP/temp_oxystaadjust_ml btl sbecoef gistat_cal sbecoefd wswitch p_equation
                
            end % while
            % save calibrated oxygen in btl0
            btl0.oxcal(btl.gibtl)=btl.oxcal;
            
            % accumulate residual for each station
            gsstdres=[gsstdres;stdres];
            
            % accumulate percent used for each station
            gpercentused=[gpercentused;100*(gndata(igistat)-noutliers)/gndata(igistat)];
            
            % accumulate slope and bias for each station
            gslope=[gslope;slp_bias(1)];
            gbias=[gbias;slp_bias(2)];
            
            %title for graphics
            ttlstr=sprintf('Individual Slp-Bias fit Station %d ',istat);
            sbecoefsb=sbecoef;
            sbecoefsb(1)=slp_bias(1);
            sbecoefsb(2)=slp_bias(2);
            
            %good bottles
            gigoodbtl=find(btl.wnklo2_flag==2);
            %bad bottles
            gibadbtl=find(btl.wnklo2_flag==3);
            
            if p_plot %intermediate step for each station
                figure(5); clf; orient landscape;
                q=get(gcf,'pos'); q(1)=100; q(2)=100; set(gcf,'pos',q)
                subplot(1,2,1);
                plot(btl.oxcal,-btl.dpres,'kx',btl.wnklo2um(gigoodbtl),-btl.dpres(gigoodbtl),'bo',btl.wnklo2um(gibadbtl),-btl.dpres(gibadbtl),'ro')
                title(ttlstr)
                subplot(1,2,2);
                plot(btl.doxres(gigoodbtl),-btl.dpres(gigoodbtl),'bo',btl.doxres(gibadbtl),-btl.dpres(gibadbtl),'ro');
                set(gca,'xlim',[-16 16],'xtick',[-16:4:16])
                grid on
                ppause
            end
            
        end % for igistat=1:length(gistat_cal)
        
        save TMP/temp_oxystaadjust_ml btl0 sbecoef sbecoefd wswitch gistat_cal  p_equation
        %calmatrix=[gistat_cal,gslope,gbias,ones(size(gistat_cal))*sbecoef(3:5),gndata,gpercentused,gsstdres,ones(size(gistat_cal))*iswitch,ones(size(gistat_cal))*wswitch];
        calmatrixv.sbecoefs=[gslope,gbias,ones(size(gistat_cal))*sbecoef(3:5)];
        calmatrixv.ndata=gndata;
        calmatrixv.percentused=gpercentused;
        calmatrixv.stdres=gsstdres;
        
        ttlstr='After individual station fit';
        btl0=diags_plots(btl0,ttlstr, gdeprange, p_final, csuf, iswitch, wswitch, direps, dirjpg);
        
end %if wswitch~=1&iswitch~=2
% calibration matrix: contains coefficients/information for all stations
calmatrixv.stat=gistat_cal;
calmatrixv.iswitch=ones(size(gistat_cal))*iswitch;
calmatrixv.wswitch=ones(size(gistat_cal))*wswitch;

end %function run_oxygen_cal_ml

function btl0=diags_plots(btl0,ttlstr, gdeprange, p_final, csuf, iswitch, wswitch, direps, dirjpg)
%calibration diagnostics and graphics

%good bottles
gigoodbtl=find(btl0.wnklo2_flag==2);
%bad bottles
gibadbtl=find(btl0.wnklo2_flag==3);
btl0.doxres=btl0.oxcal-btl0.wnklo2um;
%
%Global statistics
%
disp('**********************************************************************')
disp(' ')
disp(sprintf('Global residual is %5.3f std dev is %5.3f; rms is %5.3f',...
    mean(btl0.doxres(gigoodbtl)), std(btl0.doxres(gigoodbtl)), rms(btl0.doxres(gigoodbtl))))
disp(' ')
disp('**********************************************************************')

%
% now do some diagnostic plots
%

figure(1); orient landscape; %wysiwyg;
q=get(gcf,'pos'); q(1)=100; q(2)=100; set(gcf,'pos',q)
subplot(1,2,1)
plot(btl0.stat(gigoodbtl),btl0.doxres(gigoodbtl),'+',btl0.stat(gibadbtl),btl0.doxres(gibadbtl),'ro');
set(gca,'xlim',[min(btl0.stat)-0.5 max(btl0.stat)+0.5],'ylim',[-16 16],'ytick',[-16:4:16])
grid on
xlabel('Station number');
ylabel('CTD - Bottle Oxygen [\mu mol kg^{-1}]')
legend('keepers','outliers');
title(sprintf('Global residual is %5.3f rms', rms(btl0.doxres(gigoodbtl))))

subplot(1,2,2)
plot(btl0.doxres(gigoodbtl),-btl0.dpres(gigoodbtl),'+',btl0.doxres(gibadbtl),-btl0.dpres(gibadbtl),'ro');
set(gca,'xlim',[-16 16],'xtick',[-16:4:16])
grid on
ylabel('Pressure [dbar]');
xlabel('CTD - Bottle Oxygen [\mu mol kg^{-1}]');
title(ttlstr);

%for each depth range
figure(4);clf;orient portrait; wysiwyg; setlargefig
for idepr=1:length(gdeprange)-1
    subplot(length(gdeprange)-1,1,idepr)
    gidep=find(btl0.dpres(gigoodbtl)>gdeprange(idepr) & btl0.dpres(gigoodbtl)<=gdeprange(idepr+1));
    gidepbad=find(btl0.dpres(gibadbtl)>gdeprange(idepr) & btl0.dpres(gibadbtl)<=gdeprange(idepr+1));
    plot(btl0.stat(gigoodbtl(gidep)),btl0.doxres(gigoodbtl(gidep)),'+',btl0.stat(gibadbtl(gidepbad)),btl0.doxres(gibadbtl(gidepbad)),'ro');
    set(gca,'xlim',[min(btl0.stat)-0.5 max(btl0.stat)+0.5],'ylim',[-6 6],'ytick',[-6:2:6])
    hold on
    plot([min(btl0.stat)-0.5 max(btl0.stat)+0.5],[0 0],'k')
    grid on
    ylabel(sprintf('%4d m\n%4d m',gdeprange(idepr),gdeprange(idepr+1)))
end
xlabel('CTD-BTL; Station');
legend('keepers','outliers');

% Rajout glazeo - plot by bottle_id
figure(10)
orient landscape; %wysiwyg;
q=get(gcf,'pos'); q(1)=100; q(2)=100; set(gcf,'pos',q) ;
plot(btl0.btlid(gigoodbtl),btl0.doxres(gigoodbtl),'+',btl0.btlid(gibadbtl),btl0.doxres(gibadbtl),'ro');
set(gca,'xlim',[min(btl0.btlid)-0.5 max(btl0.btlid)+0.5],'ylim',[-16 16],'ytick',[-16:4:16])
grid on
xlabel('Bottle Id');
ylabel('CTD - Bottle Oxygen [\mu mol kg^{-1}]');
title(sprintf('Global residual is %5.3f rms', rms(btl0.doxres(gigoodbtl))))

% Rajout A.Ganachaud - plot by flask_id
figure(11)
orient landscape; %wysiwyg;
q=get(gcf,'pos'); q(1)=100; q(2)=100; set(gcf,'pos',q) ;
plot(btl0.flask(gigoodbtl),btl0.doxres(gigoodbtl),'b+',btl0.flask(gibadbtl),btl0.doxres(gibadbtl),'ro');
set(gca,'xlim',[min(btl0.flask)-0.5 max(btl0.flask)+0.5],'ylim',[-10 10],'ytick',[-10:2:10])
grid on
xlabel('Flask ID');
ylabel('CTD - Bottle Oxygen [\mu mol kg^{-1}]');
title(sprintf('Global residual is %5.3f rms', rms(btl0.doxres(gigoodbtl))))

if p_final
    
    namefig=sprintf('o2cal_residuals_stat_%s_swslp_%d_%d',csuf,iswitch,wswitch)
    set(1,'paperunits','cent','paperorient','port')
    set(1,'position',[561   460   875   496])
    nfopeps=sprintf('%s%s.eps',direps,namefig)
    exportfig(1,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);
    
    nfopjpg=sprintf('%s%s.png',dirjpg,namefig)
    exportfig(1,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);
    
    % rajout glazeo
    set(4,'paperunits','cent','paperorient','port')
    % fin rajout glazeo-21/01/2015
    set(4,'position',[929   150   639   857])
    nfopeps=sprintf('%s%s-bydepthrange.eps',direps,namefig)
    exportfig(4,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
    nfopjpg=sprintf('%s%s-bydepthrange.png',dirjpg,namefig)
    exportfig(4,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',25);
    
    % rajout glazeo - plot by bottle_id
    namefig=sprintf('o2cal_residuals_stat_%s_swslp_%d_%d',csuf,iswitch,wswitch)
    set(10,'paperunits','cent','paperorient','port')
    set(10,'position',[561   460   875   496])
    nfopeps=sprintf('%s%s-bybottleid.eps',direps,namefig)
    exportfig(10,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);
    
    nfopjpg=sprintf('%s%s-bybottleid.png',dirjpg,namefig)
    exportfig(10,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);
    
    % rajout A. Ganachaud - plot by flask_id
    namefig=sprintf('o2cal_residuals_stat_%s_swslp_%d_%d',csuf,iswitch,wswitch)
    set(11,'paperunits','cent','paperorient','port')
    set(11,'position',[561   460   875   496])
    nfopeps=sprintf('%s%s-byflaskid.eps',direps,namefig)
    exportfig(11,nfopeps,'preview','tiff','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);
    
    nfopjpg=sprintf('%s%s-byflaskid.png',dirjpg,namefig)
    exportfig(11,nfopjpg,'format','png','LineMode','scaled','color','rgb',...
        'FontMode','fixed','FontSize',10,'height',12);


end
end %function diags_plots
