function [newpiecedepths,BEFORE,AFTER] = corelogintegrate(GRA, log, pieces, percent)

% [NEWPIECEDEPTHS, BEFORE, AFTER] = corelogintegrate(GRA, LOG, PIECES, PERCENT)
% 
% This function takes GRA density data, compares it to log density data and
% readjusts the depths of the pieces of the core to minimize the density
% differences. The input variables are: GRA (piece density data), LOG
% (logging density data), PIECES (curated depths of the pieces), and
% PERCENT (decimal percent limit of density to restrict possible logging
% depths).
% -------------------------------------------------------------------------
% INPUTS: GRA, LOG, PIECES
%  GRA is a 9-column matrix of GRA density data from ODP/IODP core data, in
%  a format similar to Janus, http://www-odp.tamu.edu/database/.
%    GRA(:,1) is core number
%    GRA(:,2) is section number 
%    GRA(:,3) is depth to top of core (mbsf) 
%    GRA(:,4) is depth from top of section (cm)
%    GRA(:,5) is piece number 
%    GRA(:,6) is curated depth (mbsf) 
%    GRA(:,7) is density (g/cc) 
%    GRA(:,8) is length of core recovered from this section (m) 
%    GRA(:,9) is curated bottom of cored interval (mbsf)
%  LOG is a 2-column matrix of downhole logging data.
%    LOG(:,1) is depth (mbsf) LOG(:,2) is logging density (g/cc)  
%  PIECES is a 7-column matrix. 
%    PIECES(:,1) is core number 
%    PIECES(:,2) is section number 
%    PIECES(:,3) is piece number 
%    PIECES(:,4) is distance below top of section to top of piece (cm) 
%    PIECES(:,5) is distance below top of section to bottom of piece (cm) 
%    PIECES(:,6) is length of piece (cm) 
%    PIECES(:,7) is tide-corrected depth to top of core (mbsf)
%  PERCENT is user choice of percent match required between core and log
%  data, as a decimal (use 0.2 for 20% match), from the maximum GRA value
%  of a given piece.  Value used in Gilbert and Burke, 2008 is 0.2.
%
% OUTPUT: NEWPIECEDEPTHS,BEFORE,AFTER
%   NEWPIECEDEPTHS is a 4-column matrix. 
%    NEWPIECEDEPTHS(:,1) is core number
%    NEWPIECEDEPTHS(:,2) is section number 
%    NEWPIECEDEPTHS(:,3) is piece number 
%    NEWPIECEDEPTHS(:,4) is new depth to the top of the piece (mbsf)
%   BEFORE is a 4-column matrix representing the data before depth shifting
%    BEFORE(:,1) is curated depth (mbsf) 
%    BEFORE(:,2) is the density value for the core piece
%    BEFORE(:,3) is the nearest logging depth to the curated depth (mbsf)
%    BEFORE(:,4) is the density value at that nearest depth in the log record
%   AFTER is a 4-column matrix representing the data after depth shifting
%    AFTER(:,1) is new depth after shifting of physical property (mbsf) 
%    AFTER(:,2) is the density value for the core piece
%    AFTER(:,3) is the nearest logging depth to the new depth (mbsf)
%    AFTER(:,4) is the density value at that new nearest depth in the log record
%
% ------------------------------------------------------------------------
% *Variations*
%   Other data types (Vp, porosity, etc.) can be substituted as follows:
%   GRA(:,7) can be replaced with data obtained from the recovered core if
%   LOG(:,2) is similarly replaced with the same type of data from the log. 
%   All information in the remaining columns can remain unchanged.
% -------------------------------------------------------------------------
%
% Gilbert, L.A. and A. Burke, 2008, Depth-shifting cores
% incompletely recovered from the upper oceanic crust, IODP Hole 1256D,
% G-cubed, 2008gc002010.
%
% Revised 5/8/08, AB and LAG
% for further information contact aburke@mit.edu or lgilbert@williams.edu


finalmatches = [];
originalpieces = 0;
deleted4range = 0;
nologdata = 0;
multiple = 0;
% First want to isolate one core
j = 1;
while j <= length(GRA)
    matches= [];
    onecore = find(GRA(:,1) == GRA(j,1));
    % Now want to count the number of pieces in the core because we want to
    % know how many pieces we started with
    originalpieces = originalpieces + 1;
    for k = 1:length(onecore)-1
        if GRA(onecore(k),5) ~= GRA(onecore(k+1),5) | GRA(onecore(k),2) ~= GRA(onecore(k+1),2)
            originalpieces = originalpieces +1;
        end
    end
    % Now want to find the depths whose log value densities are w/in a certain percent of
    % the maximum GRA value of a given piece. 
    i= 1;
    while i <= length(onecore)
        pieceindex = find(GRA(onecore,2) == GRA(onecore(i),2) & GRA(onecore,5) == GRA(onecore(i),5)); %find same section and same piece
        piece = GRA(pieceindex + (j - 1), :);
        [C,I] = max(piece(:,7));
        maxdens = [piece(I, 6) piece(I,7) piece(I,1) piece(I,2) piece(I,5) piece(I,4)]; % curated depth, corrected density, core, section, piece, depth from top of section
        topbound = maxdens(1,1);
        bottombound = piece(I,9) - piece(I,8) + piece(I,6) - piece(I,3);
        velindex = find(log(:,1)>= topbound & log(:,1) < bottombound); % find possible log depths
        if isempty(velindex) % count the number of pieces that are ignored because there is no log data at the right depths
            nologdata = nologdata +1;
        end
        possiblematches = matchpercent(maxdens, log(velindex,:), percent);
        if size(possiblematches) >1 % count the number of pieces with multiple possible depths
            multiple = multiple +1;
        end
        matches = [matches; possiblematches];
        i = i + length(pieceindex);
    end
    %Get rid of all the depths that did not match up with any log values within
    %the restricted percentage range.
    badrows = find(matches(:,1) == -999.99);
    if not(isempty(badrows))
        deleted4range = deleted4range + length(badrows);
        matches(badrows,:) = [];
    end
    finalmatches = [finalmatches; matches];
    j = j+ length(onecore);
end

% Here we want to know the number of GRA data values (and hence number of
% pieces) left after all of this

remainingpieces = 1;
for p = 1:length(finalmatches)-1
    if finalmatches(p,3) ~= finalmatches(p+1,3)
        remainingpieces = remainingpieces +1;
    end
end

%Now pick out the best option for each piece by choosing the shallowest
%depth for each piece greater than the piece above it. If one doesn't exist
%than see if there is a depth greater than the depth of two pieces above
%it. If there is, then check which piece (the current piece or the prior
%piece) has a smaller difference between GRA density and log density, and
%choose that one

finaldiff = [finalmatches abs(finalmatches(:,2)-finalmatches(:,4))];

optimized = [];
deleted4order = 0;

j=1;
while j <= length(finaldiff);
    core = finaldiff(finaldiff(:,5) == finaldiff(j,5),:);%find a core
    k = 1;
    
    piecedeleted = 0;
    while k <= size(core,1)
        pieceindex = find(core(:,6) == core(k,6) & core(:,7) == core(k,7)); %find a piece
        if k == 1
            optimized = [optimized; core(k,:)]; % first piece is shallowest depth automatically
        end
        if k >1
            possible = find(core(pieceindex(:,1),1) >= optimized(end,1) + (core(k,3) - core(k-1,3))); %depths greater than previous piece
            a = size(optimized,1);
            if a > 1
                alternative = find(core(pieceindex(:,1),1) >= optimized(a-1,1) + (core(k,3)-optimized(a-1,3))); %depths greater than the pre-previous piece
            elseif a == 1
                alternative = find(core(pieceindex(:,1),1) >= optimized(a,1) + (core(k,3)-optimized(a,3)));
            end
            if not(isempty(possible)) % if there is a depth greater than previous piece
                shallowest = core(possible(1,1)+k-1,:); %take the shallowest one
                optimized = [optimized; shallowest];
            elseif not(isempty(alternative)) % if there is a depth greater than the pre-previous piece
                shallowest = core(alternative(1,1) + k - 1, :);
                if shallowest(1,9) < optimized(end,9) % check if the difference b/w piece and log densities is less than that of the previous piece
                    optimized(end,:) = shallowest; % if it is keep that depth for that piece and discard the previous piece depth
                end
                piecedeleted = piecedeleted + 1; %either way you have to delete one piece
            else % if there is no depth greater than either of the previous two pieces
                piecedeleted = piecedeleted +1; % you have to delete a piece
            end
        end
        
        k = k+length(pieceindex);
    end
    p = size(optimized,1);
    otheroptions = find(core(:,1) > optimized(p,1) & core(:,6) == optimized(p,6) & core(:,7) == optimized(p,7)); %find deeper possible depths for the last piece in the core
    if otheroptions % if deeper options exist
        bestoption = min(core(otheroptions,9)); 
        bestoptionindex = find(core(:,9) == bestoption & core(:,6) == optimized(p,6) & core(:,7) == optimized(p,7));
        optimized(p,:) = core(bestoptionindex(end,1),:); % choose the depth that minimizes the difference between piece and log density
    end
    p= p - 1;
    while p > 0 & optimized(p,5) == optimized(end,5) %same core
        otheroptions = find(core(:,1) > optimized(p,1) & core(:,1)< optimized(p+1,1) & core(:,6) == optimized(p,6) & core(:,7) == optimized(p,7)); %find other deeper depths for next piece from bottom of core
        if otheroptions %if deeper options exist
            bestoption = min(core(otheroptions,9));
            bestoptionindex = find(core(:,9) == bestoption & core(:,1) > optimized(p,1) & core(:,1)< optimized(p+1,1) & core(:,6) == optimized(p,6) & core(:,7) == optimized(p,7));
            optimized(p,:) = core(bestoptionindex(1,1),:); %choose the depth that minimizes the difference b/w piece and log density
        end
        p= p - 1;
    end
    
    deleted4order = deleted4order + piecedeleted; %total pieces deleted for all cores
    j = j + length(core);
end



% Make a matrix of the core, section, piece and new adjusted depth of the
% top of the piece
for q = 1:size(optimized,1)
    output(q,1) = optimized(q,5); %core
    output(q,2) = optimized(q,6); % section
    output(q,3) = optimized(q,7); %piece
    piecetop = pieces(pieces(:,1) == optimized(q,5) & pieces(:,2) == optimized(q,6) & pieces(:,3) == optimized(q,7), 4)/100; % meters from top of section of top of the piece
    output(q,4) = optimized(q,1) - (optimized(q,8)/100) + piecetop; %new adjusted top of the piece
end


% Now fit in the remaining pieces that have been deleted by placing them as
% shallow as they can be given the new depths of the pieces that weren't
% deleted
piecelog = pieces(pieces(:,1) <= output(end,1),:);
b=1;
newpiecedepths = [];
while b <= length(piecelog)
    core = piecelog(piecelog(:,1) == piecelog(b,1),:); %For each core
    % For the first piece
    tiedpiece = find(output(:,1) == core(1,1) & output(:,2) == core(1,2) & output(:,3) == core(1,3)); %See if it was assigned a new depth
    if tiedpiece % if it was assigned a new depth
        newpiecedepths = [newpiecedepths; output(tiedpiece,:)]; %keep the new depth
    else % if it was deleted
        newpiecedepths = [newpiecedepths; core(1,1) core(1,2) core(1,3) core(1,7)]; %assign its depth to the top of the core (the shallowest it can be)
    end
    % For pieces deeper than the first piece
    
    for c = 1:size(core,1)-1
        tiedpiece = find(output(:,1) == core(c+1,1) & output(:,2) == core(c+1,2) & output(:,3) == core(c+1,3)); %See if the piece was assigned a depth
        if tiedpiece %if it was assigned a new depth
            newpiecedepths = [newpiecedepths; output(tiedpiece,:)]; %keep the new depth
        else % if the piece was deleted
            prevpiecelength = core(c, 6);
            g = size(newpiecedepths,1);
            % assign it to the shallowest depth beneath the previous piece
            % (taking into consideration the length of the previous piece
            newpiecedepths = [newpiecedepths; core(c+1,1) core(c+1,2) core(c+1,3) newpiecedepths(g,4)+(prevpiecelength/100)];
        end
    end
    b=b+size(core,1);
end

%As a final step, calculate the new depths of all of the GRA data. The
%new matrix gradepths has density in first column, curated depth in 2nd
%column, and new depth in 3rd column.  If there is no new depth for the
%piece or if there is no record of the piece before shifting, then -999.99
%will be returned as the new depth

index1 = [];
index2 = [];
gradepths = [];

for i = 1:length(GRA)
    index1 = find(pieces(:,1) == GRA(i,1) & pieces(:,2) == GRA(i,2) & pieces(:,3) == GRA(i,5)); % find the index of the piece information before shiftiing
    index2 = find(newpiecedepths(:,1) == GRA(i,1) & newpiecedepths(:,2) == GRA(i,2) & newpiecedepths(:,3) == GRA(i,5)); %find the index of the piece information after shifting
    if size(index1) == 1 & size(index2)== 1 % if information about the piece exists
        gradepths = [gradepths; GRA(i,7) GRA(i,6) (newpiecedepths(index2,4)+(GRA(i,4)/100)-(pieces(index1,4)/100))]; %record the density, curated depth and new depth
    else
        gradepths = [gradepths; GRA(i,7) GRA(i,6) -999.99]; %otherwise return -999.99
        str = strcat('No piece information for core ', num2str(GRA(i,1)), 'section', num2str(GRA(i,2)) ,' piece ',num2str(GRA(i,5)));
        warning(str)
    end
end

%% And now find the closest log density before depth shifting and after
%% depth shifting for GRA

BEFORE = [];
AFTER = [];

for i = 1:length(gradepths)
    beforeindex = findnearest(gradepths(i,2),log(:,1));
    afterindex = findnearest(gradepths(i,3),log(:,1));
    BEFORE = [BEFORE; gradepths(i,2) gradepths(i,1) log(beforeindex(1,1),1) log(beforeindex(1,1),2)];
    AFTER = [AFTER; gradepths(i,3) gradepths(i,1) log(afterindex(1,1),1) log(afterindex(1,1),2)];
end