import ROOT
import array
import math

execfile('timeSeries_functions.py') # Load function definitions
execfile('timeSeries_functions_sol.py') #sol

# You will analyse pictures, which I obtained by filming a small
# slime mould plasmodium over a few hours. The pictures you will see
# have already been processed in various ways and only contain pixels
# with the values 0 and 1, where 1 means that in the original microscope
# picture was part of the plamodium.
# Slime moulds contract in a cyclic manner to move stuff around inside
# themselves (a slime mould consists of one huge cell). We want
# to measure the typical duration of a contraction-expansion cycle.

# TASK 1: Importing the data & understand the data format.

# Use the ImportJSON function to store the content of the data file in a variable.
# The filename is 'images.json.gz' and it contains a list of images.
# Each image consists of a list of list of int (a 2D list because we have
# a 2D coordinate system). Note that you can specify how many pictures are
# imported (optional) in the ImportJSON function, to save time (or prevent
# out of memory crashes).

#YOUR CODE HERE
images, times = ImportJSON('images.json.gz',200) #sol

# Draw the third image in the images list, using the DrawImage function

#YOUR CODE HERE
DrawImage(images[2]) #sol

# To check if you understood the data format (this will be important later),
# change the value of the pixel with the coordinates (50,100) to 1 in the
# picture you have drawn above and draw it again.

#YOUR CODE HERE
images[2][50][100]=1 #sol
DrawImage(images[2]) #sol

# Have a look at the whole dataset by calling the AnimateImages function
# with the whole list of pictures (You can change the speed if you want)

#YOUR CODE HERE
#AnimateImages(images) #sol

# TASK 2: Compute cell cross-section area. Plot an area time series.

# Fill in the missing code for the SumOfPixelValues function. Test it by
# printing the sum of values in picture 3.
# Print the total number of pixels (not values) of picture 3. Check if
# the fraction of the two numbers make sense (compare with area in drawings).
# Hint: Use the len function for finding the dimensions of the picture.
# Feel free to get rid of the animation in the code above (to save time).

#YOUR CODE HERE
summ=SumOfPixelValues(images[2]) #sol
tot=len(images[2])*len(images[2][0]) #sol
print "Sum of pixel values:",summ #sol
print "Total number of pixels",tot #sol
print "Fraction:", summ/tot #sol
print "\n" #sol

# We can use the SumOfPixels function to obtain an estimate of the
# Area of the cross-section of the cell in each picture.
# Compute this area for each picture and store the results in
# a list.

#YOUR CODE HERE
areas = [] #sol
for img in images: #sol
	areas.append(SumOfPixelValues(img)) #sol
canvas = DrawGraph(times,areas,title="Cell cross-section",xlabel="Time [min]",ylabel="Area") #sol

# TASK 3: Separate slow cell area changes (cell-growth) and fast
#         ones (contractions, expansions) using a moving average and FFT.

# Fill in the missing code for the MovingAverage function. Test it with a window
# of 22 data points in each direction on the time series of the area (45 in total).
# Note that the resulting time series should have the same length as the orignal
# one and should not display any artifacts at beginning and end.

#YOUR CODE HERE
MAareas = MovingAverage(areas, 22) #sol
AddGraph(times,MAareas) #sol

# Use the DFT function to calculate the fourier coefficients
# of the area time series. Use the Amplitude function to
# transform the coefficients into amplitudes.
# You can use the Frequencies function to generate a
# list of frequencies that matches the coefficients
# and amplitudes.
# Plot the amplitudes over frequencies (you can use
# the Log function for rescaling the amplitudes).

#YOUR CODE HERE
coefficients = DFT(areas) #sol
amplitudes = Amplitudes(coefficients) #sol
frequencies = Frequencies(len(amplitudes),times[1]) #sol
DrawGraph(frequencies,Log(amplitudes),ylabel="Log-Amplitudes",xlabel="Frequency [1/min]") #sol

# Compare the log-amplitude-frequency plot with the
# area time series. Which frequency range corresponds to
# the contractions and which part to cell growth? Find
# an appropriate threshold (it may help make a plot of
# a small part of the time series for easier frequency estimation).
# Use the LowPassFilter function to set all coefficients
# with lower frequency than the threshold to zero.
# Then transform the filtered coefficients back into a
# time series using the InverseDFT. Plot the resulting time
# series together with the original time series of the area.

#YOUR CODE HERE
DrawGraph(times[0:50],areas[0:50],title="Zoomed",ylabel="Area",xlabel="Time [min]") #sol
LPcoefficients = LowPassFilter(frequencies,coefficients,0.5) #sol
LPareas = InverseDFT(LPcoefficients) #sol
DrawGraph(times,areas,xlabel="Time [min]",ylabel="Area",title="Low Pass") #sol
AddGraph(times,LPareas) #sol
AddGraph(times,MAareas) #sol

# If you carefully tune your threshold, the low-pass
# filtered area curve should reflect the slow trends
# better than the moving average. The moving average
# always shows a small bump for each contraction cycle,
# which is undesirable. The low pass filtered curve
# also reflects quick changes in area better that are
# only slightly above the timescale of individual
# contraction cycles. Compare the moving average plot
# and the low pass plot to find these features. Adjust
# your threshold if necessary.
# Can you identify disadvantages of the low pass filter
# method, compared to the moving average?

# TASK 4: Use an autocorrelation to fit the residuals
#         of the smoothened time series.

# Calculate the difference of the original
# time series and one of the two
# smoothened time series (low pass or moving average).
# Use the Residuals function. Plot the result over
# time.

#YOUR CODE HERE
MAresiduals = Residuals(areas,MAareas) #sol
LPresiduals = Residuals(areas,LPareas) #sol
DrawGraph(times,MAresiduals,title="Residuals MA",xlabel="Time [min]",ylabel="Area") #sol
DrawGraph(times,LPresiduals,title="Residuals LP",xlabel="Time [min]",ylabel="Area") #sol

# Fill in the missing two lines in the autocorrelation function.
# Calculate and plot the autocorrelation of the residuals.
# Hint: Use the same list of times you use for your time
# series plots as x-values in the autocorrelation plot.
# 2nd Hint: Limit the autocorrelation plot to the first
# 10min (200 data points). Then you can easily see the
# length of a single contraction cycle.

#YOUR CODE HERE
MAacorr = Autocorrelation(MAresiduals) #sol
DrawGraph(times[0:200],MAacorr[0:200],title="Autocorr. MA",xlabel="Lag [min]",ylabel="Corr. (not norm.)") #sol
LPacorr = Autocorrelation(LPresiduals) #sol
DrawGraph(times[0:200],LPacorr[0:200],title="Autocorr. LP",xlabel="Lag [min]",ylabel="Corr. (not norm.)") #sol

# My lab report states an estimated contraction period
# of ~62.5 seconds. Does that fit with your own estimate?

raw_input("Press any key to exit")