Some confusions for the code in Peeling.py

[AprilYUZhang]

This function as follows didn't value a object, has it been discard?

estimateSegregationWithNorm(
                segregationTensor,
                segregationTensor_norm,
                parentsMinusChild[i, :, :, :],
                childValues,
                pointSeg[child, :, :],
            )

@jit(nopython=True, nogil=True)
def estimateSegregationWithNorm(
    segregationTensor, segregationTensor_norm, parentValues, childValues, output
):
    # pointSeg[child,:,:] = np.einsum("abcd, abi, ci-> di", segregationTensor, parentsMinusChild[i,:,:,:], childValues)
    nLoci = childValues.shape[1]
    output[:, :] = 0
    for a in range(4):
        for b in range(4):
            for c in range(4):
                for d in range(4):
                    for i in range(nLoci):
                        # Check if norm is 0. Otherwise use norm to normalize.
                        if segregationTensor_norm[a, b, c] != 0:
                            output[d, i] += (
                                segregationTensor[a, b, c, d]
                                * parentValues[a, b, i]
                                * childValues[c, i]
                                / segregationTensor_norm[a, b, c]
                            )
    return output

[AprilYUZhang]

The same issue:

createChildSegs(
            segregationTensor, currentSeg, childSegTensor[index, :, :, :, :]
        )

def createChildSegs(segregationTensor, currentSeg, output):
    # childSegs[index,:,:,:,:] = np.einsum("abcd, di -> abci", segregationTensor, currentSeg)
    nLoci = currentSeg.shape[1]
    output[:, :, :, :] = 0
    for a in range(4):
        for b in range(4):
            for c in range(4):
                for d in range(4):
                    for i in range(nLoci):
                        output[a, b, c, i] += (
                            segregationTensor[a, b, c, d] * currentSeg[d, i]
                        )

    return output

[AprilYUZhang]

 projectChildGenotypes(
            childSegTensor[index, :, :, :, :],
            childValues,
            childToParents[index, :, :, :],
        )

def projectChildGenotypes(childSegs, childValues, output):
    # childToParents[index,:,:,:] = np.einsum("abci, ci -> abi", childSegs[index,:,:,:,:], childValues)
    nLoci = childSegs.shape[3]
    output[:, :, :] = 0
    for a in range(4):
        for b in range(4):
            for c in range(4):
                for i in range(nLoci):
                    output[a, b, i] += childSegs[a, b, c, i] * childValues[c, i]

    return output

[AprilYUZhang]

projectParentGenotypes(
                childSegTensor[i, :, :, :, :],
                parentsMinusChild[i, :, :, :],
                anterior[child, :, :],
            )

@jit(nopython=True, nogil=True)
def projectParentGenotypes(childSegs, parentValues, output):
    # anterior[child,:,:] = np.einsum("abci, abi -> ci", childSegs[i,:,:,:,:], parentsMinusChild[i,:,:,:])
    nLoci = childSegs.shape[3]
    output[:, :] = 0
    for a in range(4):
        for b in range(4):
            for c in range(4):
                for i in range(nLoci):
                    output[c, i] += childSegs[a, b, c, i] * parentValues[a, b, i]

    return output

[XingerTang]

@AprilYUZhang Arrays are mutable objects which means that their values are changed even when they are not returned as the outputs of the functions.

[AprilYUZhang]

It still makes me confused because these function just create a new array and there is no object to accept this array rather than changing the original arrays.

[XingerTang]

It still makes me confused because these function just create a new array and there is no object to accept this array rather than changing the original arrays.

Actually, no new array would be created when the function is called, and this is a result of how mutable objects are constructed and stored. Here is a good reference for the mutable objects: https://realpython.com/python-mutable-vs-immutable-types/

[AprilYUZhang]

Thank you; I got it. "return output" is useless, right?

[XingerTang]

Thank you; I got it. "return output" is useless, right?

:) No return is needed.