|
| 1 | +package mascot.dynamics; |
| 2 | + |
| 3 | +import beast.base.spec.domain.Real; |
| 4 | +import beast.base.spec.inference.parameter.BoolVectorParam; |
| 5 | +import beast.base.spec.inference.parameter.RealScalarParam; |
| 6 | +import beast.base.spec.inference.parameter.RealVectorParam; |
| 7 | +import mascot.glmmodel.Covariate; |
| 8 | +import mascot.glmmodel.CovariateList; |
| 9 | +import mascot.glmmodel.LogLinear; |
| 10 | +import org.junit.jupiter.api.Test; |
| 11 | + |
| 12 | +import static org.junit.jupiter.api.Assertions.assertEquals; |
| 13 | +import static org.junit.jupiter.api.Assertions.assertTrue; |
| 14 | + |
| 15 | +public class GLMTest { |
| 16 | + |
| 17 | + /** |
| 18 | + * Regression for the single-epoch GLM intervalNr bug. |
| 19 | + * <p> |
| 20 | + * In a single-epoch model (rateShifts = [Infinity]), firstlargerzero == 0 |
| 21 | + * and dimension == 1, so the boundary branch in getCoalescentRate / getNe / |
| 22 | + * getMig is always taken for any i. The old hardcoded fallback |
| 23 | + * intervalNr = dim - 2 evaluated to -1 for dim=1, causing |
| 24 | + * ArrayIndexOutOfBoundsException when rates were looked up. The fix |
| 25 | + * uses dim - firstlargerzero - 1, which is the correct last valid |
| 26 | + * interval index for any rateShifts configuration. |
| 27 | + */ |
| 28 | + @Test |
| 29 | + public void testSingleEpochGLMRatesDoNotThrow() { |
| 30 | + int dim = 2; |
| 31 | + |
| 32 | + // single-epoch rate shifts: dim = 1 with the only value > 0, so firstlargerzero = 0 |
| 33 | + RateShifts rateShifts = new RateShifts(); |
| 34 | + rateShifts.initByName("value", "1.0"); |
| 35 | + |
| 36 | + GLM glm = buildGLM(dim, rateShifts); |
| 37 | + |
| 38 | + double[] coalRate = glm.getCoalescentRate(0); |
| 39 | + assertEquals(dim, coalRate.length); |
| 40 | + for (double v : coalRate) |
| 41 | + assertTrue(Double.isFinite(v), "coalescent rate must be finite, got " + v); |
| 42 | + |
| 43 | + for (int s = 0; s < dim; s++) { |
| 44 | + double ne = glm.getNe(s, 0); |
| 45 | + assertTrue(Double.isFinite(ne), "Ne must be finite, got " + ne); |
| 46 | + } |
| 47 | + |
| 48 | + double mig = glm.getMig(0, 1, 0); |
| 49 | + assertTrue(Double.isFinite(mig), "migration rate must be finite, got " + mig); |
| 50 | + } |
| 51 | + |
| 52 | + private GLM buildGLM(int dim, RateShifts rateShifts) { |
| 53 | + LogLinear migGLM = buildLogLinear(dim * (dim - 1)); |
| 54 | + LogLinear neGLM = buildLogLinear(dim); |
| 55 | + |
| 56 | + GLM glm = new GLM(); |
| 57 | + glm.initByName( |
| 58 | + "dimension", dim, |
| 59 | + "rateShifts", rateShifts, |
| 60 | + "migrationGLM", migGLM, |
| 61 | + "NeGLM", neGLM, |
| 62 | + "types", "a b"); |
| 63 | + return glm; |
| 64 | + } |
| 65 | + |
| 66 | + private LogLinear buildLogLinear(int covariateDim) { |
| 67 | + Double[] vals = new Double[covariateDim]; |
| 68 | + for (int i = 0; i < covariateDim; i++) |
| 69 | + vals[i] = 1.0; |
| 70 | + |
| 71 | + Covariate cov = new Covariate(vals, "cov"); |
| 72 | + cov.initAndValidate(); |
| 73 | + |
| 74 | + CovariateList covList = new CovariateList(); |
| 75 | + covList.initByName("covariates", cov); |
| 76 | + |
| 77 | + RealVectorParam<Real> scaler = new RealVectorParam<>(new double[]{0.0}, Real.INSTANCE); |
| 78 | + BoolVectorParam indicator = new BoolVectorParam(new boolean[]{true}); |
| 79 | + RealScalarParam<Real> clock = new RealScalarParam<>(); |
| 80 | + clock.initByName("value", "1.0"); |
| 81 | + |
| 82 | + LogLinear glm = new LogLinear(); |
| 83 | + glm.initByName( |
| 84 | + "covariateList", covList, |
| 85 | + "scaler", scaler, |
| 86 | + "indicator", indicator, |
| 87 | + "clock", clock); |
| 88 | + return glm; |
| 89 | + } |
| 90 | +} |
0 commit comments