Before deciding what to do about this, I think we have to ask ourselves
what we are trying to achieve here. Will making these changes allow
anything to be done which cannot currently be done? I agree that it
*would* be more beautiful if there was no code duplication in the USSC,
but given that few of us have copious quantities of spare time, I think we
need to assign a priority to this work in competition with our other
tasks. I for one have been trying to get on with STC with little success
for one reason or another recently. And then there is making ast_resample
conserve flux which everyone seemed to see as a high priority (which I've
not even started on yet), and then there is the new FluxFrame class that
I was asked to produce (again no work so far)...
If this is a show stopper (or is likely to become one) then, fine, let's
get it sorted, but otherwise we should maybe live with it until we have
justified our continued existence to the forthcoming review.
David
On Mon, 4 Oct 2004, Tim Jenness wrote:
> Okay, hot on the heels of yesterday's "how come ccdpack has forked kappa
> routines" discussions here comes the rematch of "how come kappa has loads
> of ccdpack routines".
>
> KAPPA has 75 ccdpack routines, all but one are derived from generic source
> without actually having the generic source in CVS [ccdpack does]
>
> The diffs are below, the main difference seems to be that Dave has added
> extra IMETH switches to support for unweighted mean and also support for
> IMETH=300 and an alternative way of calculating variance by estimation.
>
> ccg1_mdr3 seems to be identical except for a missing bug fix comment from
> Mark (it seems that David fixed it during the same period).
>
> David's additions seem to be plausible for CCDPACK (a couple would
> require argument tweaks but nothing that couldn't be specified as a
> constant).
>
> --
> Tim Jenness
> JAC software
> http://www.jach.hawaii.edu/~timj
>
> Fortran file: ccd1_orva Generic: NONE
>
> --- ccg1_cm1d.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_cm1dd.f 2004-09-07 10:03:24.000000000 -1000
> @@ -41,7 +41,8 @@
> * IMETH = INTEGER (Given)
> * The method to use in combining the lines. Has a code of 2 to 10
> * which represent.
> -* 2 = MEAN
> +* 1 = UNWEIGHTED MEAN
> +* 2 = WEIGHTED MEAN
> * 3 = MEDIAN
> * 4 = TRIMMED MEAN
> * 5 = MODE
> @@ -50,6 +51,7 @@
> * 8 = MINMAX MEAN
> * 9 = BROADENED MEDIAN
> * 10 = SIGMA CLIPPED MEDIAN
> +* 300 = Median, but estimating variance from mean variance.
> * MINPIX = INTEGER (Given)
> * The minimum number of pixels required to contribute to an
> * output pixel.
> @@ -75,7 +77,7 @@
> * Workspace for order statistics calculations.
> * COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given and Returned)
> * Workspace for storing ordered statistics variace-covariance
> -* matrix.
> +* matrix. Not used for IMETHs 1, 2 or 300.
> * NCON( NLINES ) = DOUBLE PRECISION (Given and Returned)
> * The actual number of contributing pixels from each input line
> * to the output line.
> @@ -123,6 +125,14 @@
> * Added MEDCLIP method.
> * 16-NOV-1998 (PDRAPER):
> * Added fast median.
> +* 14-DEC-2001 (DSB):
> +* Added unweighted mean method.
> +* 1-NOV-2002 (DSB):
> +* Added workspace argument to CCD1_ORVAR.
> +* 27-AUG-2003 (DSB):
> +* Added IMETH 300 to provide a means of estimating the output
> +* variances for large samples, for which the memory requirements
> +* for the COVEC array would be too large.
> * {enter_further_changes_here}
>
> * Bugs:
> @@ -135,6 +145,7 @@
>
> * Global Constants:
> INCLUDE 'SAE_PAR' ! Standard SAE constants
> + INCLUDE 'CNF_PAR' ! For CNF_PVAL function
>
> * Arguments Given:
> INTEGER NPIX
> @@ -166,6 +177,8 @@
> * Status:
> INTEGER STATUS ! Global status
>
> +* Local Variables:
> + INTEGER IPW1 ! Work space pointer
> *.
>
> * Check inherited global status.
> @@ -175,23 +188,33 @@
> * for the order statistics of a normal popl with up to NLINE members.
> * This also sets up the scale factor for converting mean variances to
> * median variances.
> - IF( IMETH .NE. 2 ) THEN
> - CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC, STATUS )
> + IF( IMETH .NE. 1 .AND. IMETH .NE. 2 .AND. IMETH .NE. 300 ) THEN
> + CALL PSX_CALLOC( NLINES*NLINES, '_DOUBLE', IPW1, STATUS )
> + CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC,
> + : %VAL( CNF_PVAL( IPW1 ) ),
> + : STATUS )
> + CALL PSX_FREE( IPW1, STATUS )
> END IF
>
> * Now branch for each method.
> - IF ( IMETH .EQ. 2 ) THEN
> + IF ( IMETH .EQ. 1 ) THEN
> +
> +* Forming the unweighted mean.
> + CALL CCG1_UMD1D( STACK, NPIX, NLINES, VARS, MINPIX,
> + : RESULT, RESVAR, NCON, STATUS )
> +
> + ELSE IF ( IMETH .EQ. 2 ) THEN
>
> * Forming the weighted mean.
> CALL CCG1_MED1D( STACK, NPIX, NLINES, VARS, MINPIX,
> : RESULT, RESVAR, NCON, STATUS )
>
> - ELSE IF ( IMETH .EQ. 3 ) THEN
> + ELSE IF ( IMETH .EQ. 3 .OR. IMETH .EQ. 300 ) THEN
>
> * Forming the weighted median.
> - CALL CCG1_MDD1D( STACK, NPIX, NLINES, VARS, MINPIX, COVEC,
> - : NMAT, RESULT, RESVAR, WRK1, WRK2, NCON,
> - : POINT, USED, STATUS )
> + CALL CCG1_MDD1D( ( IMETH .EQ. 3 ), STACK, NPIX, NLINES, VARS,
> + : MINPIX, COVEC, NMAT, RESULT, RESVAR, WRK1,
> + : WRK2, NCON, POINT, USED, STATUS )
>
> ELSE IF ( IMETH .EQ. 4 ) THEN
>
>
> Fortran file: ccg1_cm1r Generic: ccg1_cm1r.gen
> --- ccg1_cm1r.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_cm1rr.f 2004-09-07 10:03:24.000000000 -1000
> @@ -41,7 +41,8 @@
> * IMETH = INTEGER (Given)
> * The method to use in combining the lines. Has a code of 2 to 9
> * which represent.
> -* 2 = MEAN
> +* 1 = MEAN
> +* 2 = WEIGHTED MEAN
> * 3 = MEDIAN
> * 4 = TRIMMED MEAN
> * 5 = MODE
> @@ -50,6 +51,7 @@
> * 8 = MINMAX MEAN
> * 9 = BROADENED MEDIAN
> * 10 = SIGMA CLIPPED MEDIAN
> +* 300 = Median, but estimating variance from mean variance.
> * MINPIX = INTEGER (Given)
> * The minimum number of pixels required to contribute to an
> * output pixel.
> @@ -75,7 +77,7 @@
> * Workspace for order statistics calculations.
> * COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given and Returned)
> * Workspace for storing ordered statistics variace-covariance
> -* matrix.
> +* matrix. Not used for IMETHs 1, 2 or 300.
> * NCON( NLINES ) = DOUBLE PRECISION (Given and Returned)
> * The actual number of contributing pixels from each input line
> * to the output line.
> @@ -123,6 +125,14 @@
> * Added sigma clipped median.
> * 16-NOV-1998 (PDRAPER):
> * Added fast median.
> +* 9-SEP-2002 (DSB):
> +* Added unweighted mean method.
> +* 1-NOV-2002 (DSB):
> +* Added workspace argument to CCD1_ORVAR.
> +* 27-AUG-2003 (DSB):
> +* Added IMETH 300 to provide a means of estimating the output
> +* variances for large samples, for which the memory requirements
> +* for the COVEC array would be too large.
> * {enter_further_changes_here}
>
> * Bugs:
> @@ -135,6 +145,7 @@
>
> * Global Constants:
> INCLUDE 'SAE_PAR' ! Standard SAE constants
> + INCLUDE 'CNF_PAR' ! For CNF_PVAL function
>
> * Arguments Given:
> INTEGER NPIX
> @@ -166,6 +177,8 @@
> * Status:
> INTEGER STATUS ! Global status
>
> +* Local Variables:
> + INTEGER IPW1 ! Work space pointer
> *.
>
> * Check inherited global status.
> @@ -175,23 +188,33 @@
> * for the order statistics of a normal popl with up to NLINE members.
> * This also sets up the scale factor for converting mean variances to
> * median variances.
> - IF( IMETH .NE. 2 ) THEN
> - CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC, STATUS )
> + IF( IMETH .NE. 1 .AND. IMETH .NE. 2 .AND. IMETH .NE. 300 ) THEN
> + CALL PSX_CALLOC( NLINES*NLINES, '_DOUBLE', IPW1, STATUS )
> + CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC,
> + : %VAL( CNF_PVAL( IPW1 ) ),
> + : STATUS )
> + CALL PSX_FREE( IPW1, STATUS )
> END IF
>
> * Now branch for each method.
> - IF ( IMETH .EQ. 2 ) THEN
> + IF ( IMETH .EQ. 1 ) THEN
> +
> +* Forming the unweighted mean.
> + CALL CCG1_UMR1R( STACK, NPIX, NLINES, VARS, MINPIX,
> + : RESULT, RESVAR, NCON, STATUS )
> +
> + ELSE IF ( IMETH .EQ. 2 ) THEN
>
> * Forming the weighted mean.
> CALL CCG1_MER1R( STACK, NPIX, NLINES, VARS, MINPIX,
> : RESULT, RESVAR, NCON, STATUS )
>
> - ELSE IF ( IMETH .EQ. 3 ) THEN
> + ELSE IF ( IMETH .EQ. 3 .OR. IMETH .EQ. 300 ) THEN
>
> * Forming the weighted median.
> - CALL CCG1_MDR1R( STACK, NPIX, NLINES, VARS, MINPIX, COVEC,
> - : NMAT, RESULT, RESVAR, WRK1, WRK2, NCON,
> - : POINT, USED, STATUS )
> + CALL CCG1_MDR1R( ( IMETH .EQ. 3 ), STACK, NPIX, NLINES, VARS,
> + : MINPIX, COVEC, NMAT, RESULT, RESVAR, WRK1,
> + : WRK2, NCON, POINT, USED, STATUS )
>
> ELSE IF ( IMETH .EQ. 4 ) THEN
>
>
> Fortran file: ccg1_cm3d Generic: ccg1_cm3d.gen
> --- ccg1_cm3d.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_cm3dd.f 2002-09-10 00:12:55.000000000 -1000
> @@ -39,7 +39,8 @@
> * IMETH = INTEGER (Given)
> * The method to use in combining the lines. Has a code of 2 to 10.
> * which represent.
> -* 2 = MEAN
> +* 1 = MEAN
> +* 2 = WEIGHTED MEAN
> * 3 = MEDIAN
> * 4 = TRIMMED MEAN
> * 5 = MODE
> @@ -100,6 +101,8 @@
> * Added sigma clipped median.
> * 16-NOV-1998 (PDRAPER):
> * Added fast median.
> +* 9-SEP-2002 (DSB):
> +* Added unweighted mean method.
> * {enter_further_changes_here}
>
> * Bugs:
> @@ -145,7 +148,13 @@
> IF ( STATUS .NE. SAI__OK ) RETURN
>
> * Branch for each method.
> - IF ( IMETH .EQ. 2 ) THEN
> + IF ( IMETH .EQ. 1 ) THEN
> +
> +* Forming the unweighted mean.
> + CALL CCG1_UMD3D( STACK, NPIX, NLINES, VARS, MINPIX,
> + : RESULT, NCON, STATUS )
> +
> + ELSE IF ( IMETH .EQ. 2 ) THEN
>
> * Forming the weighted mean.
> CALL CCG1_MED3D( STACK, NPIX, NLINES, VARS, MINPIX,
> @@ -220,4 +229,4 @@
> END IF
>
> END
> -* $Id: ccg1_cm3d.gen,v 1.3 1998/11/18 19:12:48 pdraper Exp $
> +* $Id: ccg1_cm3dd.f,v 1.2 2002/09/10 10:12:55 dsb Exp $
> Fortran file: ccg1_cm3r Generic: ccg1_cm3r.gen
> --- ccg1_cm3r.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_cm3rr.f 2002-09-10 00:12:55.000000000 -1000
> @@ -39,7 +39,8 @@
> * IMETH = INTEGER (Given)
> * The method to use in combining the lines. Has a code of 2 to 10.
> * which represent.
> -* 2 = MEAN
> +* 1 = MEAN
> +* 2 = WEIGHTED MEAN
> * 3 = MEDIAN
> * 4 = TRIMMED MEAN
> * 5 = MODE
> @@ -101,6 +102,8 @@
> * Added sigma clipped median.
> * 16-NOV-1998 (PDRAPER):
> * Added fast median.
> +* 9-SEP-2002 (DSB):
> +* Added unweighted mean method.
> * {enter_further_changes_here}
>
> * Bugs:
> @@ -146,7 +149,13 @@
> IF ( STATUS .NE. SAI__OK ) RETURN
>
> * Branch for each method.
> - IF ( IMETH .EQ. 2 ) THEN
> + IF ( IMETH .EQ. 1 ) THEN
> +
> +* Forming the unweighted mean.
> + CALL CCG1_UMR3R( STACK, NPIX, NLINES, VARS, MINPIX,
> + : RESULT, NCON, STATUS )
> +
> + ELSE IF ( IMETH .EQ. 2 ) THEN
>
> * Forming the weighted mean.
> CALL CCG1_MER3R( STACK, NPIX, NLINES, VARS, MINPIX,
>
> Fortran file: ccg1_mdd1 Generic: ccg1_mdd1.gen
> --- ccg1_mdd1.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_mdd1d.f 2003-08-27 08:20:06.000000000 -1000
> @@ -1,6 +1,6 @@
> - SUBROUTINE CCG1_MDD1D( STACK, NPIX, NLINES, VARS, MINPIX,
> - : COVEC, NMAT, RESULT, RESVAR, WRK1,
> - : WRK2, NCON, POINT, USED, STATUS )
> + SUBROUTINE CCG1_MDD1D( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
> + : COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
> + : NCON, POINT, USED, STATUS )
> *+
> * Name:
> * CCG1_MDD1
> @@ -12,9 +12,9 @@
> * Starlink Fortran 77
>
> * Invocation:
> -* CALL CCG1_MDD1( STACK, NPIX, NLINES, VARS, MINPIX, COVEC, NMAT,
> -* RESULT, RESVAR, WRK1, WRK2, NCON, POINT, USED,
> -* STATUS )
> +* CALL CCG1_MDD1D( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
> +* COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
> +* NCON, POINT, USED, STATUS )
>
> * Description:
> * This routine accepts an array consisting a series of (vectorised)
> @@ -24,6 +24,11 @@
> * in the array RESULT.
>
> * Arguments:
> +* CALCMV = LOGICAL (Given)
> +* If .FALSE. then the output variances are estimated by scaling
> +* the variance on the weighted mean (rather than the weighted median)
> +* by Pi/2. Otherwise, the output variances are calculated using the
> +* COVEC array.
> * STACK( NPIX, NLINES ) = DOUBLE PRECISION (Given)
> * The array of lines which are to be combined into a single line.
> * NPIX = INTEGER (Given)
> @@ -38,7 +43,7 @@
> * COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given)
> * The packed variance-covariance matrix of the order statistics
> * from a normal distribution of sizes up to NLINES, produced by
> -* CCD1_ORVAR.
> +* CCD1_ORVAR. Not used if CALCMV is .FALSE.
> * NMAT = INTEGER (Given)
> * Size of the first dimension of COVEC.
> * RESULT( NPIX ) = DOUBLE PRECISION (Returned)
> @@ -67,6 +72,8 @@
> * History:
> * 21-MAY-1992 (PDRAPER):
> * Original version.
> +* 27-AUG-2003 (DSB):
> +* Added argument CALCMV.
> * {enter_changes_here}
>
> * Bugs:
> @@ -82,6 +89,7 @@
> INCLUDE 'PRM_PAR' ! PRIMDAT constants
>
> * Arguments Given:
> + LOGICAL CALCMV
> INTEGER NPIX
> INTEGER NLINES
> INTEGER MINPIX
> @@ -178,8 +186,13 @@
> CALL CCG1_IS3D( WRK1, WRK2, POINT, NGOOD, STATUS )
>
> * Find the weighted median.
> - CALL CCG1_WTM3D( WRK1, WRK2, SVAR, NGOOD, USED,
> - : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
> + IF( CALCMV ) THEN
> + CALL CCG1_WTM3D( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
> + : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
> + ELSE
> + CALL CCG1_WTM3D( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
> + : COVEC, VAL, VAR, STATUS )
> + END IF
>
> * Update the used line counters.
> DO 3 J = 1, NGOOD
> @@ -211,4 +224,4 @@
> CALL NUM_REVRT
>
> END
> -* $Id: ccg1_mdd1.gen,v 1.3 1998/11/20 16:05:26 pdraper Exp $
> +* $Id: ccg1_mdd1d.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
> Fortran file: ccg1_mdd3 Generic: ccg1_mdd3.gen
> --- ccg1_mdd3.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_mdd3d.f 2000-09-05 03:02:58.000000000 -1000
> @@ -53,15 +53,11 @@
>
> * Authors:
> * PDRAPER: Peter Draper (STARLINK)
> -* MBT: Mark Taylor (STARLINK)
> * {enter_new_authors_here}
>
> * History:
> * 19-MAY-1992 (PDRAPER):
> * Original version.
> -* 27-NOV-2000 (MBT):
> -* Fixed bug; variance values were not sorted corresponding to
> -* order of data values being stacked.
> * {enter_changes_here}
>
> * Bugs:
> @@ -157,7 +153,7 @@
> IF ( NGOOD .GT. 0 ) THEN
>
> * Sort these values into increasing order.
> - CALL CCG1_IS3D( WRK1, WRK2, POINT, NGOOD, STATUS )
> + CALL CCG1_IS2D( WRK1, POINT, NGOOD, STATUS )
>
> * Find the weighted median.
> CALL CCG1_WTM2D( WRK1, WRK2, NGOOD, USED, VAL, STATUS )
> @@ -189,4 +185,4 @@
> CALL NUM_REVRT
>
> END
> -* $Id: ccg1_mdd3.gen,v 1.5 2000/11/28 13:12:14 mbt Exp $
> +* $Id: ccg1_mdd3d.f,v 1.1 2000/09/05 13:02:58 dsb Exp $
> Fortran file: ccg1_mdr1 Generic: ccg1_mdr1.gen
> --- ccg1_mdr1.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_mdr1r.f 2003-08-27 08:20:06.000000000 -1000
> @@ -1,6 +1,6 @@
> - SUBROUTINE CCG1_MDR1R( STACK, NPIX, NLINES, VARS, MINPIX,
> - : COVEC, NMAT, RESULT, RESVAR, WRK1,
> - : WRK2, NCON, POINT, USED, STATUS )
> + SUBROUTINE CCG1_MDR1R( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
> + : COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
> + : NCON, POINT, USED, STATUS )
> *+
> * Name:
> * CCG1_MDR1
> @@ -12,9 +12,9 @@
> * Starlink Fortran 77
>
> * Invocation:
> -* CALL CCG1_MDR1( STACK, NPIX, NLINES, VARS, MINPIX, COVEC, NMAT,
> -* RESULT, RESVAR, WRK1, WRK2, NCON, POINT, USED,
> -* STATUS )
> +* CALL CCG1_MDR1R( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
> +* COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
> +* NCON, POINT, USED, STATUS )
>
> * Description:
> * This routine accepts an array consisting a series of (vectorised)
> @@ -24,6 +24,11 @@
> * output means are returned in the array RESULT.
>
> * Arguments:
> +* CALCMV = LOGICAL (Given)
> +* If .FALSE. then the output variances are estimated by scaling
> +* the variance on the weighted mean (rather than the weighted median)
> +* by Pi/2. Otherwise, the output variances are calculated using the
> +* COVEC array.
> * STACK( NPIX, NLINES ) = REAL (Given)
> * The array of lines which are to be combined into a single line.
> * NPIX = INTEGER (Given)
> @@ -38,7 +43,7 @@
> * COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given)
> * The packed variance-covariance matrix of the order statistics
> * from a normal distribution of sizes up to NLINES, produced by
> -* CCD1_ORVAR.
> +* CCD1_ORVAR. Not used if CALCMV is .FALSE.
> * NMAT = INTEGER (Given)
> * Size of the first dimension of COVEC.
> * RESULT( NPIX ) = REAL (Returned)
> @@ -67,6 +72,8 @@
> * History:
> * 20-MAY-1992 (PDRAPER):
> * Original version.
> +* 27-AUG-2003 (DSB):
> +* Added argument CALCMV.
> * {enter_changes_here}
>
> * Bugs:
> @@ -82,6 +89,7 @@
> INCLUDE 'PRM_PAR' ! PRIMDAT constants
>
> * Arguments Given:
> + LOGICAL CALCMV
> INTEGER NPIX
> INTEGER NLINES
> INTEGER MINPIX
> @@ -178,8 +186,13 @@
> CALL CCG1_IS3R( WRK1, WRK2, POINT, NGOOD, STATUS )
>
> * Find the weighted median.
> - CALL CCG1_WTM3R( WRK1, WRK2, SVAR, NGOOD, USED,
> - : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
> + IF( CALCMV ) THEN
> + CALL CCG1_WTM3R( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
> + : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
> + ELSE
> + CALL CCG1_WTM3R( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
> + : COVEC, VAL, VAR, STATUS )
> + END IF
>
> * Update the used line counters.
> DO 3 J = 1, NGOOD
> @@ -211,4 +224,4 @@
> CALL NUM_REVRT
>
> END
> -* $Id: ccg1_mdr1.gen,v 1.3 1998/11/20 16:05:26 pdraper Exp $
> +* $Id: ccg1_mdr1r.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
> Fortran file: ccg1_mdr3 Generic: ccg1_mdr3.gen
> --- ccg1_mdr3.f 2004-10-04 16:25:59.000000000 -1000
> +++ ccg1_mdr3r.f 2001-11-13 03:32:18.000000000 -1000
> @@ -53,15 +53,11 @@
>
> * Authors:
> * PDRAPER: Peter Draper (STARLINK)
> -* MBT: Mark Taylor (STARLINK)
> * {enter_new_authors_here}
>
> * History:
> * 19-MAY-1992 (PDRAPER):
> * Original version.
> -* 27-NOV-2000 (MBT):
> -* Fixed bug; variance values were not sorted corresponding to
> -* order of data values being stacked.
> * {enter_changes_here}
>
> * Bugs:
> @@ -189,4 +185,4 @@
> CALL NUM_REVRT
>
> END
> -* $Id: ccg1_mdr3.gen,v 1.5 2000/11/28 13:12:14 mbt Exp $
> +* $Id: ccg1_mdr3r.f,v 1.2 2001/11/13 13:32:18 dsb Exp $
>
> --- ccg1_wtm3.f 2004-10-04 16:26:01.000000000 -1000
> +++ ccg1_wtm3d.f 2003-08-27 08:20:06.000000000 -1000
> @@ -1,5 +1,5 @@
> - SUBROUTINE CCG1_WTM3D( ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
> - : RESULT, RESVAR, STATUS )
> + SUBROUTINE CCG1_WTM3D( CALCMV, ORDDAT, WEIGHT, VAR, NENT, USED,
> + : COVAR, RESULT, RESVAR, STATUS )
> *+
> * Name:
> * CCG1_WTM3D
> @@ -13,8 +13,8 @@
> * Starlink Fortran 77
>
> * Invocation:
> -* CALL CCG1_WTM3D( ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
> -* RESULT, RESVAR, STATUS )
> +* CALL CCG1_WTM3D( CALCMV, ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
> +* RESULT, RESVAR, STATUS )
>
> * Description:
> * This routine finds a value which can be associated with the half-
> @@ -33,6 +33,11 @@
> * variance.
>
> * Arguments:
> +* CALCMV = LOGICAL (Given)
> +* If .FALSE. then the output variances are estimated by scaling
> +* the variance on the weighted mean (rather than the weighted median)
> +* by Pi/2. Otherwise, the output variances are calculated using the
> +* COVEC array.
> * ARR( NENT ) = DOUBLE PRECISION (Given)
> * The list of ordered data for which the weighted median is
> * required
> @@ -47,7 +52,7 @@
> * true in this array, otherwise the array is set to false.
> * COVAR( * ) = DOUBLE PRECISION (Given)
> * The packed variance-covariance matrix of the order statistics
> -* from a normal distribution of size NENT.
> +* from a normal distribution of size NENT. Not used if CALCMV is .FALSE.
> * RESULT = DOUBLE PRECISION (Returned)
> * The weighted median
> * RESVAR = DOUBLE PRECISION (Returned)
> @@ -77,6 +82,8 @@
> * Added used array.
> * 27-MAY-1992 (PDRAPER):
> * Added direct variance estimates.
> +* 27-AUG-2003 (DSB):
> +* Added argument CALCMV.
> * {enter_further_changes_here}
>
> * Bugs:
> @@ -91,6 +98,7 @@
> INCLUDE 'SAE_PAR' ! Standard SAE constants
>
> * Arguments Given:
> + LOGICAL CALCMV
> INTEGER NENT
> DOUBLE PRECISION ORDDAT( NENT )
> DOUBLE PRECISION WEIGHT( NENT )
> @@ -105,8 +113,13 @@
> * Status:
> INTEGER STATUS ! Global status
>
> +* Local Constants:
> + DOUBLE PRECISION PIBY2 ! 90 degs in radians
> + PARAMETER ( PIBY2 = 1.57079632679489661923132 )
> +
> * Local Variables:
> DOUBLE PRECISION TOTWT ! The total value of weights
> + DOUBLE PRECISION TRGWT ! The target value of weights
> INTEGER I ! Loop variable
> INTEGER J ! Loop variable
> INTEGER K ! Loop variable
> @@ -150,7 +163,7 @@
> 1 CONTINUE
>
> * Search for median weight.
> - TOTWT = TOTWT * 0.5D0
> + TRGWT = TOTWT * 0.5D0
> WTSUM = 0.0D0
> DO 2 I = 1, NENT
> IF ( I .EQ. 1 ) THEN
> @@ -160,7 +173,7 @@
> : * 0.5D0
> END IF
> WTSUM = WTSUM + WTINC
> - IF ( WTSUM .GT. TOTWT ) GO TO 66
> + IF ( WTSUM .GT. TRGWT ) GO TO 66
> 2 CONTINUE
> 66 I = MIN( I, NENT )
>
> @@ -177,43 +190,51 @@
> USED( I - 1 ) = .TRUE.
>
> * Set weights factors
> - W1 = ( WTSUM - TOTWT ) / MAX( WTINC, 1.0D-20 )
> + W1 = ( WTSUM - TRGWT ) / MAX( WTINC, 1.0D-20 )
> W2 = 1.0D0 - W1
>
> * Interpolate between data values
> RESULT = D1 * W1 + D2 * W2
>
> -* Sum the relevant ordered statistic variances and covariances.
> -* weighting accordingly.
> - VSUM = 0.0D0
> - DO 3 K = LBND, UBND
> - IF( K .EQ. LBND ) THEN
> - IW = W1
> - ELSE
> - IW = W2
> - END IF
> - DO 4 J = K, UBND
> - IF( J .EQ. LBND ) THEN
> - JW = W1
> +* If required, sum the relevant ordered statistic variances and
> +* covariances weighting accordingly.
> + IF( CALCMV ) THEN
> + VSUM = 0.0D0
> + DO 3 K = LBND, UBND
> + IF( K .EQ. LBND ) THEN
> + IW = W1
> ELSE
> - JW = W2
> + IW = W2
> END IF
> + DO 4 J = K, UBND
> + IF( J .EQ. LBND ) THEN
> + JW = W1
> + ELSE
> + JW = W2
> + END IF
>
> * Sum variances and twice covariances ( off diagonal elements ).
> - IF( K .EQ. J ) THEN
> - VSUM = VSUM + IW * JW * COVAR( K + J * ( J - 1 )/ 2 )
> - ELSE
> - VSUM = VSUM +
> - : 2.0D0 * IW * JW * COVAR( K + J * ( J - 1 )/ 2 )
> - END IF
> - 4 CONTINUE
> - 3 CONTINUE
> + IF( K .EQ. J ) THEN
> + VSUM = VSUM + IW * JW * COVAR( K + J*( J - 1 )/2 )
> + ELSE
> + VSUM = VSUM +
> + : 2.0D0 * IW * JW * COVAR( K + J*( J - 1 )/2 )
> + END IF
> + 4 CONTINUE
> + 3 CONTINUE
>
> * Right make the new variance estimate. Use the sum of variances
> * and covariances of the order statistic of the `trimmed' sample size
> * Sample variance changes to NENT * VAR to represent total variance
> * of original data.
> - RESVAR = VAR * NENT * VSUM
> + RESVAR = VAR * NENT * VSUM
> +
> +* If we are estimating the variance on the basis of the variance of the
> +* mean, calculate the estimate.
> + ELSE
> + RESVAR = PIBY2/TOTWT
> + END IF
> +
> END IF
> END
> -* $Id: ccg1_wtm3.gen,v 1.1 1997/06/27 09:01:41 pwd Exp $
> +* $Id: ccg1_wtm3d.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
>
|