PeriodicityAnalysis.R

# The kf parameter of functions below should be a frame with
# Time and n_flux variables, like data from http://wherestheflux.com/

# amplitudePhaseAnalysis
# Produces amplitude and phase series by fitting sine/cosine 
# functions to a light curve.
amplitudePhaseAnalysis = function(kf, windowSize = 700, stepSize = 21, period = 0.88) {
    sequence = seq(from = 1, to = nrow(kf) - windowSize + 1, by = stepSize);
    lenseq = length(sequence);
    amplitudes = numeric(lenseq);
    phases = numeric(lenseq);
    times = numeric(lenseq);
    for(seqIdx in 1:lenseq) {
        start = sequence[seqIdx];
        range = start:(start + windowSize - 1);
        kfr = kf[range,];
        x = kfr$Time;
        y = kfr$n_flux;
        minX = min(x);
        maxX = max(x);
        sf = fitSinusoid(x, y, period);
        amplitudes[seqIdx] = sf$amp;
        phases[seqIdx] = sf$phase;
        times[seqIdx] = (minX + maxX) / 2;
    }    
    return(data.frame(
        Time = times,
        Amplitude = amplitudes,
        Phase = phases
    ));
}

periodAnalysis = function(kf, windowSize = 2000, stepSize = 200, minPeriod = 0.8, maxPeriod = 1.0, minAmpFactor = 2.0) {
    sequence = seq(from = 1, to = nrow(kf) - windowSize + 1, by = stepSize);
    lenseq = length(sequence);
    times = numeric(lenseq);
    periods = numeric(lenseq);
    for(seqIdx in 1:lenseq) {
        start = sequence[seqIdx];
        range = start:(start + windowSize - 1);
        kfr = kf[range,];
        x = kfr$Time;
        y = kfr$n_flux;
        minX = min(x);
        maxX = max(x);
        bestPeriod = periodSearch(x, y, minPeriod, maxPeriod);
        ampBest = fitSinusoid(x, y, period = bestPeriod)$amp;
        ampTwice = fitSinusoid(x, y, period = bestPeriod * 2)$amp;
        ampHalf = fitSinusoid(x, y, period = bestPeriod / 2)$amp;
        if(ampBest < ampTwice * minAmpFactor || ampBest < ampHalf * minAmpFactor) {
            bestPeriod = NA;            
        }
        times[seqIdx] = (minX + maxX) / 2;
        periods[seqIdx] = bestPeriod;
    }    
    return(data.frame(
        Time = times,
        Period = periods
    ));
}

minInfo = function(frame, fromTime, toTime, fieldName) {
    range = frame$Time >= fromTime & frame$Time <= toTime;
    paRange = frame[range,];
    paRangeField = paRange[[fieldName]];
    minAmplitude = min(paRangeField);
    minIndex = which(minAmplitude == paRangeField)[1];
    minTime = paRange[minIndex,]$Time;
    return(list(
       minAmp = minAmplitude,
       timeAtMin = minTime
    ));    
}

estimatePeriod = function(kf, fromTime, toTime, minPeriod = 0.6, maxPeriod = 1.0) {
    kfr = kf[kf$Time >= fromTime & kf$Time <= toTime,];
    return(optimizePeriod(kfr$Time, kfr$n_flux, minPeriod, maxPeriod));
}

optimizePeriod = function(x, y, minPeriod, maxPeriod) {
    bestPeriod = periodSearch(x, y, minPeriod, maxPeriod);
    sf = fitSinusoid(x, y, bestPeriod);
    return(list(
       period = bestPeriod,
       amp = sf$amp,
       phase = sf$phase
    ));    
}

periodSearch = function(x, y, minPeriod, maxPeriod, steps = 10, tolerance = 0.0001) {
    sequence = seq(from = minPeriod, to = maxPeriod, length.out = steps);
    bestPeriodIndex = 0;
    bestAmplitude = 0;
    bestPeriod = 0;
    lenseq = length(sequence);
    for(seqIdx in 1:lenseq) {
        period = sequence[[seqIdx]];
        sf = fitSinusoid(x, y, period);
        if(sf$amp > bestAmplitude) {
            bestAmplitude = sf$amp;
            bestPeriodIndex = seqIdx;
            bestPeriod = period;
        }
    }
    if(abs(maxPeriod - minPeriod) <= tolerance) {
        return(bestPeriod);
    }
    newFrom = ifelse(bestPeriodIndex <= 1, 1, bestPeriodIndex - 1);
    newTo = ifelse(bestPeriodIndex >= lenseq , lenseq, bestPeriodIndex + 1);
    return(periodSearch(x, y, sequence[newFrom], sequence[newTo], steps = steps, tolerance = tolerance));
}

produceSyntheticDip = function(fromTime = 120, toTime = 1580, interval = 0.024, period = 0.88, amplitude = 0.0001, noiseSD = 0.0001, pivot = 792.8, CLeft = 0.4, CRight = 0.2, plateauWidth = 0.20, minFlux = 0.85) {
    timeSeries = seq(from = fromTime, to = toTime, by = interval);
    leftBaseTime = pivot - plateauWidth / 2;
    rightBaseTime = pivot + plateauWidth / 2;
    declineSeriesLeft = CLeft / (CLeft + (leftBaseTime - timeSeries));
    declineSeriesRight = CRight / (CRight + (timeSeries - rightBaseTime));
    fluxFunction = function(time) {
        return(ifelse(time < leftBaseTime, 
            1 - declineSeriesLeft * (1 - minFlux),
            ifelse(time > rightBaseTime, 1 - declineSeriesRight * (1 - minFlux), minFlux)
        ));
    }; 
    baseFlux = fluxFunction(timeSeries);
    fluxWithSinusoid = baseFlux + (sin(timeSeries * 2 * pi / period) - 0.5) * amplitude;
    fluxWithNoise = fluxWithSinusoid + rnorm(n = length(fluxWithSinusoid)) * noiseSD;
    return(data.frame(
        Time = timeSeries,
        n_flux = fluxWithNoise
    ));    
}

produceSyntheticLightCurveWithDips = function(fromTime = 120, toTime = 1580, interval = 0.024, period = 0.88, amplitude = 0.0001, noiseSD = 0.0001, numDips = 20) {
    timeSeries = seq(from = fromTime, to = toTime, by = interval);
    fluxWithSinusoid = 1.0 + (sin(timeSeries * 2 * pi / period) - 0.5) * amplitude;
    fluxSeries = fluxWithSinusoid + rnorm(n = length(fluxWithSinusoid)) * noiseSD;
    dipTimes = sample(timeSeries, size = numDips, replace = TRUE);
    dipMinFlux = runif(n = numDips, min = 0.997, max = 0.999);
    dipSDs = runif(n = numDips, min = 0.5, max = 3.0);
    for(dipIndex in 1:numDips) {
        variance = dipSDs[dipIndex] ^ 2;
        fluxSeries = fluxSeries - (1 - dipMinFlux[dipIndex]) * exp(-((timeSeries - dipTimes[dipIndex])^2) / (2 * variance));
    }
    return(data.frame(
        Time = timeSeries,
        n_flux = fluxSeries
    ));    
}

testSinusoidFit = function(fromTime = 1000, toTime = 1100, interval = 0.024, period = 0.88, amplitude = 1.0, noiseSD = 0, trendFunction = NULL) {
    timeSeries = seq(from = fromTime, to = toTime, by = interval);
    sinusoid = (sin(timeSeries * 2 * pi / period) - 0.5) * amplitude;
    trendSeries = 0;
    if(!is.null(trendFunction)) {
        trendSeries = trendFunction(timeSeries);
    }
    y = sinusoid + rnorm(n = length(timeSeries)) * noiseSD + trendSeries;
    return(fitSinusoid(timeSeries, y, period));
}

# fitSinusoid
# This function fits a series given by (x,y) using a pair of sine/cosine
# functions and other monotonic functions, asuming the given period.
fitSinusoid = function(x, y, period) {
    angles = x * 2 * pi / period;
    minAngle = min(angles);
    maxAngle = max(angles);
    sinParams = sin(angles);
    cosParams = cos(angles);
    trendVar1 = angles;
    trendVar2 = (angles - minAngle) ^ 2;
    trendVar3 = log(angles - minAngle + 0.01);
    trendVar4 = log(maxAngle - angles + 0.01);
    model = lm(y ~ cbind(sinParams, cosParams, trendVar1, trendVar2, trendVar3, trendVar4));
    a = model$coefficients[[2]];
    b = model$coefficients[[3]];
    amp = sqrt(a^2 + b^2);
    phaseAngle = atan2(b, a);
    while(phaseAngle < 0) {
        phaseAngle = phaseAngle + 2 * pi;
    }
    phase = phaseAngle * period / (2 * pi);
    return(list(
        a = a,
        b = b,
        amp = amp,
        phase = phase
    ));    
}

displayFluxSeries = function(kf, fromTime, toTime, ylim = NULL, outPngFileName = NULL, title = NULL) {
    kfr = kf[kf$Time >= fromTime & kf$Time <= toTime,];
    par(cex=0.7, pch=19);
    plot(kfr$Time, kfr$n_flux, type="n", 
         xlab = "Day of Kepler Mission", 
         ylab = "Normalized Flux", 
         ylim = ylim);
    axisPoints = seq(fromTime, toTime, by=1);
    axis(1, at = axisPoints);
    abline(v = axisPoints, col="blue", lty="dashed");
    points(kfr$Time, kfr$n_flux, col="black");
    if(!is.null(title)) {
        title(main=title);        
    }
    if(!is.null(outPngFileName)) {
        dev.copy(png, paste0("./img/pa/", outPngFileName), width = 800, height = 400);
        dev.off();    
    }
}

displayAmplitudePhaseAnalysis = function(amplitudePhaseAnalysis, fromTime = 0, toTime= 1600, axisPoints = c(), ylim = c(0, 5E-4), outPngFileName = NULL, title = NULL, fieldName = "Amplitude") {
    paRange = amplitudePhaseAnalysis[amplitudePhaseAnalysis$Time >= fromTime & amplitudePhaseAnalysis$Time <= toTime,];
    par(cex=0.7, pch=19);
    plot(paRange$Time, paRange[[fieldName]], type="n", 
         xlab = "Day of Kepler Mission", 
         ylab = "Amplitude of signal",
         ylim = ylim);
    abline(v = axisPoints, col="blue", lty="dashed");
    points(paRange$Time, paRange[[fieldName]], col="black");
    if(!is.null(title)) {
        title(main=title);        
    }
    if(!is.null(outPngFileName)) {
        dev.copy(png, paste0("./img/pa/", outPngFileName), width = 800, height = 400);
        dev.off();    
    }
}