-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP: Sub DeletePoint
[The following is a generalised function, used and reused by SQUID 2.50 on at least a couple of different occasions during a 'typical' data-reduction. As such, the inputs and outputs of the functions have generalised names. The identities of the specific inputs and outputs are defined when the function is called.]
The function performs the deletion of a single point from a vector of data-points, errors, and error-correlations. For the sigma-rho matrix, this essentially amounts to the deletion of the row and column of the matrix corresponding to the point to be deleted.
DeletePoint(LinReg, N, y1, y2, SigRho1, SigRho2, RejPoint, x1, x2)
LinReg: Boolean defining the type of calculation to be performed. TRUE = error-weighted linear regression vs time, with the regression-error calculated at analysis midtime; FALSE = time-independent error-weighted average and its error.
N: Integer defining the number of data-points submitted to the calculation.
y1: Vector array of length N, containing the y-values of the data-points submitted to the calculation.
y2: Vector array of length (N - 1), containing the y-values of the data-points returned by the calculation.
SigRho1: Array of size N x N, defining the sigma-rho matrix associated with the (y1) data-points submitted to the calculation. The diagonal elements are the sigma values associated with each Y value, and the (symmetric) immediate off-diagonal elements are the correlation coefficients between time-adjacent Y values. All other elements are zero.
SigRho2: Array of size (N - 1) x (N - 1), defining the sigma-rho matrix associated with the (y2) data-points returned by the calculation. The diagonal elements are the sigma values associated with each Y value, and the (symmetric) immediate off-diagonal elements are the correlation coefficients between time-adjacent Y values. All other elements are zero.
RejPoint: Integer representing the index of the point to be rejected.
[These apply ONLY when LinReg = True]
x1: Vector array of length N, containing the x-values of the data-points submitted to the calculation.
x2: Vector array of length (N - 1), containing the x-values of the data-points returned by the calculation.
Values of type Boolean
LinReg
Values of type Integer
j, m, Nn, p, RejPoint
The first step is to redefine the shortened array SigRho2 (and, if required, the shortened vector x2). Note that there does not seem to be any formal shortening/redefinition of the vector y2; as yet, I can't tell whether this is merely poor housekeeping, or deliberate/intentional.
Nn = N - 1
--Redefine as an array of size Nn x Nn (addressed [1 to Nn, 1 to Nn] ): SigRho2
If LinReg = True
--Redefine as a vector of length Nn (addressed 1 to Nn): x2
End If
The remainder of the subroutine defines vectors and arrays (x2, y2 and SigRho2) shortened by the omission of the deleted point, relative to unshortened 'reference' vectors and arrays (x1, y1 and SigRho1). Within the 'For... Next', the first If defines the values of x2, y2 and the diagonal elements of SigRho2, and the second If defines the off-diagonal elements of SigRho2. The latter accounts for the fact that when RejPoint is NOT 1 or N, the omission of RejPoint juxtaposes data-points that are no longer index-adjacent for the purpose of evaluating correlated errors, and for which the correlation coefficients need to be zero.
For j = 1 To N
m = j + 1
p = j + 2
If j < RejPoint
SigRho2[j, j] = SigRho1[j, j]
y2[j] = y1[j]
If LinReg = True
x2[j] = x1[j]
End If
ElseIf j < N
SigRho2[j, j] = SigRho1[m, m]
y2[j] = y1[m]
If LinReg = True
x2[j] = x1[m]
End If
End If
If j < (RejPoint - 1)
SigRho2[j, m] = SigRho1[j, m]
SigRho2[m, j] = SigRho1[m, j]
ElseIf j = (RejPoint - 1) And m < N
SigRho2[j, m] = 0
SigRho2[m, j] = 0
ElseIf j < (N - 1)
SigRho2[j, m] = SigRho1[m, p]
SigRho2[m, j] = SigRho1[p, m]
End If
Next j
End Sub