classdef Measurement
    % MEASUREMENT Represents a value with an uncertainty
    %   Written by Grant Elliott, June 2011
    %   Uses nonlinear fitting function written by the MIT 8.13/8.14 staff
    %
    %   A measurement is constructed using any of the following syntax: 
    %   Measurement()     Constructs 0±0
    %   Measurement(a)    Constructs a±0 if a is a singleton
    %   Measurement(a)    Constructs a   if a is a Measurement
    %   Measurement(a,b)  Constructs a±b
    %
    %   Can be assembled into matrices, for which concatenation, reference, index, assignment,
    %   transpose, and display all work as expected.  Overloaded operators assume perfect indepedence 
    %   (zero covariance).  Use functional operations for covariance control.
    % 
    %   Warning: These calculations are not necessarily associative, particularly if covariance
    %   matrices are not used carefully. Consider using the highest level function available which
    %   suits your needs.
    %
    %   Measurement Properties: (a.PROPERTY)
    %   Value                    - Gets the value of a
    %   Uncertainty              - Gets the uncertainty of a
    %
    %   Measurement Methods:    (a.METHOD(ARGS))
    %   toString                 - Format as string
    %   distribution             - Generate monte carlo distribution of points about a Measurement
    %   plus                     - Binary addition        +
    %   uplus                    - Unary plus             +
    %   minus                    - Binary subtraction     -
    %   uminus                   - Unary minus            -
    %   times                    - Array multiplication   .*
    %   mtimes                   - Matrix multiplication  *  (Singletons only)
    %   rdivide                  - Right array division   ./
    %   ldivide                  - Left array division    .\
    %   mrdivide                 - Right matrix division  /  (Singletons only)
    %   mldivide                 - Left matrix division   \  (Singletons only)
    %   power                    - Array power            .^ (Numeric exponents only)
    %   mpower                   - Matrix power           ^  (Singletons with numeric exponents only)
    %
    %   Measurement Static Methods: (Measurement.METHOD(ARGS))
    %   add                      - Array add with a given covariance
    %   addCorrelated            - Array add assuming maximal correlation               (Worst case)
    %   addAnticorrelated        - Array add assuming maximal negative correlation
    %   addUncorrelated          - Array add assuming zero correlation
    %   subtract                 - Array subtract with a given covariance
    %   subtractCorrelated       - Array subtract assuming maximal correlation
    %   subtractAnticorrelated   - Array subtract assuming maximal negative correlation (Worst case)
    %   subtractUncorrelated     - Array subtract assuming zero correlation
    %   multiply                 - Array multiply with a given covariance
    %   multiplyCorrelated       - Array multiply assuming maximal correlation          (Worst case)
    %   multiplyAnticorrelated   - Array multiply assuming maximal negative correlation
    %   multiplyUncorrelated     - Array multiply assuming zero correlation
    %   divide                   - Array divide with a given covariance
    %   divideCorrelated         - Array divide assuming maximal correlation
    %   divideAnticorrelated     - Array divide assuming maximal negative correlation   (Worst case)
    %   divideUncorrelated       - Array divide assuming zero correlation
    %   sum                      - Sum over array with a given covariance matrix
    %   sumCorrelated            - Sum over array assuming maximal covariance matrix    (Worst case)
    %   sumUncorrelated          - Sum over array assuming zero covariance
    %   prod                     - Prod over array with a given covariance matrix
    %   prodCorrelated           - Prod over array assuming maximal covariance matrix   (Worst case)
    %   prodUncorrelated         - Prod over array assuming zero covariance
    %   meanUncorrelated         - Mean of array assuming zero covariance
    %   meanUncorrelatedWeighted - Weighted Mean of array assuming zero covariance
    %   functionDistribution     - Propagates uncertainties through a function by Monte Carlo
    %   plotFunctionDistribution - Plots a region associated with a function's uncertainty
    %   errorBar                 - Plot error bars in one or two dimensions
    %   fit                      - Non-linear curve fitting with error bars both axes
    %
    %   For help on a specific property or method, type:
    %               help Measurement.(METHOD OR PROPERTY)
	
    properties (Access=protected)
        value = 0;
        uncertainty = 0;
    end
    
    methods
        function obj = Measurement(varargin)
            % MEASUREMENT
            %   Measurement()     Constructs 0±0
            %   Measurement(a)    Constructs a±0 if a is a singleton
            %   Measurement(a)    Constructs a   if a is a Measurement
            %   Measurement(a,b)  Constructs a±b
            if (nargin==0)
                obj.value = 0;
                obj.uncertainty = 0;
            elseif (nargin==1 && isa(varargin{1},'Measurement'))
                a=varargin{1};
                obj=repmat(Measurement, size(a));
                for i=1:numel(a) %There has to be a way to do this that doesn't require an explicit iteration
                    obj(i).value = a(i).Value;
                    obj(i).uncertainty = a(i).Uncertainty;
                end
            elseif (nargin==1 && isnumeric(varargin{1}))
                v=varargin{1};
                obj=repmat(Measurement, size(v));
                for i=1:numel(v)
                    obj(i).value = v(i);
                    obj(i).uncertainty = 0;
                end
            elseif (nargin==2 && isnumeric(varargin{1}) && isnumeric(varargin{2}))
                v=varargin{1};
                u=varargin{2};
                if numel(u)==1
                    u=u*ones(size(v));
                end
                if numel(v)==1
                    v=v*ones(size(v));
                end
                if ~all(size(u)==size(v))
                    error('measurement:constructor:invalidArgument', 'Values and uncertainties must have same length of be singleton');
                end
                obj=repmat(Measurement, size(v));
                for i=1:numel(v)
                    obj(i).value = v(i);
                    obj(i).uncertainty = u(i);
                end
            else
                error('measurement:constructor:invalidArgument', 'Cannot construct a measurement from these arguments');
            end
        end
        
        function v = Value(a)
            % VALUE  The center of the distribution represented by the measurement
            v=[a.value];
        end
        
        function u = Uncertainty(a)
            % UNCERTAINTY  The standard deviation of the distribution represented by the measurement
            u=[a.uncertainty];
        end 
        
        function u = RelativeUncertainty(a)
            % RELATIVEUNCERTAINTY  The Uncertainty represented as a fraction of the Value.
            u=a.Uncertainty./a.Value;
        end
        
        function display(a)
            % DISPLAY   Display array of measurements
            %   Displays an array of measurements using ± notation
            if ~isequal(get(0,'FormatSpacing'),'compact')
                disp(' ')
            end
            disp([inputname(1) ' =']);
            if ~isequal(get(0,'FormatSpacing'),'compact')
                disp(' ')
            end
            if (numel(a)==1)
                fprintf('     %s\n',a.toString);
            else
                dim=sprintf('%dx',size(a));
                dim=dim(1:end-1);
                fprintf('     [%s Matrix of Measurements]\n',dim);
            end
            if ~isequal(get(0,'FormatSpacing'),'compact')
                disp(' ')
            end
        end
        
        function a=toString(a)
            % TOSTRING    Format as string
            a=sprintf('%g%s%g',a.Value,char(177),a.Uncertainty);
        end
        
        function c = plus(a,b)
            % + Plus
            %   A+B adds Measurement arrays A and B assuming zero covariance and is equivalant to
            %   Measurement.addUncorrelated(A,B). For covariance control, use Measurement.add
            %   instead.
            % 
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also add, addCorrelated, addAnticorrelated, addUncorrelated.
            c = Measurement.addUncorrelated(a,b);
        end
        
        function c = minus(a,b)
            % - Plus
            %   A-B subtracts Measurement array B from Measurement array A assuming zero covariance
            %   and is equivalant to Measurement.subtractUncorrelated(A,B). For covariance control,
            %   use Measurement.subtract instead.
            % 
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also subtract, subtractCorrelated, subtractAnticorrelated, subtractUncorrelated.
            c = Measurement.subtractUncorrelated(a,b);
        end
        
        function c = uminus(a)
            % + Unary minus 
            %   -A negates the Value of Measurement array A
            c = 0-Measurement(a);
        end
        
        function c = uplus(a)
            % + Unary plus 
            %   +A for Measurement arrays is A
            c = Measurement(a);
        end
       
        function c = times(a,b)
            % .* Array multiply
            %   A.*B denotes element-by-element multiplication of Measurement arrays A and B
            %   assuming zero covariance and is equivalant to Measurement.multiplyUncorrelated(A,B).
            %   For covariance control, use Measurement.multiply instead.  A and B must have the
            %   same dimensions unless one of them is singleton.
            % 
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also multiply, multiplyCorrelated, multiplyAnticorrelated, multiplyUncorrelated.
            c = Measurement.multiplyUncorrelated(a,b);
        end
        
        function c = mtimes(a,b)
            %   * Matrix multiply
            %   A*B is equivalent to A.*B if one of A or B is singleton.  True matrix multiplication
            %   is not supported.
            %
            %   See also times, multiply, multiplyCorrelated, multiplyAnticorrelated, 
            %   multiplyUncorrelated.
            if (numel(a)==1 || numel(b)==1)
                c = times(a,b);
            else
                error('Measurement:mtimes:IllegalArgument','Matrix multiplication only supported for singleton multiplication');
            end
        end
        
        function c = rdivide(a,b)
            % ./ Right array divide
            %   A./B denotes element-by-element division of Measurement array A by B assuming
            %   zero covariance and is equivalant to Measurement.divideUncorrelated(A,B). For
            %   covariance control, use Measurement.divide instead.  A and B must have the same
            %   dimensions unless one of them is singleton.
            % 
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also ldivide, divide, divideCorrelated, divideAnticorrelated,
            %   divideUncorrelated.
            c = Measurement.divideUncorrelated(a,b);
        end
        
        function c = ldivide(a,b)
            % .\ Left array divide
            %   A.\B denotes element-by-element division of Measurement array B by A assuming
            %   zero covariance and is equivalant to Measurement.divideUncorrelated(B,A). For
            %   covariance control, use Measurement.divide instead.  A and B must have the same
            %   dimensions unless one of them is singleton.
            % 
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also rdivide, divide, divideCorrelated, divideAnticorrelated,
            %   divideUncorrelated.
            c = rdivide(b,a);
        end
        
        function c = mrdivide(a,b)
            % / Right matrix divide
            %   A/B is equivalent to A./B if one of A or B is singleton.  True matrix division is
            %   not supported.
            %
            %   See also rdivide, ldivide, divide, divideCorrelated, divideAnticorrelated,
            %   divideUncorrelated.
            if (numel(a)==1 || numel(b)==1)
                c = rdivide(a,b);
            else
                error('Measurement:mrdivide:IllegalArgument','Matrix division only supported for singleton division');
            end
        end
        
        function c = mldivide(a,b)
            % \ Left matrix divide
            %   A\B is equivalent to A.\B if one of A or B is singleton.  True matrix division is
            %   not supported.
            %
            %   See also rdivide, ldivide, divide, divideCorrelated, divideAnticorrelated,
            %   divideUncorrelated.
            if (numel(a)==1 || numel(b)==1)
                c = ldivide(a,b);
            else
                error('Measurement:mldivide:IllegalArgument','Matrix division only supported for singleton division');
            end
        end
        
        function c = power(a,b)
            % .^ Array power
            %   A.^B denotes element-by-element powers of Measurement arrays A to the power of
            %   numeric array B.  A and B must have the same dimensions unless one of them is
            %   singleton.
            a = Measurement(a);
            if (~isnumeric(b))
                error('Measurement:power:IllegalArgument','Only numeric exponents are supported');
            end
            val = a.Value.^b;
            c = Measurement(val,val.*b.*a.Uncertainty./a.Value);
        end
        
        function c = mpower(a,b)
            % ^ Matrix power
            %   A^B is equivalent to A.^B if A is a singleton Measurement and B is a numeric 
            %   singleton.  True matrix power is not supported.
            %
            %   See also power.
            if (numel(a)~=1 || numel(b)~=1)
                error('Measurement:mpower:IllegalArgument','Only singleton arguments are supported');
            end
            if (~isnumeric(b))
                error('Measurement:mpower:IllegalArgument','Only numeric exponents are supported');
            end
            a = Measurement(a);
            val = a.Value.^b;
            c = Measurement(val,val.*b.*a.Uncertainty./a.Value);
        end
        
        function dist=distribution(a,N)
            % DISTRIBUTION    Generate monte carlo distribution of points about a Measurement
            %   dist=a.distribution(N) returns a random vector of N values from the normal distribution
            %   characterized by a.  a must be a row vector.
            %
            %   See also functionDistribution.
            a=Measurement(a);
            if size(a,1)~=1
                a=a.';
            end
            if size(a,1)~=1
                error('Measurement:distribution:IllegalArgument','Must act on a vector');
            end
            dist=repmat(a.Value,N,1)+repmat(a.Uncertainty,N,1).*randn(N,numel(a));
        end
            
    end
    
    methods(Access=public, Static=true)        
        function c = add(a,b,covariance)
            % ADD   Sum of Measurement arrays with covariance
            %   add(A,B,COVAR) adds two Measurement arrays with a given covariance.  A, B, and COVAR
            %   must have the same size or be singleton.  COVAR must be numeric.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also plus, addCorrelated, addAnticorrelated, addUncorrelated.
            [a,b,covariance] = Measurement.arithmaticArgumentCheck(a,b,covariance,'add');
            c = Measurement(a.Value+b.Value,sqrt(a.Uncertainty.^2+b.Uncertainty.^2+2*covariance));
        end
        
        function c = addCorrelated(a,b)
            % ADDCORRELATED   Sum of Measurement arrays with maximal covariance
            %   addCorrelated(A,B) adds two Measurement arrays assuming maximal covariance.  A and B
            %   must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   For addition, assuming maximal covariance yields the worst case error.
            %
            %   See also plus, add, addAnticorrelated, addUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.add(a,b,a.Uncertainty.*b.Uncertainty);
        end
        
        function c = addAnticorrelated(a,b)
            % ADDANTICORRELATED   Sum of Measurement arrays with maximal negative covariance
            %   addAnticorrelated(A,B) adds two Measurement arrays assuming maximal negative
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also plus, add, addCorrelated, addUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.add(a,b,-a.Uncertainty.*b.Uncertainty);
        end
        
        function c = addUncorrelated(a,b)
            % ADDUNCORRELATED   Sum of Measurement arrays with zero covariance
            %   add(A,B) adds two Measurement arrays assuming zero covariance.  A and B must have 
            %   the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also plus, add, addCorrelated, addAnticorrelated.
            c = Measurement.add(a,b,0);
        end
        
        function c = subtract(a,b,covariance)
            % SUBTRACT   Difference of Measurement arrays with covariance
            %   subtract(A,B,COVAR) subtracts Measurement array B from A with a given covariance.
            %   A, B, and COVAR must have the same size or be singleton.  COVAR must be numeric.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also minus, subtractCorrelated, subtractAnticorrelated, subtractUncorrelated.
            [a,b,covariance] = Measurement.arithmaticArgumentCheck(a,b,covariance,'subtract');
            c = Measurement(a.Value-b.Value,sqrt(a.Uncertainty.^2+b.Uncertainty.^2-2*covariance));
        end
        
        function c = subtractCorrelated(a,b)
            % SUBTRACTCORRELATED   Difference of Measurement arrays with maximal covariance
            %   subtractCorrelated(A,B) subtracts Measurement array B from A assuming maximal
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also minus, subtract, subtractAnticorrelated, subtractUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.subtract(a,b,a.Uncertainty.*b.Uncertainty);
        end
        
        function c = subtractAnticorrelated(a,b)
            % SUBTRACTANTICORRELATED   Difference of Measurement arrays with maximal negative 
            %   subtractAntiCorrelated(A,B) subtracts Measurement array B from A assuming maximal
            %   negative covariance.  A and B must have the same size or be singleton.
            %
            %   For subtraction, assuming maximal negative covariance yields the worst case error.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also minus, subtract, subtractCorrelated, subtractUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.subtract(a,b,-a.Uncertainty.*b.Uncertainty);
        end
        
        function c = subtractUncorrelated(a,b)
            % SUBTRACTUNCORRELATED   Difference of Measurement arrays with maximal negative 
            %   subtractUncorrelated(A,B) subtracts Measurement array B from A assuming zero
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also minus, subtract, subtractCorrelated, subtractAnticorrelated.
            c = Measurement.subtract(a,b,0);
        end
        
        function c = multiply(a,b,covariance)
            % MULTIPLY   Product of Measurement arrays with covariance
            %   multiply(A,B,COVAR) multiplies two Measurement arrays with a given covariance.  A,
            %   B, and COVAR must have the same size or be singleton.  COVAR must be numeric.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also times, multiplyCorrelated, multiplyAnticorrelated, multiplyUncorrelated.
            [a,b,covariance] = Measurement.arithmaticArgumentCheck(a,b,covariance,'multiply');
            val = a.Value.*b.Value;
            c = Measurement(val, val.*sqrt((a.Uncertainty./a.Value).^2+(b.Uncertainty./b.Value).^2+2*covariance./(a.Value.*b.Value)));
        end
        
        function c = multiplyCorrelated(a,b)
            % MULTIPLYCORRELATED   Product of Measurement arrays with maximal covariance
            %   multiplyCorrelated(A,B) multiplys Measurement arrays A and B assuming maximal
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   For multiplication, assuming maximal covariance yields the worst case error.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also times, multiply, multiplyAnticorrelated, multiplyUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.multiply(a,b,a.Uncertainty.*b.Uncertainty);
        end
        
        function c = multiplyAnticorrelated(a,b)
            % MULTIPLYANTICORRELATED  Product of Measurement arrays with maximal negative covariance
            %   multiplyAnticorrelated(A,B) multiplys Measurement arrays A and B assuming maximal
            %   negative covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also times, multiply, multiplyCorrelated, multiplyUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.multiply(a,b,-a.Uncertainty.*b.Uncertainty);
        end
        
        function c = multiplyUncorrelated(a,b)
            % MULTIPLYANTICORRELATED   Product of Measurement arrays with zero covariance
            %   multiplyAnticorrelated(A,B) multiplys Measurement arrays A and B assuming zero
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also times, multiply, multiplyCorrelated, multiplyAnticorrelated.
            c = Measurement.multiply(a,b,0);
        end
        
        function c = divide(a,b,covariance)
            % DIVIDE   Quotient of Measurement arrays with covariance
            %   divide(A,B,COVAR) divides Measurement array A by B with a given covariance.  A, B,
            %   and COVAR must have the same size or be singleton.  COVAR must be numeric.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also rdivide, ldivide, divideCorrelated, divideAnticorrelated,
            %   divideUncorrelated.
            [a,b,covariance] = Measurement.arithmaticArgumentCheck(a,b,covariance,'divide');
            val = a.Value./b.Value;
            c = Measurement(val, val.*sqrt((a.Uncertainty./a.Value).^2+(b.Uncertainty./b.Value).^2-2*covariance./(a.Value.*b.Value)));
        end
        
        function c = divideCorrelated(a,b)
            % DIVIDECORRELATED   Quotient of Measurement arrays with maximal covariance
            %   divideCorrelated(A,B) divides Measurement array A by B assuming maximal
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also rdivide, ldivide, divide, divideAnticorrelated, divideUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.divide(a,b,a.Uncertainty.*b.Uncertainty);
        end
        
        function c = divideAnticorrelated(a,b)
            % DIVIDEANTICORRELATED   Quotient of Measurement arrays with maximal negative covariance
            %   divideAnticorrelated(A,B) divides Measurement array A by B assuming maximal
            %   negative covariance.  A and B must have the same size or be singleton.
            %
            %   For division, assuming maximal negative covariance yields the worst case error.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also rdivide, ldivide, divide, divideCorrelated, divideUncorrelated.
            a = Measurement(a);
            b = Measurement(b);
            c = Measurement.divide(a,b,-a.Uncertainty.*b.Uncertainty);
        end
        
        function c = divideUncorrelated(a,b)
            % DIVIDEUNCORRELATED   Quotient of Measurement arrays with zero covariance
            %   divideUncorrelated(A,B) divides Measurement array A by B assuming zero
            %   covariance.  A and B must have the same size or be singleton.
            %
            %   If A or B is numeric, it will be treated as a Measurement with zero uncertainty.
            %
            %   See also rdivide, ldivide, divide, divideCorrelated, divideAnticorrelated.
            c = Measurement.divide(a,b,0);
        end
        
        function c = sum(a,covariance)
            % SUM   Sum of elements
            %   sum(V,COVAR) sums a Measurement Array using a covariance matrix.
            %   Covariance must be square with size given by the number of elements in V.  Only the
            %   upper triangle is used.  The diagonal is implicitly unity and the lower triangle
            %   implicitly mirrors the upper triangle.
            %
            %   See also sumCorrelated, sumUncorrelated.
            [a,covariance] = Measurement.accumulatorArgumentCheck(a,covariance,'sum');
            val=sum(a.Value);
            var=sum(a.Uncertainty.^2);
            for i=1:numel(a)
                for j=i+1:numel(a)
                    var=var+2*covariance(i,j);
                end
            end
            c=Measurement(val,sqrt(var));
        end
        
        function c = sumCorrelated(a)
            % SUMCORRELATED     Sum of elements
            %   sumCorrelated(V,COVAR) sums a Measurement Array assuming maximal covariance.
            %
            %   For addition, assuming maximal covariance yields the worst case error.
            %
            % 	See also sum, sumUncorrelated.
            covariance = zeros(numel(a));
            a=reshape(a,[1,numel(a)]);
            a=Measurement(a);
            for i=1:numel(a)
                for j=i+1:numel(a)
                    covariance(i,j)=a(i).Uncertainty.*a(j).Uncertainty;
                end
            end
            c = Measurement.sum(a,covariance);
        end
        
        function c = sumUncorrelated(a)
            % SUMUNCORRELATED     Sum of elements
            %   sumUncorrelated(V,COVAR) sums a Measurement Array assuming zero covariance.
            %
            %   See also sum, sumCorrelated.
            c = Measurement.sum(a,zeros(numel(a)));
        end
        
        function c = prod(a,covariance)
            % PROD  Product of elements
            %   prod(V,COVAR) sums a Measurement Array using a given covariance matrix.
            %   Covariance must be square with size given by the number of elements in V.  
            %   Only the upper triangle is used.  The diagonal is implicitly unity and the lower
            %   triangle implicitly mirrors the upper triangle.
            %
            %   See also prodCorrelated, prodUncorrelated.
            [a,covariance] = Measurement.accumulatorArgumentCheck(a,covariance,'prod');
            val=prod(a.Value);
            normvar=sum((a.Uncertainty./a.Value).^2);
            for i=1:numel(a)
                for j=i+1:numel(a)
                    normvar=normvar+2*covariance(i,j)/(a(i).Value*a(j).Value);
                end
            end
            c=Measurement(val,val*sqrt(normvar));
        end
        
        function c = prodCorrelated(a)
            %   PRODCORRELATED     Sum of elements
            %   prodCorrelated(V,COVAR) sums a Measurement Array assuming maximal covariance.
            %
            %   For multiplication, assuming maximal covariance yields the worst case error.
            %
            %   See also prod, prodUncorrelated.
            covariance = zeros(numel(a));
            a=reshape(a,[1,numel(a)]);
            a=Measurement(a);
            for i=1:numel(a)
                for j=i+1:numel(a)
                    covariance(i,j)=a(i).Uncertainty.*a(j).Uncertainty;
                end
            end
            c = Measurement.prod(a,covariance);
        end
        
        function c = prodUncorrelated(a)
            % PRODUNCORRELATED     Sum of elements
            %   prodUncorrelated(V,COVAR) sums a Measurement Array assuming zero covariance.
            %
            %   See also prod, and prodCorrelated.
            c = Measurement.prod(a,zeros(numel(a)));
        end
        
        function [x, y] = meanUncorrelated(a)
            % MEANUNCORRELATED    Mean of elements
            %   [X,Y] = meanUncorrelated(A) takes the arithmetic mean of the values of A assuming no
            %   correlation.  X has an uncertainty representing the distribution of samples
            %   (standard deviation), while Y has an uncertainty representing the distribution of
            %   sampling means (standard error).  A must be either entirely absolute values
            %   (Uncertainty equals zero) or entirely uncertain values (Uncertainty does not equal
            %   zero).
            %
            %   For absolute values:
            %   n        = numel(A)
            %   w        = 1/n
            %   mu       = sum(w.*A.Value)         = mean(a.Value)
            %   sigma    = sqrt(n/(n-1)*sum(w.*(A.Value-mu).^2)) = std(a.Value)
            %   sigma_mu = sigma/sqrt(n)
            %   X        = mu ± sigma
            %   Y        = mu ± sigma_mu
            %
            %   For uncertain Values:
            %   n        = numel(A)
            %   lambda   = (1./A.Uncertainty^2) / sum(/A.Uncertainty^2)
            %   mu       = sum(lambda.*A.Value)
            %   sigma    = sqrt(n/(n-1)*sum(lambda.*(A.Value-mu).^2))
            %   sigma_mu = sqrt(1/sum(1/A.Uncertainty^2))
            %   X        = mu ± sqrt(sigma^2+sigma_mu^2)
            %   Y        = mu ± sigma_mu
            %
            %   See also meanUncorrelatedWeighted.
            [x,y] = Measurement.meanUncorrelatedWeighted(a, ones(size(a)));
        end
        
        function [x, y] = meanUncorrelatedWeighted(a,w)
            % MEANUNCORRELATEDWEIGHTED    Weighted mean of elements
            %   [X,Y] = meanUncorrelatedWeighted(A,W) takes the arithmetic mean of the values of A
            %   with weighting factors W assuming no correlation.  X has an uncertainty representing
            %   the distribution of samples (standard deviation), while Y has an uncertainty
            %   representing the distribution of sampling means (standard error).  A must be either
            %   entirely absolute values (Uncertainty equals zero) or entirely uncertain values 
            %   (Uncertainty does not equal zero).  A and w must have the same number of elements.
            %
            %   For absolute values:
            %   n        = numel(A)
            %   w        = W/sum(W)
            %   mu       = sum(w.*A.Value)
            %   sigma    = sqrt(n/(n-1)*sum(w.*(A.Value-mu).^2))
            %   sigma_mu = sigma/sqrt(n)
            %   X        = mu ± sigma
            %   Y        = mu ± sigma/sqrt(n)
            %
            %   For uncertain Values:
            %   n        = numel(a)
            %   w        = W/sum(W)
            %   lambda   = (1./A.Uncertainty^2) / sum(w./A.Uncertainty^2)
            %   mu       = sum(w.*lambda.*A.Value)
            %   sigma    = sqrt(n/(n-1)*sum(w.*lambda.*(A.Value-mu).^2))
            %   sigma_mu = sqrt(1/sum(n*w./A.Uncertainty^2))
            %   X        = mu ± sqrt(sigma^2+sigma_mu^2)
            %   Y        = mu ± sigma_mu
            %
            %   See also meanUncorrelated.
            a=Measurement(a);
            if ~all(a.Uncertainty==0) && any(a.Uncertainty==0)
                error('Measurement:mean:IllegalArgument','Cannot calculate mean of mixed uncertain and absolute values');
            end
            if numel(a) ~= numel(w)
                error('Measurement:mean:IllegalArgument','Measurement array and weight array must have the same number of elements');
            end
            absolute = all(a.Uncertainty==0);
            n = numel(a);
            a = reshape(a,[1,n]);
            w = reshape(w,[1,n]);
            w = w/sum(w);
            if absolute
                lambda = 1;
            else
                lambda = (1./a.Uncertainty.^2) / sum(w./a.Uncertainty.^2);
            end
            mu = sum(w.*lambda.*a.Value);
            sigma = sqrt(n/(n-1)*sum(w.*lambda.*(a.Value-mu).^2));
            if absolute
                sigma_mu = sigma/sqrt(n);
            	x = Measurement(mu,sigma);
            else
                sigma_mu = sqrt(1/sum(n*w./a.Uncertainty.^2));
                x = Measurement(mu,sqrt(sigma^2+sigma_mu^2));
            end
            y = Measurement(mu,sigma_mu);
        end
        
        function [l,p]=errorBar(x,y,varargin)
            % ERRORBAR    Plot error bars in one or two dimensions for Measurements
            %   [l,p] = errorBar(x,y) plots one or two dimensional error bars where x and y are vectors of
            %   Measurements or numbers.  Returns handles for the error bar lines and markers.
            %
            %   Accepts the following optional parameters:
            %   Marker                 - Marker to use at central values (Defaults to 'o')
            %   MarkerSize             - Size of Markers in points (Defaults to 6)
            %   MarkerEdgeColor        - Length 3 Array specifying color of marker edges (Defaults to [1,1,1])
            %   MarkerFaceColor        - Length 3 Array specifying color of marker faces (Defaults to [1,1,1])
            %   ErrorColor             - Length 3 Array specifying color of error lines (Defaults to [1,1,1])
            %   ErrorLineWidth         - Width of error lines in pixels (Defaults to 1)
            %
            %   See also fit, plotFunctionDistribution
            if numel(x)~=numel(y)
                error('Measurement:ErrorBar:ArgumentMismatch','X and Y must have same length');
            end
            p=inputParser;
            p.addParamValue('Marker','o', @(x) ischar(x) && numel(x)==1);
            p.addParamValue('MarkerSize',6, @(x) isnumeric(x) && numel(x)==1 && x>=0);
            p.addParamValue('MarkerEdgeColor',[0,0,0], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('MarkerFaceColor',[0,0,0], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('ErrorColor',[0,0,0], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('ErrorLineWidth',1, @(x) isnumeric(x) && numel(x)==1 && x>0);
            p.parse(varargin{:});
            settings = p.Results;
            x=Measurement(x);
            y=Measurement(y);
            p=plot(x.Value,y.Value,'linestyle','none','Marker',settings.Marker,'MarkerSize',settings.MarkerSize,'MarkerEdgeColor',settings.MarkerEdgeColor,'MarkerFaceColor',settings.MarkerFaceColor);
            linex=zeros(2,2*numel(x));
            liney=zeros(2,2*numel(x));
            linex(1,1:2:end) = x.Value-x.Uncertainty;
            linex(2,1:2:end) = x.Value+x.Uncertainty;
            liney(1,1:2:end) = y.Value;
            liney(2,1:2:end) = y.Value;
            liney(1,2:2:end) = y.Value-y.Uncertainty;
            liney(2,2:2:end) = y.Value+y.Uncertainty;
            linex(1,2:2:end) = x.Value;
            linex(2,2:2:end) = x.Value;
            l=line(linex,liney,'linewidth',settings.ErrorLineWidth,'Color',settings.ErrorColor);
        end
        
        function [y] = functionDistribution(f,x,param,N)
            % FUNCTIONDISTRIBUTION    Propagates uncertainties of parameters and values through a function by Monte Carlo
            %   [y] = functionDistribution(f,x,param,N) evaluates f(x) at param's nominal value and
            %   with lower and upper 1-sigma values given uncertainties in param using Monte Carlo methods
            %   f(x,param) is a handle to a function taking a vector of indepedent variable values x and
            %   a vector of parameters param as arguments.
            %   x is a vector of independent variables
            %   param is a vector of parameters
            %   N is the number of Monte Carlo points to generate.
            %   If not specified, N defaults to 10000).
            %   y is a vector of Measurements
            %
            %   See also distribution, plotFunctionDistribution, fit, errorBar
            if nargin<4
                N=10000;
            end
            x=Measurement(x);
            x=reshape(x,1,numel(x));
            param = Measurement(param);
            param=reshape(param,1,numel(param));
            v=f(x.Value,param.Value);
            xdist = x.distribution(N);
            paramdist = param.distribution(N);
            ydist = zeros(N,numel(x));
            for i=1:N;
                ydist(i,:)=f(xdist(i,:),paramdist(i,:));
            end
            u=std(ydist,0,1);
            y=Measurement(v,u);
        end
        
        function [l,e,f] = plotFunctionDistribution(x,y,varargin)
            % PLOTFUNCTIONDISTRIBUTION    Plot a region corresponding to distribution of a function
            %   [l,e,f] = plotFunctionDistribution(x,y,paramname,paramval,...) plots a measurement y against
            %   a value x with a line, 1-sigma lines, and shading between them
            %
            %   Returns handles to the central line, uncertainty lines, and filled region
            %
            %   Accepts the following optional parameters:
            %   LineStyle         - Style to use for central line (Defaults to '-')
            %   LineColor         - Length 3 Array specifying color of central line (Defaults to [0,0,1])
            %   LineWidth         - Width of central line (Defaults to 1)
            %   BorderLineStyle   - Style to use for 1-sigma line (Defaults to '--')
            %   BorderLineColor   - Length 3 Array specifying color of 1-sigma line (Defaults to [0,0,1])
            %   BorderLineWidth   - Width of 1-sigma line (Defaults to 1)
            %   FillColor         - Length 3 Array specifying color of filled region (Defaults to [0,0,1])
            %   FillAlpha         - Alpha value to apply to filled region (Defaults to .1)
            %   LineVisible       - Weather to plot central line (Defaults to true)
            %   BorderLineVisible - Whether to plot 1-sigma lines (Defaults to true)
            %   FillVisible       - Whether to fill region between 1-sigma lines (Defaults to true)
            %
            %   See also functionDistribution, fit, errorBar
            p=inputParser;
            p.addParamValue('LineStyle','-', @(x) ischar(x) && numel(x)>0);
            p.addParamValue('BorderLineStyle','--', @(x) ischar(x) && numel(x)>0);
            p.addParamValue('LineColor',[0,0,1], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('BorderLineColor',[0,0,1], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('FillColor',[0,0,1], @(x) isnumeric(x) && numel(x)==3 && max(x)<=1 && min(x)>=0);
            p.addParamValue('LineWidth',1, @(x) isnumeric(x) && numel(x)==1 && x>0);
            p.addParamValue('BorderLineWidth',1, @(x) isnumeric(x) && numel(x)==1 && x>0);
            p.addParamValue('LineVisible',true, @(x) islogical(x) && numel(x)==1);
            p.addParamValue('BorderLineVisible',true, @(x) islogical(x) && numel(x)==1);
            p.addParamValue('FillVisible',true, @(x) islogical(x) && numel(x)==1);
            p.addParamValue('FillAlpha',.1, @(x) isnumeric(x) && numel(x)==1 && x>0 && x<1);
            p.parse(varargin{:});
            settings = p.Results;
            if isa(x,'Measurement')
                x=x.Value;
            end
            y=Measurement(y);
            x=reshape(x,1,numel(x));
            y=reshape(y,1,numel(y));
            l=0;
            e=0;
            f=0;
            if settings.FillVisible
                xs=[x,fliplr(x)];
                ys=[y.Value-y.Uncertainty,fliplr(y.Value+y.Uncertainty)];
                ind= and(~isnan(xs), ~isnan(ys));
                xs=xs(ind);
                ys=ys(ind);
                %This is a workaround for a matlab bug that affects rendering with opengl.  Actual transparency is problematic.
                %Unfortunately, this won't look great if regions overlap
                %c=settings.FillColor*settings.FillAlpha+[1,1,1]*(1-settings.FillAlpha);
                %f=patch(xs,ys,-eps*ones(size(xs)),c,'edgecolor',c);
                %Here's how it should be done
                f=patch(xs,ys,settings.FillColor,'edgeAlpha',0,'faceAlpha',settings.FillAlpha);
            end
            if settings.BorderLineVisible
                e=plot([x;x].',[y.Value+y.Uncertainty;y.Value-y.Uncertainty].','linestyle',settings.BorderLineStyle,'linewidth',settings.BorderLineWidth,'color',settings.BorderLineColor);
            end
            if settings.LineVisible
                l=plot(x,y.Value,'linestyle',settings.LineStyle,'linewidth',settings.LineWidth,'color',settings.LineColor);
            end    
        end
        
        function [param,chisq] = fit(x,y,fitfun,param0)
            % FIT    Non-linear curve fitting with error bars in dependent and independent variables
            %   [param,chisq] = fit(x,y,fitfun,param0) fits data in x and y to a function given by fitfun.
            %   x and y are vectors of measurements of equal length.  Every point must have uncertainty
            %   at least one direction
            %
            %   fitfun(x,param) is a handle to a function taking a vector of indepedent variable values x and
            %   a vector of parameters param as arguments.
            %   param0 is an initial guess as to parameter values (For many fitfun, 0 is not recommended
            %   as an initial guess even if it seems most likely)
            %   param is a vector of final parameter values as Measurements with uncertainties
            %   chisq is the Chi^2 value associated with the fit
            %
            %   Uncertainties in the independent variable are handled using the estimated variances (EV) method
            %
            %   See also functionDistribution, plotFunctionDistribution, errorBar.
            
            x=Measurement(x);
            y=Measurement(y);
            
            if numel(x)~=numel(y)
                error('Measurement:ArgumentMismatch','X and Y must have same length');
            end
            if any(and(x.Uncertainty==0,y.Uncertainty==0))
                error('Measurement:IllegalArgument','All points must have uncertainty in either x or y');
            end
            
            if all(x.Uncertainty==0)
                % No uncertainty in independent variable
                [a,aerr,chisq] = Measurement.fitnonlin(x.Value,y.Value,y.Uncertainty,fitfun,param0);
            else
                % Uncertainty in the independent variable is handled by EV method
                % Ow, my brain
                a=fminsearch(@(a)Measurement.fitHelper(x,y,fitfun,a),param0);
                [~,a,aerr,chisq]=Measurement.fitHelper(x,y,fitfun,a);
            end
            param=Measurement(a,aerr);
        end
        
    end
    
    methods (Static=true, Access=protected)
        function [a,b,covar] = arithmaticArgumentCheck(a,b,covar,funcStr)
            a = Measurement(a);
            b = Measurement(b);
            s=sprintf('Measurement:%s:IllegalArgument',funcStr);
            if (~isnumeric(covar))
                error(s,'Covariance must be numeric');
            end
            if ((numel(a)~=1 && numel(b)~=1) && any(size(a)~=size(b))) || ((numel(covar)~=1 && numel(b)~=1) && any(size(covar)~=size(b))) || ((numel(a)~=1 && numel(covar)~=1) && any(size(a)~=size(covar)))
                error(s,'All arguments must be the same size or be singleton');
            end
            if any(abs(covar)>a.Uncertainty.^2) || any(abs(covar)>b.Uncertainty.^2)
                error(s,'Covariance must be less than variance of both terms to be summed');
            end 
        end
        
        function [a,covar] = accumulatorArgumentCheck(a,covar,funcStr)
            s=sprintf('Measurement:%s:IllegalArgument',funcStr);
            if (numel(covar)==1)
                covar = covar * ones(numel(covar));
            elseif (any(size(covar)~=[numel(a),numel(a)]))
                error(s,'Covariance must be a singleton or a square matrix with size given by the number of elements in the array to sum');
            end
            if (~isnumeric(covar))
                error(s,'Covariance must be numeric');
            end
            a=reshape(a,[1,numel(a)]);
            a=Measurement(a);
            for i=1:numel(a)
                if ~isnan(covar(i,i)) && covar(i,i)~=0 && covar(i,i)~=a(i).Uncertainty.^2 %#ok<BDSCA>
                    error(s,'Diagonal entries of covariance matrix must be NaN, zero, or the square of the corresponding uncertainty');
                end
                for j=i+1:numel(a)
                    if abs(covar(i,j)) > a(i).Uncertainty.*a(j).Uncertainty; %#ok<BDSCA>
                        error(s,'Covariance matrix inconsistent with uncertainties');
                    end
                    if ~isnan(covar(j,i)) && covar(j,i)~=0 && covar(i,j)~=covar(j,i);
                        error(s,'Lower triangle entries of covariance matrix must be NaN, zero, or equal to the corresponding upper triangle entry');
                    end
                end
            end
        end
        
        function [err, anext, anexterr, chisq] = fitHelper(x,y,fitfun,a)
            % returns an error based on how much the guessed a differs from the a calculated a after using it to propagate x uncertainties into y uncertainties
            N=numel(x);
            Nsamp=N*20;
            dx=(max(x.Value)-min(x.Value))/Nsamp;
            xsamp=[x.Value-dx, x.Value+dx];
            ysamp=fitfun(xsamp,a);
            dydxsamp = (ysamp(N+1:end)-ysamp(1:N))/(2*dx);
            yerreff = sqrt(y.Uncertainty.^2 + (dydxsamp.*x.Uncertainty).^2);
            [anext, anexterr, chisq]=Measurement.fitnonlin(x.Value, y.Value, yerreff, fitfun, a);
            err = sum(((a-anext)./anexterr).^2);
        end
        
        function [a,aerr,chisq,yfit] = fitnonlin(x,y,sig,fitfun,a0)
            % FITNONLIN Fit a nonlinear function to data.
            %    Written by the MIT 8.13/8.14 staff
            %    [a,aerr,chisq,yfit] = fitnonlin(x,y,sig,fitfun,a0)
            %
            %    Inputs:  x -- the x data to fit
            %             y -- the y data to fit
            %             sig -- the uncertainties on the data points
            %             fitfun -- handle of function to fit to
            %             a0 -- the initial guess at the parameters
            %
            %    Outputs: a -- the best fit parameters
            %             aerr -- the errors on these parameters
            %             chisq -- the final value of chi-squared
            %             yfit -- the value of the fitted function
            %                     at the points in x
            
            % First set up the parameters needed by the algorithm
            % These parameters can be varied if you have reason to believe your fit is
            % converging to quickly or that you are in a local minima of the chi
            % square.
            stepdown = 0.1;
            stepsize = abs(a0)*0.01; % the amount each parameter will be varied by in each iteration
            chicut = 0.00001;  % maximum differential allowed between successive chi sqr values
            a = a0;
            chi2 = Measurement.calcchi2(x,y,sig,fitfun,a);
            chi1 = chi2+chicut*2;
            % keep looking while the value of chi^2 is changing
            while (abs(chi2-chi1))>chicut
                [anew,stepsum] = Measurement.gradstep(x,y,sig,fitfun,a,stepsize,stepdown);
                a = anew;
                stepdown = stepsum;
                chi1 = chi2;
                chi2 = Measurement.calcchi2(x,y,sig,fitfun,a);
            end
            % calculate the returned values
            aerr = Measurement.sigparab(x,y,sig,fitfun,a,stepsize);
            chisq = Measurement.calcchi2(x,y,sig,fitfun,a);
            yfit = feval(fitfun,x,a);
        end

        function chi2 = calcchi2(x,y,sig,fitfun,a)
            %    Written by the MIT 8.13/8.14 staff
            chi2 = sum( ((y-feval(fitfun,x,a)) ./sig).^2);
        end
        
        function [anew,stepsum] = gradstep(x,y,sig,fitfun,a,stepsize, stepdown)
            %    Written by the MIT 8.13/8.14 staff
            % calculates the (negative) chi^2 gradient at the current point in parameter space, and moves in that direction until a minimum is found
            % returns the new value of the parameters and the total length travelled
            chi2 = Measurement.calcchi2(x,y,sig,fitfun,a);
            grad = Measurement.calcgrad(x,y,sig,fitfun,a,stepsize);
            chi3 = chi2*1.1;
            chi1 = chi3;
            % cut down the step size until a single step yields a decrease in chi^2
            stepdown = stepdown*2;
            j = 0;
            while chi3>chi2
                stepdown = stepdown/2;
                anew = a+stepdown*grad;
                chi3 = Measurement.calcchi2(x,y,sig,fitfun,anew);
                j=j+1;
            end
            stepsum = 0;
            % keep going until a minimum is passed
            while chi3<chi2
                stepsum = stepsum+stepdown;
                chi1 = chi2;
                chi2 = chi3;
                anew = anew+stepdown*grad;
                chi3 = Measurement.calcchi2(x,y,sig,fitfun,anew);
            end
            % approximate the minimum as a parabola
            step1 = stepdown*((chi3-chi2)/(chi1-2*chi2+chi3)+.5);
            anew = anew - step1*grad;
        end
        
        function grad = calcgrad(x,y,sig,fitfun,a, stepsize)
            %    Written by the MIT 8.13/8.14 staff
            % calculates the (negative) gradient at a point in parameter space
            f = 0.01;
            [~, nparm] = size(a);
            grad = a;
            chisq2 = Measurement.calcchi2(x,y,sig,fitfun,a);
            for i=1:nparm
                a2 = a;
                da = f*stepsize(i);
                a2(i) = a2(i)+da;
                chisq1 = Measurement.calcchi2(x,y,sig,fitfun,a2);
                grad(i) = chisq2-chisq1;
            end
            t = sum(grad.^2);
            grad = stepsize.*grad/sqrt(t);
        end
        
        function err=sigparab(x,y,sig,fitfun,a,stepsize)
            %    Written by the MIT 8.13/8.14 staff
            % calculates the errors on the final fitted parameters by approximating the minimum as parabolic in each parameter.
            [~, nparm] = size(a);
            err = a;
            chisq2 = Measurement.calcchi2(x,y,sig,fitfun,a);
            for i=1:nparm
                a2 = a;
                da = stepsize(i);
                a2(i) = a2(i)+da;
                chisq3 = Measurement.calcchi2(x,y,sig,fitfun,a2);
                a2(i) = a2(i)-2*da;
                chisq1 = Measurement.calcchi2(x,y,sig,fitfun,a2);
                err(i)=da*sqrt(2/(chisq1-2*chisq2+chisq3));
            end
        end
    end
end