Source code for MSQuant: PILstatistics2.cs, MSQlib1/src/utility/PILstatistics2.cs.

Table of contents page.

Home page for MSQuant.

/****************************************************************************
 * Copyright (C) 2008 Peter Mortensen and Matthias Mann                     *
 * This file is part of MSQuant.                                            *
 *                                                                          *
 * MSQuant is distributed under the terms of                                *
 * the GNU General Public License. See src/COPYING.TXT or                   *
 * <http://www.gnu.org/licenses/gpl.txt> for details.                       *
 *                                                                          *
 * MSQuant is free software; you can redistribute it                        *
 * and/or modify it under the terms of the GNU                              *
 * General Public License as published by the Free                          *
 * Software Foundation; either version 2 of the                             *
 * License, or (at your option) any later version.                          *
 *                                                                          *
 * MSQuant is distributed in the hope that it will be                       *
 * useful, but WITHOUT ANY WARRANTY; without even the                       *
 * implied warranty of MERCHANTABILITY or FITNESS FOR                       *
 * A PARTICULAR PURPOSE.  See the GNU General Public                        *
 * License for more details.                                                *
 *                                                                          *
 * You should have received a copy of the GNU General                       *
 * Public License along with MSQuant; if not, write to                      *
 * the Free Software Foundation, Inc., 59 Temple                            *
 * Place, Suite 330, Boston, MA  02111-1307  USA                            *
 *                                                                          *
 * Purpose: General statistics routines (that are not dependent on a        *
 *          particular application.).                                       *
 *                                                                          *
 ****************************************************************************/

/****************************************************************************
 *                               CEBI                                       *
 *                    Software Development Group                            *
 *                         Peter Mortensen                                  *
 *                E-mail: NUKESPAMMERSdrmortensen@get2netZZZZZZ.dk          *
 *                 WWW: http://www.cebi.sdu.dk/                             *
 *                                                                          *
 *  Program for post-processing of result from search in mass               *
 *    spectrometric data.                                                   *
 *                                                                          *
 *    FILENAME:   PILstatistics2.cs                                         *
 *    TYPE:       CSHARP                                                    *
 *                                                                          *
 * CREATED: PM 2008-03-14   Vrs 1.0. Translated from PILstatistics2.vb      *
 * UPDATED: PM 2008-xx-xx                                                   *
 *                                                                          *
 *                                                                          *
 ****************************************************************************/

using System; //For Math.
using System.Collections.Generic;
using System.Diagnostics; //For Trace


//****************************************************************************
//d$ <summary>
//d$ Purpose: Namespace for application independent and domain independent
//d$ utility classes (that could be reused in any application,
//d$ not just mass spectrometric applications.)
//d$ <see cref="T:VBXMLDoc.CVBXMLDoc" />.
//d$ <isUnitTest></isUnitTest>
//d$ <applicationname>XYZ</applicationname>
//d$ <author>Peter Mortensen</author>
//d$ <seealso>http://www.cebi.sdu.dk/</seealso>
//d$ <codetype>PLATFORM independent</codetype>
//d$ </summary>
namespace SDUPutility
{

    //Changed PM_REFACTOR 2007-10-04
    //Is this the right level? Should it be within Class SDUPstatistics?
    //Make this structure into a class??
    public struct statsStructure
    {
        //These two are only updated if different from null.
        //Better with built-in arrays?
        public List<double> xValues;
        //Better with built-in arrays?
        public List<double> yValues;
        //Dim xValues() As Double
        //Dim yValues() As Double

        public double minX;
        public double maxX;

        public double minY;
        public double maxY;

        public double sumOfX;

        public double sumOfabsY;
        public double sumOfY;
        public double sumOfSqY;


        //Changed PM_REFACTOR 2007-12-10
        public double sumOfSqX7;
        //Dim sumOfSqY As Double not yet. May not be needed.
        public double sumOfXYproduct;

        //Note:
        // linear regression can be computed:
        //
        // tri = N * sumOfSqX - sumOfX * sumOfX
        //
        // slope = (N * sumOfXYproduct - sumOfX * sumOfY ) /
        // tri 
        //
        // offset = (sumOfSqX * sumOfY - sumOfX * sumOfXYproduct ) /
        // tri 


        public int N;

        public int tag;
        // Primarily for debugging. If
        // we want to track values for a particular 
        // instance we can have if structures testing on
        // value and set breakpoints.


        //For convenience initialisation.
        public statsStructure(int aTag, bool aInitXValues, bool aInitYValues)
        {
            //aInitXValues and aInitYValues: whether to 
            //remember values or not. If only the sums or min/max
            //are needed then they should be false.

            tag = aTag;

            N = 0;
            sumOfabsY = 0.0;

            sumOfY = 0.0;
            sumOfSqY = 0.0;

            //Changed PM_CSHARP_DETECTED_PROBLEM 2008-03-13
            sumOfX = 0.0; //Discovered by C# !!!!!
            sumOfSqX7 = 0.0; //Discovered by C# !!!!!

            sumOfXYproduct = 0.0; //Discovered by C# !!!!!

            minX = 1E+19;
            maxX = -1E+19;

            minY = 1E+19;
            maxY = -1E+19;

            xValues = null;
            if (aInitXValues)
            {
                xValues = new List<double>();
            }

            yValues = null;
            if (aInitYValues)
            {
                yValues = new List<double>();
            }
        } //Convenience constructor.


        public void update(double anInX, double anInY)
        {
            N += 1;

            sumOfX += anInX;
            sumOfSqX7 += anInX * anInX;

            sumOfXYproduct += anInX * anInY;

            sumOfabsY += Math.Abs(anInY);

            sumOfY += anInY;
            sumOfSqY += anInY * anInY;


            if (minX > anInX)
            {
                minX = anInX;
            }
            if (maxX < anInX)
            {
                maxX = anInX;
            }

            if (minY > anInY)
            {
                minY = anInY;
            }
            if (maxY < anInY)
            {
                maxY = anInY;
            }

            if ((xValues != null))
            {
                xValues.Add(anInX);
            }

            if ((yValues != null))
            {
                yValues.Add(anInY);
            }
        } //update()


        //Changed PM_REFACTOR 2007-12-10
        public void linearRegress(out double anOutSlope, out double anOutOffset)
        {
            double tri = N * sumOfSqX7 - sumOfX * sumOfX;
            anOutSlope = (N * sumOfXYproduct - sumOfX * sumOfY) / tri;
            anOutOffset = (sumOfSqX7 * sumOfY - sumOfX * sumOfXYproduct) / tri;
        } //linearRegress()


        //Changed PM_PCP_CONSENSUS_PROFILE 2008-11-05
        public void average(out double anOutXaverage, out double anOutYaverage)
        {
            anOutXaverage = sumOfX / N;
            anOutYaverage = sumOfY / N;
        } //average()

        //Changed PM_CONSENSUSSCORE 2008-11-30
        public void stddev(out double anOutXstddev, out double anOutYstddev)
        {
            //From <http://en.wikipedia.org/wiki/Standard_deviation>,
            //  "Rapid calculation methods".

            //Change of square of (small) negative number? (due to round-off error?)
            anOutXstddev = Math.Sqrt( N * sumOfSqX7 - sumOfX * sumOfX);
            anOutYstddev = Math.Sqrt(N * sumOfSqY - sumOfY * sumOfY); 
        } //stddev()

    
    } //statsStructure


    //Changed PM_REFACTOR 2005-02-02. Moved here...
    public class SDUPstatistics
    {
        //Changed PM_MS3SCORING_HIGHN_ZERO 2004-10-27. Changed return type
        //to Double as Integer can not represent the range of numbers we need.
        //****************************************************************************
        //* Note: this is permutations (not combinations). Is this true?????         *
        //* E.g. (17 2) returns 136. *
        //* *
        //****************************************************************************
        public static double NoverK(int aN, int aK)
        {
            //Use integer multiplication for low values of aN and aK?

            double toReturn = 0.0;

            if (aN == 196)
            {
                int peter2 = 2;
            }
            
            //Zero in case of numeric error exception, etc.

            //If n > 30 Then
            // MsgBox("n too large")
            // Exit Function
            //End If
            //test

            // If k = 0 OR n = 0 Then
            // Return 0 'Score should be zero in this case.

            try
            {
                //Note: for high numbers we risk overflow. Divisions and
                //multiplications could be done alternately.

                //We currently use the straightforward way. It would be better
                //to use the one from <http://en.wikipedia.org/wiki/Combination>,
                //e.g.:
                //
                //    52     52   51   50   49   48   
                //        =  ---  ---  ---  ---  --- = 2598960  
                //    5       5    4    3    2    1    

                if (true) //Normally true, false for old way - e.g. for testing.
                {
                    double numerator   = aN;
                    double denominator = aK;

                    double result = 1.0;
                    int i;
                    for (i = 1; i <= aK; i++)
                    {
                        double frac = numerator / denominator;
                        result *= frac;

                        numerator -= 1.0;
                        denominator -= 1.0;
                    }
                    toReturn = result;
                }
                else
                {
                    //Old way
                    double result1 = 1.0;
                    int i;
                    for (i = 1; i <= aN; i++)
                    {
                        result1 *= i;
                    }

                    double result2 = 1;
                    for (i = 1; i <= aK; i++)
                    {
                        result2 *= i;
                    }
                    if (result2 < 1)
                    {
                        result2 = 1;
                    }

                    double result3 = 1;
                    for (i = 1; i <= aN - aK; i++)
                    {
                        result3 *= i;
                    }
                    if (result3 < 1)
                    {
                        result3 = 1;
                    }
                    toReturn = result1 / (result2 * result3);
                } //Old way.

                if (toReturn > 2150000000.0)
                {
                    //32 bit signed integer limit. Note:
                    // we are using doubles now.
                    int peter9 = 9;
                }
            }
            catch
            {
                Trace.Assert(false, 
                    "PIL ASSERT. Numerical error in NoverK(). This should never happen...");
            }
            return toReturn;
        } //NoverK()


        //Changed PM_REFACTOR 2005-07-28
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //*                                                                          *
        //****************************************************************************
        public static double absAverage(ref double[] aValues)
        {
            double sum = 0.0;
            int numbDataPts = aValues.Length;

            int i;
            for (i = 0; i <= aValues.Length - 1; i++)
            {
                double absVal = Math.Abs(aValues[i]);
                sum += absVal;
            }
            double toReturn = sum / numbDataPts;
            return toReturn;
        } //absAverage()


        //Changed PM_REFACTOR 2005-07-28
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //****************************************************************************
        public static void convertToRelativeError(
            ref double[] anInOutErrorValues, 
            ref double[] anInAbsoluteValues)
        {
            //double sum = 0.0;
            int numbDataPts = anInOutErrorValues.Length;
            int lastIndex = numbDataPts - 1;

            Trace.Assert(
                numbDataPts == anInAbsoluteValues.Length, "PIL ASSERT. <message>.");

            int i;
            for (i = 0; i <= lastIndex; i++)
            {
                double relError =
                    1000000.0 * anInOutErrorValues[i] / anInAbsoluteValues[i];
                anInOutErrorValues[i] = relError; //Replace absolute error
                                                  //with relative error.
            }
        } //convertToRelativeError()


        //Changed PM_REFACTOR 2007-10-08
        //****************************************************************************
        //* <placeholder for header> *
        //****************************************************************************
        public static double medianValue(ref double[] anInOutVector)
        {
            //Note: side-effects. The order of elements in the input
            //      vector may change...
            //
            //Use sorted keys to avoid changing the input?

            double toReturn = -1E+20;

            //Assume all elements are in use (would not be the case for pre-dim'ed).

            int numbDataPts = anInOutVector.Length;
            int lastIndex = numbDataPts - 1;
            Array.Sort(anInOutVector);

            int indexForMedian = numbDataPts / 2;
            //Note: integer
            // division. For even number of data points: integer division 
            // will round down, but as the result is used as a zero
            // based index the higher value will be picked.

            Trace.Assert(
                indexForMedian >= 0, 
                "PIL ASSERT. indexForMedian, " + 
                indexForMedian + 
                ", out of range [0; " + 
                lastIndex + 
                "].");

            toReturn = anInOutVector[indexForMedian];
            return toReturn;
        } //medianValue()


        //Changed PM_REFACTOR 2007-10-08
        //****************************************************************************
        //* <placeholder for header> *
        //****************************************************************************
        public static double medianValue(ref List<double> anInOutVector)
        {
            //Note: side-effects. The order of elements in the input
            // vector may change.
            //
            //Use sorted keys to avoid changing the input?

            double toReturn = -1E+20;

            //Assume all elements are in use (would not be the case for pre-dim'ed).

            int numbDataPts = anInOutVector.Count;
            int lastIndex = numbDataPts - 1;
            anInOutVector.Sort();

            int indexForMedian = numbDataPts / 2;
            Trace.Assert(
                indexForMedian >= 0, 
                "PIL ASSERT. indexForMedian, " + 
                indexForMedian + 
                ", out of range [0; " + 
                lastIndex + 
                "].");

            toReturn = anInOutVector[indexForMedian];

            return toReturn;
        } //medianValue()


        //Changed PM_REFACTOR 2005-02-02
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //* anOutMedianAbsoluteError3: half of the numerical errors are              *
        //* above this value.                                                        *
        //****************************************************************************
        public static void linearRegression(
            ref double[] aXvalues, 
            ref double[] aYvalues, 
            out double   anOutSlope,
            out double   anOutOffset,
            out double[] anOutErr2,
            out double   anOutMinAbsError2,
            out double   anOutMaxAbsError2,
            out double   anOutMedianAbsoluteError3, 
            ref bool[]   aUseErrors,
            out double   anOutMinimumXused,
            out double   anOutMaximumXused)
        {
            //Changed PM_RECALIB_OUTLIERDETECTION 2007-12-12
            //aUseErrors tells which points should be used for the
            //regression. aUseErrors can be null. Then it uses all points.

            //Note1 : anMedianAbsoluteError is used in ParsingCompleted(),
            // near "markOutLiers(yErrors, xySet1, medianErrorToUse)". To
            // not use the worst half of the errors.
            //
            //Note2: anOutMinAbsError and anOutMaxAbsError are for ***absolute***
            // error (no sign).

            int numbDataPts_Max = aYvalues.Length;
            //ss
            int lastIndex_Max = numbDataPts_Max - 1;
            anOutMinAbsError2 = 1000000000.0;
            anOutMaxAbsError2 = -1000000000.0;

            //Changed PM_REFACTOR 2007-12-10
            // 'calcMCRval(i) = anOutB + anOutA * aXvalues(i) so for
            // ' correction we mulitply the measured values like that.
            //
            // ' sx , sy ,sxoss , t ,st2
            // Dim sumMeasMCRval As Double = 0.0
            // Dim sumCalcMCRval As Double = 0.0
            // Dim meanMeasX As Double = 0.0
            // Dim devMeasX As Double = 0.0
            // Dim sumDevMeasXSquare As Double = 0.0
            //
            // Dim i As Integer
            // For i = 0 To lastIndex
            // sumMeasMCRval += aXvalues(i) 'sx += x[i];
            // sumCalcMCRval += aYvalues(i) 'sy += y[i]
            // Next
            //
            // meanMeasX = sumMeasMCRval / numbDataPts 'sxoss=sx/ss
            //
            // For i = 0 To lastIndex
            // devMeasX = aXvalues(i) - meanMeasX 't=x[i]-sxoss;
            // sumDevMeasXSquare += devMeasX * devMeasX 'st2 += t*t;
            // anOutA += devMeasX * aYvalues(i) '*anOutA += t*y[i];
            // Next i
            // anOutA /= sumDevMeasXSquare '*anOutA /= st2;
            // anOutB = (sumCalcMCRval - sumMeasMCRval * anOutA) / numbDataPts '*anOutB=(sy-sx*(*anOutA))/ss;
            statsStructure massValues = new statsStructure(1, false, false);
            //2 x false: no saving of
            // values - we already have 'em !

            //Changed PM_RECALIB_OUTLIERDETECTION 2007-12-12
            bool useAll = aUseErrors == null;

            int i;
            for (i = 0; i <= lastIndex_Max; i++)
            {
                //Changed PM_RECALIB_OUTLIERDETECTION 2007-12-12
                bool doUse = useAll || aUseErrors[i];
                //Note: relies
                // on Boolean short-circuit (must use it and order must
                // be preserved), else an exception would result when aUseErrors 
                // is null.

                if (doUse)
                {
                    double xVal = aXvalues[i];
                    double yVal = aYvalues[i];
                    massValues.update(xVal, yVal);
                }
                else
                {
                    int peter2 = 2; //For breakpoints.
                }
            }
            massValues.linearRegress(out anOutSlope, out anOutOffset);
            int lastIndex_Out = massValues.N - 1;

            anOutMinimumXused = massValues.minX;
            anOutMaximumXused = massValues.maxX;

            //Note: the direct formulas used in linearRegress() are susceptible to 
            // roundoff error (for large N?). The old way is better, but 
            // requires two loops. Should we revert to the old method?

            //Block. Post-recalibration calculations.
            {
                //Changed PM_CSHARP_NO_REDIM 2008-03-13
                double[] errVector = new double[lastIndex_Out + 1];

                //Changed PM_RECALIB_OUTLIERDETECTION 2007-12-12
                int outIndex = 0;

                // ReDim anOutErr2(lastIndex_Out)
                // ERROR: Not supported in C#: ReDimStatement 
                //Note: actual size will be less than 
                      //      input size if some points are not used
                      //      due content of aUseErrors.
                anOutErr2 = new double[lastIndex_Out+1];

                //Now test this. Get absolute average mass error after recalibration.
                double aveAbsPPMErr = 0.0;
                double aveAbsPPMErrBefore = 0.0;
                for (i = 0; i <= lastIndex_Max; i++)
                {
                    //Changed PM_RECALIB_OUTLIERDETECTION 2007-12-12
                    bool doUse = useAll || aUseErrors[i];
                    //Note: relies
                    // on Boolean short-circuit (must use it and order must
                    // be preserved), else an exception would result when aUseErrors 
                    // is null.

                    if (doUse)
                    {
                        double curX = aXvalues[i];
                        double curY = aYvalues[i];
                        double fx = anOutSlope * curX + anOutOffset;
                        double errY = curY - fx;
                        double absErrY = Math.Abs(errY);

                        anOutErr2[outIndex] = errY;
                        errVector[outIndex] = absErrY;

                        if (absErrY < anOutMinAbsError2)
                        {
                            anOutMinAbsError2 = absErrY;
                        }
                        if (absErrY > anOutMaxAbsError2)
                        {
                            anOutMaxAbsError2 = absErrY;
                        }

                        aveAbsPPMErrBefore += 1000000.0 * Math.Abs(curY - curX) / curY;
                        aveAbsPPMErr += 1000000.0 * absErrY / curY;

                        outIndex += 1;
                    }
                    else
                    {
                        int peter2 = 2;
                        //For breakpoints.
                    }
                }
                aveAbsPPMErrBefore /= aYvalues.Length;
                aveAbsPPMErr /= aYvalues.Length;

                //Changed PM_REFACTOR 2007-10-08
                // 'Changed PM_VS2005_SERIOUS 2006-06-15
                // ' Warning 58 Access of shared member, constant member, enum 
                // ' member or nested type through an instance; qualifying 
                // ' expression will not be evaluated. 
                // ' \src\main\massbase\PILstatist
                // ' ics.vb 252 13 
                // 'errVectorToSort.Sort(errVectorToSort)
                // Array.Sort(errVectorToSort)
                // 
                // Dim indexForMedian As Integer = numbDataPts \ 2
                // anMedianNumericalError = errVectorToSort(indexForMedian)

                //Note: side-effects in medianValue(): order of errVector is changed.
                anOutMedianAbsoluteError3 = medianValue(ref errVector);
            } //Block. Post-recalibration calculations.

        } //linearRegression()


        //Does not seem to be used at the moment... Outcommented.
        // '****************************************************************************
        // '* <placeholder for header> *
        // '* Convenience function for callers that do not require the error *
        // '* information. *
        // '****************************************************************************
        // Public Shared Sub linearRegression_NoErrorInfo2( _
        // ByRef aXvalues() As Double, ByRef aYvalues() As Double, _
        // ByRef anOutA As Double, ByRef anOutB As Double)
        // 
        // Dim dummyErr() As Double = null 'Keep compiler happy.
        // Dim dummyMinErr As Double
        // Dim dummyMaxErr As Double
        // Dim dummyMedianErr As Double
        // 
        // SDUPstatistics.linearRegression(aXvalues, aYvalues, anOutA, anOutB, _
        // dummyErr, dummyMinErr, dummyMaxErr, dummyMedianErr) 'B before A because.
        // End Sub 'linearRegression_NoErrorInfo


        //Changed PM_REFACTOR 2005-07-28. Moved to here.
        //****************************************************************************
        //* SUBROUTINE NAME: dumpXandYandYerrors *
        //d$ <summary>
        //d$ Purpose: xyz.
        //d$ 
        //d$ <see cref="T:VBXMLDoc.CVBXMLDoc" />.
        //d$ </summary>
        public static void dumpXandYandYerrors(
            ref string aHeadline, 
            ref double[] aXVals, 
            double aSlope, 
            double anOffset, 
            ref double[] anYVals, 
            ref string[] anInTags, 
            ref string anInXUnit, 
            ref string anInYUnit, 
            ref System.Text.StringBuilder anInOutStringBuffer)
        {
            anInOutStringBuffer.Append(aHeadline + PILInputOutput.LINEEND);

            //Changed PM_CORR_EXPORT_UNITS 2006-10-18
            anInOutStringBuffer.Append(
                "X " + anInXUnit + "\t" + 
                "Y " + anInYUnit + "\t" + 
                "X " + anInXUnit + "\t" + 
                "yErrors " + anInYUnit + "\t" + 
                "yPred " + anInYUnit + 
                PILInputOutput.LINEEND);

            int lastIndex = aXVals.Length - 1;
            int n;
            for (n = 0; n <= lastIndex; n++)
            {
                double X = aXVals[n];
                double Y = anYVals[n];
                double Ypred = aSlope * X + anOffset;
                //Dim errY As Double = anInYErrors(n)
                double errY = Y - Ypred;

                //anInOutStringBuffer.Append( _
                // X & vbTab & Y & vbTab & _
                // X & vbTab & errY & vbTab & _
                // Ypred & vbTab & PILInputOutput.LINEEND)

                anInOutStringBuffer.Append(X);
                anInOutStringBuffer.Append("\t");

                anInOutStringBuffer.Append(Y);
                anInOutStringBuffer.Append("\t");

                anInOutStringBuffer.Append(X);
                anInOutStringBuffer.Append("\t");

                anInOutStringBuffer.Append(errY);
                anInOutStringBuffer.Append("\t");

                anInOutStringBuffer.Append(Ypred);
                anInOutStringBuffer.Append("\t");

                if ((anInTags != null))
                {
                    anInOutStringBuffer.Append(anInTags[n]);
                }
                else
                {
                    int peter8 = 8;
                }
                anInOutStringBuffer.Append("\t");

                anInOutStringBuffer.Append(PILInputOutput.LINEEND);
            } //for loop.
        } //dumpXandYandYerrors()


        //Changed PM_REFACTOR 2008-03-25. binningSetup() made private and 2 new
        //  external functions added to reduce confusing on the client side.
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //****************************************************************************
        public static void binningSetup_computeBinSize(
            double anInAbsDeviationSum, 
            int anInSumCount, 
            double anInMin, 
            double anInMax, 
            out int[] anOutBinningArray, 
            out double anOutBinSize, 
            out int anOutBins)
        {
            //Hide this complexity from the client side:
            double someInOutBinSize = -17.0; //Flag that the internal function
            //  should compute the bin, assuming Gaussian like distribution.
            //  It will set the bin size so the peak in the distribution is
            //  sufficient resolved, e.g. approx 5 bins.

            binningSetup2(
                anInAbsDeviationSum, anInSumCount, 
                anInMin, anInMax, 
                out anOutBinningArray,
                ref someInOutBinSize,
                out anOutBins);
            anOutBinSize = someInOutBinSize;
        } //binningSetup_computeBinSize()


        //Changed PM_REFACTOR 2008-03-25. binningSetup() made private and 2 new
        //  external functions added to reduce confusing on the client side.
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //****************************************************************************
        public static void binningSetup_inputBinSize(
            double anInMin,
            double anInMax,
            double anInBinSize,
            out int[] anOutBinningArray,
            out int anOutBins)
        {
            //Hide this complexity from the client side:
            double dummyDeviationSum = 0.0;
            int sumCount = 1;

            binningSetup2(
                dummyDeviationSum, sumCount,
                anInMin, anInMax,
                out anOutBinningArray,
                ref anInBinSize,
                out anOutBins);
        } //binningSetup_computeBinSize()


        //Changed PM_RECALIB_VISUAL_BINNING 2007-10-01
        //****************************************************************************
        //* <placeholder for header>                                                 *
        //****************************************************************************
        private static void binningSetup2(
            double anInAbsDeviationSum, 
            int anInSumCount, 
            double anInMin, 
            double anInMax, 
            out int[] anOutBinningArray, 
            ref double anInOutBinSize, 
            out int anOutBins)
        {
            //Note: anInDeviationSum and anInSumCount only makes sense if
            // values are clustered close to zero, like gaussian
            // distributed mass errors. The two parameters are only 
            // used if anInOutBinSize is negative, indicating that 
            // the binsize should be determined in this function.

            int binsInt = 1;
            //Safe value.
            double binSize2 = 1.0;
            //Safe value.

            Trace.Assert(anInSumCount >= 0, "PIL ASSERT. No data! In binningSetup().");

            if (anInSumCount > 0)
            {
                binSize2 = -2E+20;
                //Value: for error detection.

                if (anInOutBinSize < 0.0)
                {
                    double averageAbs = anInAbsDeviationSum / anInSumCount;
                    binSize2 = averageAbs * 0.1 * 4 * 1.5;
                    anInOutBinSize = binSize2;
                    //Returned computed value.
                }
                else
                {
                    binSize2 = anInOutBinSize;
                    //Use passed value for bin size.
                }

                Trace.Assert(binSize2 > -1E+20, "PIL ASSERT. <message>.");

                double range = anInMax - anInMin;
                if (range < 0.0)
                {
                    range = 1.0;
                    //Safe value.
                }

                //Changed PM_BINNING_INFINITY 2008-05-26
                //Avoid infinities...
                double bins;
                if (binSize2 != 0.0) //Use more robust comparison?
                {
                    bins = range / binSize2;
                }
                else
                {
                    bins = 1; //Some safe number.
                }

                if (bins < 1E-06)
                {
                    int peter5 = 5;
                }

                //Changed PM_BINNINGBUG_RECALIB_VISUAL 2008-03-04. For exactly
                // integer result in the above (e.g. 0.0) it should round up.
                //binsInt = CInt(bins + 0.5) 'Nearest higher integer.
                //binsInt = CInt(bins + 0.5001) 'Nearest higher integer.
                //binsInt = (int)(bins + 0.5001); //Nearest higher integer.
                //  Incorrect by translation-tool !!!!!!! -
                //  <http://labs.developerfusion.co.uk/convert/vb-to-csharp.aspx>
                //
                binsInt = (int)(bins + 1.0001); //Nearest higher integer.

                //Changed PM_CINT_TROUBLE 2008-03-19
                //asdasdadsasd; //CInt(), a Visual Basic thing, is ***not***
                //  the same as an int cast. CInt() rounds to nearest
                //  whereas (int) rounds down.


                //Changed PM_BINNINGBUG_RECALIB_VISUAL 2008-03-04
                Trace.Assert(binsInt > 0, "PIL ASSERT. No bins!.");
            }
            else
            {
                int peter2 = 2;
                //No data. E.g. if a peptide
                // filter let no peptides through.
            }

            int lastBinIndex = binsInt - 1;

            // ReDim anOutBinningArray(lastBinIndex)
            //   Why wasn't it ReDim Preserve??????????????? Because there 
            //   is one and only one dimension? No. Explanation: we are
            //   returning an array for binning with the proper size and
            //   containing all zeroes.
            // ERROR: Not supported in C#: ReDimStatement 
            anOutBinningArray = new int[lastBinIndex + 1]; //Implicit
            //  initialisation???

            anOutBins = binsInt;
        } //binningSetup()

    } //SDUPstatistics

} //SDUPutility


    

    

Generated by script codePublish.pl at 2009-01-05T15:20:59.