IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4568


Ignore:
Timestamp:
Jul 15, 2005, 2:06:33 PM (21 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psFunctions.c

    r4567 r4568  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2005-07-15 23:19:41 $
     9 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2005-07-16 00:06:32 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    18181818    must exist n+1 points in "knots".
    18191819 
    1820 XXX: Ensure that domain[i+1] != domain[i]
     1820XXX: Is this really needed anymore?
    18211821 
    1822 XXX: What should be the defualty type for knots be?  psF32 is assumed.
     1822XXX: Ensure that knots[i+1] != knots[i]
     1823 
     1824XXX: What should be the default type for knots be?  psF32 is assumed.
    18231825 *****************************************************************************/
    1824 psSpline1D *psSpline1DAlloc( int numSplines,
    1825                              int order,
    1826                              float min,
    1827                              float max)
    1828 {
     1826psSpline1D *psSpline1DAlloc(int numSplines,
     1827                            int order,
     1828                            float min,
     1829                            float max)
     1830{
     1831    printf("HEY: psSpline1DAlloc()\n");
    18291832    PS_ASSERT_INT_NONNEGATIVE(numSplines, NULL);
    18301833    PS_ASSERT_INT_NONNEGATIVE(order, NULL);
    18311834    PS_ASSERT_FLOAT_NON_EQUAL(max, min, NULL);
    18321835
    1833     psSpline1D *tmp = NULL;
    1834     psS32 i;
    1835     psF32 tmpDomain;
    1836     psF32 width;
    1837 
    1838     tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    1839     tmp->n = numSplines;
    1840 
    1841     tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    1842     for (i=0;i<numSplines;i++) {
    1843         (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    1844     }
    1845 
    1846     // This should be set by the psVectorFitSpline1D()
    1847     tmp->p_psDeriv2 = NULL;
    1848 
    1849     tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32);
    1850     width = (max - min) / ((psF32) numSplines);
    1851 
    1852     tmp->knots->data.F32[0] = min;
    1853     tmpDomain = min+width;
    1854     for (i=1;i<numSplines+1;i++) {
    1855         tmp->knots->data.F32[i] = tmpDomain;
    1856         tmpDomain+= width;
    1857     }
    1858     tmp->knots->data.F32[numSplines] = max;
    1859 
    1860     psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree);
    1861     return(tmp);
    1862 }
    1863 
     1836    psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
     1837    tmpSpline->n = numSplines;
     1838
     1839    //
     1840    // XXX: We might have to allocate single or double polynomials depending on the type
     1841    // of the psVector bounds.  For now, all knots and spline polynomials are 32-bit.
     1842    //
     1843    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
     1844    for (psS32 i=0;i<numSplines;i++) {
     1845        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     1846    }
     1847
     1848    // This will be computed by psVectorFitSpline1D()
     1849    tmpSpline->p_psDeriv2 = NULL;
     1850
     1851    //
     1852    // XXX:Ensure that the knots are distinct, and monotonic.
     1853    //
     1854    tmpSpline->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32);
     1855    if (tmpSpline->knots == NULL) {
     1856        printf("BAD\n");
     1857    } else {
     1858        printf("GOOD\n");
     1859    }
     1860    psF32 width = (max - min) / ((psF32) numSplines);
     1861    tmpSpline->knots->data.F32[0] = min;
     1862    for (psS32 i=1;i<numSplines;i++) {
     1863        tmpSpline->knots->data.F32[i] = min + (width * (psF32) i);
     1864    }
     1865    tmpSpline->knots->data.F32[numSplines] = max;
     1866
     1867    psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree);
     1868    return(tmpSpline);
     1869}
    18641870
    18651871/*****************************************************************************
    1866 XXX: What should be the defualty type for knots be?  psF32 is assumed.
     1872XXX: Is there a psLib function for this?
     1873 *****************************************************************************/
     1874psVector *PsVectorDup(psVector *in)
     1875{
     1876    psVector *out = psVectorAlloc(in->n, in->type.type);
     1877
     1878    if (in->type.type == PS_TYPE_F32) {
     1879        for (psS32 i = 0 ; i < in->n ; i++) {
     1880            out->data.F32[i] = in->data.F32[i];
     1881        }
     1882    } else if (in->type.type == PS_TYPE_F64) {
     1883        for (psS32 i = 0 ; i < in->n ; i++) {
     1884            out->data.F64[i] = in->data.F64[i];
     1885        }
     1886    } else {
     1887        printf("XXX: Generate an error here.\n");
     1888        return(NULL);
     1889    }
     1890    return(out);
     1891}
     1892
     1893/*****************************************************************************
     1894XXX: What should be the default type for knots, spline polys?  psF32 is assumed.
    18671895 *****************************************************************************/
    18681896psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
    18691897                                   int order)
    18701898{
     1899    printf("HEY: psSpline1DAllocGeneric()\n");
    18711900    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
    18721901    PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL);
     1902    PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
    18731903    PS_ASSERT_INT_NONNEGATIVE(order, NULL);
    18741904
    1875     psSpline1D *tmpSpline = NULL;
    1876     unsigned int i;
    1877     unsigned int numSplines;
    1878 
    1879     tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    1880 
    1881     numSplines = bounds->n - 1;
     1905    psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
     1906    psS32 numSplines = bounds->n - 1;
    18821907    tmpSpline->n = numSplines;
    18831908
     1909    //
     1910    // XXX: We might have to allocate single or double polynomials depending on the type
     1911    // of the psVector bounds.  For now, all knots and spline polynomials are 32-bit.
     1912    //
    18841913    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    1885     for (i=0;i<numSplines;i++) {
     1914    for (psS32 i=0;i<numSplines;i++) {
    18861915        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    18871916    }
    18881917
    1889     // This should be set by the psVectorFitSpline1D()
     1918    // This will be computed by psVectorFitSpline1D()
    18901919    tmpSpline->p_psDeriv2 = NULL;
    18911920
    1892     tmpSpline->knots = psVectorAlloc(bounds->n, PS_TYPE_F32);
    1893 
    1894     for (i=0;i<bounds->n;i++) {
    1895         tmpSpline->knots->data.F32[i] = bounds->data.F32[i];
    1896         if (i<(bounds->n-1)) {
    1897             if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) {
    1898                 psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n");
    1899             }
    1900         }
    1901     }
    1902 
    1903     psMemSetDeallocator(tmpSpline,(psFreeFunc)spline1DFree);
     1921    //
     1922    // Ensure that all knots are distinct.
     1923    // XXX:Ensure that the knots are monotonic.
     1924    //
     1925    for (psS32 i=0;i<bounds->n-1;i++) {
     1926        if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) {
     1927            psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n");
     1928            return(NULL);
     1929        }
     1930    }
     1931    tmpSpline->knots = PsVectorDup((psVector *) bounds);
     1932    if (tmpSpline->knots == NULL) {
     1933        printf("BAD\n");
     1934    } else {
     1935        printf("GOOD\n");
     1936    }
     1937
     1938    psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree);
    19041939    return(tmpSpline);
    19051940}

  • trunk/psLib/src/math/psFunctions.h

    r4540 r4568  
    1212 *  @author GLG, MHPCC
    1313 *
    14  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    15  *  @date $Date: 2005-07-12 19:12:01 $
     14 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2005-07-16 00:06:32 $
    1616 *
    1717 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    434434typedef struct
    435435{
    436     unsigned int n;                    ///< The number of spline polynomials
     436    int n;                             ///< The number of spline pieces
    437437    psPolynomial1D **spline;           ///< An array of n pointers to the spline polynomials
     438    psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
    438439    psF32 *p_psDeriv2;                 ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    439     psF32 *domains;                    ///< The boundaries between each spline piece.  Size is n+1.
    440     psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
     440    psF32 *p_psDomains;                ///< The boundaries between each spline piece.  Size is n+1.
    441441}
    442442psSpline1D;

  • trunk/psLib/src/math/psMinimize.c

    r4540 r4568  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.124 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-07-12 19:12:01 $
     11 *  @version $Revision: 1.125 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-07-16 00:06:32 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    166166 
    167167XXX: spline->knots must be psF32
     168 
     169XXXX: Remove this for next code shipment.
    168170 *****************************************************************************/
    169171/*
    … …  
    262264        PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
    263265    }
    264 
    265266    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
    266267            "---- psVectorFitSpline1D() begin ----\n");
    … …  
    425426
    426427    }
     428
     429    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     430            "---- psVectorFitSpline1D() end ----\n");
     431    return(mySpline);
     432}
     433
     434
     435
     436
     437
     438
     439
     440
     441
     442/*****************************************************************************
     443psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
     444 
     445xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
     446of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
     447the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
     448simply be psFree() at the end.
     449 
     450XXX: nKnots makes no sense.  This number is always equal to the size of the x
     451an y vectors.
     452 
     453XXX: How do we specify the spline order?  For now, order=3.
     454 *****************************************************************************/
     455#define PS_XXX_SPLINE_ORDER 3
     456psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
     457                                   const psVector* y,        ///< Coordinates
     458                                   int nKnots)
     459{
     460    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
     461    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     462    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     463    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
     464
     465    //
     466    // The following code ensures that xPtr points to a psF32 version of the
     467    // ordinate data.
     468    //
     469    psVector *xF32 = NULL;
     470    psVector *xPtr = NULL;
     471    if (x != NULL) {
     472        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
     473        if (PS_TYPE_F64 == x->type.type) {
     474            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     475            for (psS32 i = 0 ; i < x->n ; i++) {
     476                xF32->data.F32[i] = (psF32) x->data.F64[i];
     477            }
     478            xPtr = xF32;
     479        } else if (PS_TYPE_F32 == x->type.type) {
     480            xPtr = (psVector *) x;
     481        } else {
     482            printf("XXX: Gen Error message: x is wrong type.\n");
     483            return(NULL);
     484        }
     485    } else {
     486        // If x==NULL, create an x32 vector with x values set to (0:n).
     487        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     488        for (psS32 i = 0 ; i < x->n ; i++) {
     489            xF32->data.F32[i] = (psF32) i;
     490        }
     491        xPtr = xF32;
     492    }
     493
     494    //
     495    // If y is of type psF64, then create a new vector yF32 and convert the
     496    // y elements.  Regardless of y's type, we create a yPtr which will be
     497    // used in the remainder of this function.
     498    //
     499    psVector *yF32 = NULL;
     500    psVector *yPtr = NULL;
     501    if (PS_TYPE_F64 == y->type.type) {
     502        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     503        for (psS32 i = 0 ; i < y->n ; i++) {
     504            yF32->data.F32[i] = (psF32) y->data.F64[i];
     505        }
     506        yPtr = yF32;
     507    } else {
     508        yPtr = (psVector *) y;
     509    }
     510
     511    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
     512
     513    psS32 numSplines = nKnots - 1;
     514    psF32 tmp;
     515    psF32 H;
     516    psS32 i;
     517    psF32 slope;
     518    // XXX: get rid of x32 and y32 (this is from old code)
     519    psVector *x32 = xPtr;
     520    psVector *y32 = yPtr;
     521
     522
     523    // If these are linear splines, which means their polynomials will have
     524    // two coefficients, then we do the simple calculation.
     525    if (1 == PS_XXX_SPLINE_ORDER) {
     526        for (i=0;i<mySpline->n;i++) {
     527            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
     528                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
     529            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
     530                                              (slope * mySpline->knots->data.F32[i]);
     531
     532            (mySpline->spline[i])->coeff[1] = slope;
     533            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     534                    "---- mySpline %d coeffs are (%f, %f)\n", i,
     535                    (mySpline->spline[i])->coeff[0],
     536                    (mySpline->spline[i])->coeff[1]);
     537        }
     538        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     539                "---- Exiting psVectorFitSpline1D()()\n");
     540        return((psSpline1D *) mySpline);
     541    }
     542
     543    //
     544    // Check if these are cubic splines (n==4).  If not, psError.
     545    //
     546    if (3 != PS_XXX_SPLINE_ORDER) {
     547        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     548                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
     549        return(NULL);
     550    }
     551
     552    //
     553    // If we get here, then we know these are cubic splines.  We first
     554    // generate the second derivatives at each data point.
     555    //
     556    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
     557    for (i=0;i<y32->n;i++)
     558        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     559                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
     560
     561    //
     562    // We generate the coefficients of the spline polynomials.  I can't
     563    // concisely explain how this code works.  See above function comments
     564    // and Numerical Recipes in C.
     565    //
     566    for (i=0;i<numSplines;i++) {
     567        H = x32->data.F32[i+1] - x32->data.F32[i];
     568        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     569                "x data (%f - %f) (%f)\n",
     570                x32->data.F32[i],
     571                x32->data.F32[i+1], H);
     572        //
     573        // ******** Calculate 0-order term ********
     574        //
     575        // From (1)
     576        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
     577        // From (2)
     578        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
     579        // From (3)
     580        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     581        tmp-= (x32->data.F32[i+1] / H);
     582        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
     583        ((mySpline->spline[i])->coeff[0])+= tmp;
     584        // From (4)
     585        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     586        tmp+= (x32->data.F32[i] / H);
     587        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
     588        ((mySpline->spline[i])->coeff[0])+= tmp;
     589
     590        //
     591        // ******** Calculate 1-order term ********
     592        //
     593        // From (1)
     594        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
     595        // From (2)
     596        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
     597        // From (3)
     598        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     599        tmp+= (1.0 / H);
     600        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
     601        ((mySpline->spline[i])->coeff[1])+= tmp;
     602        // From (4)
     603        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     604        tmp-= (1.0 / H);
     605        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
     606        ((mySpline->spline[i])->coeff[1])+= tmp;
     607
     608        //
     609        // ******** Calculate 2-order term ********
     610        //
     611        // From (3)
     612        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
     613        // From (4)
     614        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
     615
     616        //
     617        // ******** Calculate 3-order term ********
     618        //
     619        // From (3)
     620        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
     621        // From (4)
     622        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     623
     624        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     625                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     626        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     627                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     628        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     629                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     630        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     631                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     632
     633    }
     634
     635    psFree(xF32);
     636    psFree(yF32);
    427637
    428638    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,

  • trunk/psLib/src/sys/psErrorCodes.c

    r4556 r4568  
    11/** @file  psErrorCodes.c
    2  *
    3  *  @brief Contains the error codes for the error classes
     2*  @brief Contains the error codes for the error classes
    43 *
    54 *  @ingroup ErrorHandling
    … …  
    76 *  @author Robert DeSonia, MHPCC
    87 *
    9  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2005-07-15 02:33:54 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-07-16 00:06:33 $
     10=======
     11 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-07-16 00:06:33 $
    1113 *
    1214 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
Note: See TracChangeset for help on using the changeset viewer.