function [corrected, details] = pin(year,d13c,alpha,varargin)

% PIN - PreINdustrial Correction for Carbon Isotope Tree-Ring Series for pC02
% 
%     [corrected, [details]] = pin(year,c,span)
%
%  A user-written Matlab function that performs a correction on carbon
%  isotope series from tree-rings to account for changes in atmospheric
%  pCO2 during the last ~150 years.
% 
%  The function is based on the physiological mechanisms and logical
%  constraints discussed in detail in McCarroll et al. [2009].  
%
%
% Inputs:   year = vector (m x 1) of years (required)
%           d13c    = vector (m x n) of atmospheric-corrected d13C values
%                  [there can be missing values (designated as 'NaN')]
%                  (required)  
%           alpha = [optional, default = 0.65]. This is the fraction of the
%                   time series used in the local (LOESS) regression (see
%                   Cleveland 1993). Must be between 0 and 1.
%
%           'verbose' = will print diagnostic figures and save results for debugging
%
%           'noplot' = do not create the final plot (batch mode)
%
% Outputs:  corrected = vector of PIN corrected d13C values
%           details   = optional structure for debugging and supporting information
%

% Revision History
% v0.1 - Original coding of loess regression
% v0.2 - full code for pin correction from initial d13C data
% v0.3 - heavy debugging to match R/Excel output and diagnose differences
% v0.4 - modified to accept multiple carbon isotope series
% v0.5 - minor corrections, extend ca data to 2005
% v0.9 - peer-review version 
% v0.91 - changes based on reviews and new testing
% v0.92 - heavy changes based on peer review
% v0.93 - heavy debugging to match R package
% v1.00 - final version for publication
%
% References: 
%
% McCarroll, D., Gagen, M. H., Loader, N. J., Robertson, I., Anchukaitis, K. J.,
% Los, S., Young, G. H. F., Jalkanen, R., Kirchhefer, A., and Waterhouse, J. S.
% (2009). Correction of tree ring stable carbon isotope chronologies for changes
% in the carbon dioxide content of the atmosphere. Geochimica et Cosmochimica 
% Acta, 73(6):1539?1547.
%
% Cleveland, W.S. and S. J. Devlin, 1988, Locally-Weighted Regression:
% An Approach to Regression Analysis by Local Fitting. Journal of the 
% American Statistical Association, 83:596-610.
%
% Cleveland, W.S. and E. Grosse, 1991, Computational Methods for Local
% Regression. Statistics and Computing, 1:47-62, 1991
%
% Cleveland, W.S. 1993, Visualizing Data, Hobart Press.

% Original Code Author: Kevin J Anchukaitis (kanchuka@ltrr.arizona.edu)
% Loess smoothing based on data visualization toolbox code by Cleveland 1993

%% Screen display
disp([' '])
disp(['PreINdustrial Correction for d13C tree-ring time series [PIN]'])
disp(['-------------------------------------------------------------'])
disp(['Version 1.00  March 2009   McCarroll et al. 2009'])
disp([' '])
disp(['McCarroll et al. 2009, Correction of tree ring stable carbon'])
disp(['isotope chronologies for changes in the carbon dioxide'])
disp(['content of the atmosphere, Geochim. Cosmochim. Acta, 73(6):1539?1547'])
disp([' '])

%% Check for verbose mode
if nargin > 3;
 if strcmp(varargin{1}, 'verbose')
  verboseMode = 1; % Extra graphics and output for debugging
 else
  verboseMode = 0; % No verbose debugging
 end
end

%% Check for plotting mode
if nargin > 3;
 if strcmp(varargin{1}, 'noplot')
  plotMode = 0; % Final graphic file
 else
  plotMode = 1; % no final graphic file
 end
end

%% Error Checking 
% Is the number of years equal to the number of isotope records?
if length(d13c) ~= length(year)
    error('Length of year vector and length of carbon isotope data must be the same')
end

% get the span value, or set to the default span 
% some intermingling here of of 'span' and 'alpha' is for consistency with R
if nargin < 3; alpha = 0.65; else alpha = alpha; end % a relatively wide 'local' span (alpha)
if any(size(alpha))~=1 | alpha < 0 | alpha > 1; 
  error('Alpha should be a single number between 0 and 1')
end

%% STEP 1. Begin Correction Procedure
% Sort seamlessly such that the data is by ascending year and arrange data vectors
data = sortrows([year d13c],1); year = data(:,1); d13c = data(:,2:end);
disp('Input data are automatically reordered prior to correction if necessary')

% save the original data before we start to operate on it
original.d13c.year = year; original.d13c.ratio = d13c;
d13co = d13c;

% How many individual series are there (nseries)?
[mseries,nseries] = size(d13c);

% Preallocate the output structure for speed and error catching
corrected.d13c.year = year;  corrected.d13c.ratio = NaN * ones(mseries,nseries);

% loop over for the number of individual series
for iter = 1:nseries

 graphname = ['pinDiagnostics',num2str(iter),'.ps'];
 disp([' '])
 disp(['Correcting carbon isotope series ', num2str(iter,0), ' ... '])
 full = ~isnan(d13co(:,iter)); % find the continous series
 yyear = year(full); % this is span of years for the series we'll work on
 cc = d13co(full,iter); % this is the series we'll work on

 
%% STEP 1 Retrieve Annual Values of ca 
if min(yyear) >= 1850  % do not operate on preindustrial values, result sensitive to this date
 firstyear = min(yyear);
else
 firstyear = 1850;
end ; 

if max(yyear) > 2005;
 lastyear = 2005;
else
 lastyear = max(yyear);
end

ca = co2(firstyear,lastyear); 
ic = find(yyear >= firstyear & yyear <= lastyear);
iyear = yyear(ic);
c = cc(ic);  % c and iyear go together now

% STEP 2 convert d13c to smoothed d13c using loess
[m,n] = size(iyear);

 % ... using loess regression (Cleveland and Grosse 1991)
 % Create the loess weights
 span = alpha*(lastyear-firstyear)/(2002-1850); % modify alpha but should be close to input 
 q = floor(span*m);  % alpha expressed as percentage of series length
 q = max(q,1); q = min(q,m);
 
 datails(iter).alpha = alpha;
 details(iter).span = span; % store the span size
 details(iter).spanN = q; % store the sample size

 for i = 1:m ; % for every year ...
   delta = abs(iyear(i) - iyear(:)); % get the distance from current point ...
   dsorted = sort(delta); % ... sort the values from lowest to highest (closest to furthest) ...
   dmax = dsorted(q); % ... find the maximum value within the desired span ...
   w = (1-abs(delta./(dmax)).^3).^3; % ... and calculate the tricubic weights.
   
  % index into the span of interest only
  gtz = find(w > 0); % find only the weights > 0, which will limit the least squares regression to the alpha proportion    
  xc = iyear(gtz); yc = c(gtz); wc = w(gtz); % cull the span from the data
  WC = spdiags(wc,0,length(wc),length(wc));
  
  % Construct the Vandermonde matrix for a second order (degree 2) polynomial
  V = [xc.^2 xc ones(length(xc),1)]; % simpler than built in function vender() for this operation 
  
  % QR decomposition for polynomial coefficients
  warning off;
  [Q,R] = qr((V'*WC*V),0); % qr is a function in both Matlab and Octave
  p = R\(Q'*(V'*WC*yc));  % this is most accurate result
    
  % use polyval.m (in both Matlab and Octave) to get predicted ci values from ca
  yhat(i) = polyval(p,iyear(i));
  cLow(i) = yhat(i);
  pp(i,:) = p;
 end


d13c.low = cLow(:); 
d13c.high = c - cLow(:);

warning on
 
if verboseMode == 1
    figure(1); clf;
    plot(iyear,c,'k'); ylabel('d13c'); xlabel('year'); hold on;
    plot(iyear,d13c.low,'b','linewidth',1.5);
    title('high and low frequency d13C')
    legend('Observed d13C','Low Frequency d13C','location','best')
    legend boxoff
    print('-dpsc',graphname)
end

clear m n p


%% STEP 3 Convert corrected low frequency d13C into ci
ci.actual = ca.ratio .* ((-10.8 - d13c.low)/22.6);

%% STEP 4 calculate ci for constant ci/ca ratio
ca.initial = ca.ratio(1); % set the initial ca
ci.initial = ci.actual(1); % set the initial ci to the actual ci
d13c.initial = d13c.low(1); 
ci.rcica.constant = ca.ratio * (ci.initial/ca.initial); % ci expected for constant ratio (contraint 1)

%% STEP 5 calculate ci for constant difference between ci - ca 
ci.dcica.constant = ca.ratio + (ci.initial - ca.initial); % ci expected for constant difference (constrait 2)

%% STEP 6 convert ci for constant ci/ca ratio and constant ci-ca difference to d13C units
d13c.rcica.ts = -10.8 - (22.6.*(ci.rcica.constant./ca.ratio)); % this should to be (nearly) constant
d13c.dcica.ts = -10.8 - (22.6.*(ci.dcica.constant./ca.ratio));

%% STEP 7 plot low frequency d13 c and constraints
if verboseMode
 figure(2); clf;
 plot(iyear,d13c.low,'k','linewidth',1.5); hold on; % ylim=range(c(d13c.low, d13c.rcica, d13c.dcica)) , type="l" , xlab="Year", ylab="d13C")
 plot(iyear,d13c.rcica.ts,'y','linewidth',1.5)
 plot(iyear,d13c.dcica.ts,'b','linewidth',1.5)
 legend('low frequency d13C','Constant Ratio Constraint','Constant Difference Contraint','location','best')
 legend boxoff
 grid on
 title('Comparison of low frequency d13c and constraints')
 print('-dpsc','-append',graphname)
end

%% STEP 8 Calculate increments for low frequency d13c (d13c.low) and for its 2 constraints
inc.d13c.low = [0; diff(d13c.low)];
inc.d13c.rcica = [0; diff(d13c.rcica.ts)];
inc.d13c.dcica = [0; diff(d13c.dcica.ts)];

% Preallocated final increment and set to zero
fin.inc1 = zeros(size(inc.d13c.low)); 
fin.inc2 = zeros(size(inc.d13c.low));

%% STEP 9 compare d13c increments against logical contraints; retain d13c increments only if within contraints
fin.inc1(find(inc.d13c.low > 0)) = inc.d13c.low(find(inc.d13c.low > 0)); % Constraint 1
fin.inc2(find((inc.d13c.low - inc.d13c.dcica) < 0)) = inc.d13c.low(find((inc.d13c.low - inc.d13c.dcica) < 0)) - inc.d13c.dcica(find((inc.d13c.low - inc.d13c.dcica) < 0)); % Constraint 2

%% STEP 10 add d13c increments after contraints are applied
fin.inc = fin.inc1 + fin.inc2;

if verboseMode
    figure(3); clf
    plot(iyear,inc.d13c.low,'k--','linewidth',2); hold on;
    plot(iyear,inc.d13c.rcica,'y--','linewidth',2);
    plot(iyear,inc.d13c.dcica,'b--','linewidth',2);
end

%% STEP 11 calculate cumulative correction and add to low frequency d13c; implement no change prior to 1850
fin.inc(find(iyear<1850)) = 0;

if verboseMode
  figure(3); 
  plot(iyear,fin.inc,'color',[0.75 0.75 0.75]);
  grid on
  legend('increment d13C low','increment d13c constant ratio constraint','increment d13c constant difference restraint','increment final','location','best');
  legend boxoff
  title('final, total, and cumulative correction')
 print('-dpsc','-append',graphname)
end

corrected.d13c.low = d13c.initial + cumsum(fin.inc);
details(iter).cumsum = cumsum(fin.inc);

% corrected.d13c.low(find(iyear<1850)) = d13c.low(find(iyear < 1850))

%% STEP 12 add high frequency back
[cyear,ai,bi] = intersect(corrected.d13c.year,iyear);
corrected.d13c.ratio(:,iter) = d13co(:,iter); % the full thing
corrected.d13c.ratio(ai,iter) = corrected.d13c.low + d13c.high; % now overwrite with corrected values

if verboseMode
 figure(4);
 plot (iyear,-c+corrected.d13c.ratio(ai,iter),'k','linewidth',1)
 legend('difference between original and corrected d13C series','location','best')
 title('Magnitude of PIN correction')
 print('-dpsc','-append',graphname)
end  
  
if plotMode & verboseMode
    figure(5); clf;
    plot(iyear,c,'k'); hold on
    plot(iyear,corrected.d13c.low,'b','linewidth',2)
    plot(iyear,d13c.low,'k','linewidth',2)
    plot(corrected.d13c.year(ai,1), corrected.d13c.ratio(ai,iter),'b');
    xlim([min(corrected.d13c.year(ai,1)) max(corrected.d13c.year(ai,1))]);
    legend('Measured','Corrected','location','best')
    legend boxoff
else
    figure; clf;
    plot(iyear,c,'k'); hold on
    plot(iyear,corrected.d13c.low,'b','linewidth',2)
    plot(iyear,d13c.low,'k','linewidth',2)
    plot(corrected.d13c.year(ai,1), corrected.d13c.ratio(ai,iter),'b');
    legend('Measured','Corrected','location','best')
    legend boxoff
end
 print('-dpsc','-append',graphname)

% add the first, uncorrected part back to the time series
if (firstyear < min(year))
 corrected.d13c.ratio = [cc(1:max(ic)-1)' corrected.d13c.ratio];
end

if verboseMode
    figure(6); clf
    plot(corrected.d13c.year,corrected.d13c.ratio(:,iter),'b'); hold on;
    plot(original.d13c.year,original.d13c.ratio(:,iter),'k');
    xlabel('year'); ylabel('d13C');
    xlim([min(corrected.d13c.year) max(corrected.d13c.year)]);
    legend('Full corrected','Full Original','location','best')
    legend boxoff
end
 print('-dpsc','-append',graphname)

if verboseMode
    details(iter).original = original;
    details(iter).fin = fin;
    details(iter).d13c = d13c;
    details(iter).p = pp;
    details(iter).ci = ci;
    details(iter).ca = ca;
end

clear delta cLow yhat pp

end  % termination of looping over different series, in columns




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Nested Subfunction for carbon dioxide concentration of the atmosphere

function [ca] = co2(iyear,jyear);

% CO2 The smoothed \emph{ca} curve for the appropriate time period
%
% Data are from Robertson et al. (2001) and McCarroll et al. (2007) and 
% based on Mauna Loa data from http://cdiac.ornl.gov/ftp/trends/co2/maunaloa.co2

caData = [...
1820	285	
1821	285	
1822	285	
1823	285	
1824	285	
1825	285	
1826	285	
1827	285	
1828	285	
1829	285	
1830	285	
1831	285	
1832	285	
1833	285	
1834	285	
1835	285	
1836	285	
1837	285	
1838	285	
1839	285	
1840	285	
1841	285	
1842	285	
1843	285	
1844	285	
1845	285	
1846	285.038334	
1847	285.102512	
1848	285.1726304	
1849	285.2487062	
1850	285.3307562	
1851	285.4187974	
1852	285.5128465	
1853	285.6129206	
1854	285.7190366	
1855	285.8312112	
1856	285.9490543	
1857	286.0722687	
1858	286.2010108	
1859	286.3354374	
1860	286.4757049	
1861	286.6219698	
1862	286.7743887	
1863	286.9331181	
1864	287.0983146	
1865	287.2701347	
1866	287.448735	
1867	287.637783	
1868	287.8398015	
1869	288.0532294	
1870	288.2765058	
1871	288.5080696	
1872	288.74636	
1873	288.9898159	
1874	289.2368764	
1875	289.4859805	
1876	289.7355671	
1877	289.9840754	
1878	290.2299444	
1879	290.480424	
1880	290.7425957	
1881	291.0146454	
1882	291.2947592	
1883	291.5811232	
1884	291.8719234	
1885	292.1653459	
1886	292.4595766	
1887	292.7528017	
1888	293.0432072	
1889	293.3289792	
1890	293.6160173	
1891	293.9106839	
1892	294.2118587	
1893	294.5184218	
1894	294.829253	
1895	295.1432322	
1896	295.4592392	
1897	295.776154	
1898	296.0928564	
1899	296.4082263	
1900	296.7211436	
1901	297.0304882	
1902	297.3414448
1903	297.6590506	
1904	297.9819645	
1905	298.3088451	
1906	298.6383513	
1907	298.9691419	
1908	299.2998756	
1909	299.6292112	
1910	299.9558075	
1911	300.2783233	
1912	300.5954172	
1913	300.9096966	
1914	301.2245886	
1915	301.539982	
1916	301.8557655	
1917	302.1718278	
1918	302.4880576	
1919	302.8043437	
1920	303.1205749	
1921	303.4366397	
1922	303.7524271	
1923	304.0678256	
1924	304.3827241	
1925	304.6864203	
1926	304.9711949	
1927	305.2414103	
1928	305.5014292	
1929	305.7556141	
1930	306.0083276	
1931	306.2639323	
1932	306.5267908	
1933	306.8012656	
1934	307.0917194	
1935	307.4025146	
1936	307.7380141	
1937	308.0644948	
1938	308.3520276	
1939	308.6106653	
1940	308.8504606	
1941	309.0814663	
1942	309.3137352	
1943	309.55732	
1944	309.8222737	
1945	310.1186489	
1946	310.4564984	
1947	310.845875	
1948	311.2968316	
1949	311.7636189	
1950	312.2012998	
1951	312.621145	
1952	313.034425	
1953	313.4524106	
1954	313.8863723	
1955	314.3475809	
1956	314.847307	
1957	315.3968211	
1958	316.007394	
1959	316.6902964	
1960	317.4253819	
1961	318.1854243	
1962	318.9710742	
1963	319.7829819	
1964	320.6217981	
1965	321.4881731	
1966	322.3827576	
1967	323.306202	
1968	324.2591568	
1969	325.2422725	
1970	326.2561997	
1971	327.3015888	
1972	328.375894	
1973	329.4764591	
1974	330.6037693	
1975	331.7583097	
1976	332.9405655	
1977	334.1510219	
1978	335.3901641	
1979	336.6584771	
1980	337.9564461	
1981	339.2845563	
1982	340.6432929	
1983	342.031623	
1984	343.4482469	
1985	344.8932493	
1986	346.3667151	
1987	347.868729	
1988	349.3993757	
1989	350.9587401	
1990	352.5469069	
1991	354.1639609	
1992	355.8099868	
1993	357.4850694	
1994	359.1892934	
1995	360.9227222	
1996	362.685303	
1997	364.476947	
1998	366.2975655	
1999	368.1470695	
2000	370.0253702	
2001	371.9323788	
2002	373.8680065
2003    375.83
2004    377.82
2005    379.85];

iyear = find(caData(:,1) >= iyear & caData(:,1) <= jyear); 
warning off
ca.year = caData(iyear,1); ca.ratio = caData(iyear,2);