#========================================================================================#
#	Laboratory 9
#
# Working with noisy data.
#
# Author: Niall Palfreyman, 17/03/2022
#========================================================================================#
[
	Activity(
		"""
		GAs showed us how useful random processes can be for finding the solution to problems.
		Indeed, the entire field of biology is concerned with the way living organisms use
		random processes to solve problems. Notice this does not mean biological process are
		random: the breakdown of starch by ptyalin in your mouth follows a random path, yet it
		is absolutely certain that the starch will be broken down. However, we think of organisms
		as CHOOSING what to do. That is, they respond creatively to conditions by guiding their
		own random processes into EXPLORATORY processes that are neither entirely random nor
		entirely deterministic.

		To see how they do this, consider the bacterial virus lambda-phage (λ-phage). When
		λ-phage infects an E.coli bacterium, the infection can go in either of two different
		ways: lysis (rupturing the bacterium) or lysogeny (virus becomes integrated into the
		bacterium's DNA). Now comes the interesting part: Observed in a laboratory, it seems as
		if the lysis/lysogeny fate of the bacterium is purely random, driven by statistical
		variations in the concentration of certain transcription factors in the cell. This is
		very different from non-biochemical reactions, which either happen or they don't!

		Let's do some quick calculations. The E.coli cell is a cylinder with length 2μm and
		diameter 1μm. Use Julia to calculate the volume of this cell in litres:
		""",
		"Remember the formula for cylinder volume, and use the float notation 1e-6",
		x -> (1e-12 < x < 2e-12)
	),
	Activity(
		"""
		Transcription factors influnce dynamics in the cell at concentrations of around 1nM
		(i.e., nanomoles per litre). Use this value to calculate the total number of
		transcription factor molecules inside the E.coli cell:
		""",
		"Remember that you will need the Avogadro constant: 6e23 particles per mole",
		x -> (900 < x < 1000)
	),
	Activity(
		"""
		As you can see, there are only under 1000 molecules of the transcription factor inside the
		E.coli cell, so statistical variations play an important role in determining the
		decisions that the cell makes. Organisms form themselves in sizes and ways that can
		make use of random noise to explore multiple possible solutions.

		To model exploratory processes on a computer, we use PSEUDORANDOM numbers: numbers
		that are not truly random, because they are generated by a computer algorithm, but
		which are difficult to predict and uniformly distributed over some interval of
		numbers. What is the Julia function for producing random values?
		""",
		"",
		x -> x[1:4] == "rand"
	),
	Activity(
		"""
		Now we'll look at how uniformly distributed random numbers can generate patterns.
		Load the packages ElectronDisplay and CairoMakie, then create a list (Vector) of
		10 000 integers between 1 and 100. Now use the Makie function hist(data,bins=100) to
		create a histogram showing the distribution of theses numbers. Notice that the heights
		of the bars vary quite widely, but the whole histogram approximates a particular kind of
		statistical distribution. What word describes this kind of distribution?
		""",
		"Notice that although the heights vary, all bars are approximately the same height",
		x -> occursin("uniform",lowercase(x))
	),
	Activity(
		"""
		The heights of these bars vary around a mean value - what is this value?
		""",
		"",
		x -> x==100
	),
	Activity(
		"""
		Now create a list of 10 000 numbers between 1 and 100 by taking the mean of 100 integers
		between 1 and 100. The following line of code does this:

			data = [sum(rand(1:100,100))./100 for _ in 1:10000]

		Now draw a histogramm of this data: hist(data,bins=100). What kind of distribution is this?
		""",
		"",
		x -> occursin("normal",lowercase(x)) || occursin("gauss",lowercase(x))
	),
	Activity(
		"""
		So average values of a uniform distribution form a gaussian distribution. This is the
		content of the Central Limit Theorem of statistics, and is a technique that your own body
		uses to recognise accurately the signal of an antigen from amongst random noise.

		In a computer, we can use statistics to reduce noise. We shall now create a noise-free
		experiment by defining a simple linear function and a range of time and signal values:
		
			t = sort(rand(100));
			signal = 1.3t .- 0.3;

		Now plot this data and tell me the vertical intercept of the resulting plot:

			lines(t,signal)
		""",
		"Don't forget to load ElectronDisplay and CairoMakie first!",
		x -> x == -0.3
	),
	Activity(
		"""
		Now imagine we measured the data x and y, but didn't know the nature of the function f.
		If we know that f is a polynomial function, we can simply test the data against lots of
		possible polynomials to see which of them best fits our data. The CurveFit package does
		this job using least-squares fitting:

			using CurveFit
			fit = curve_fit(Polynomial, t, signal);

		What is the type of the curve_fit() result "fit"?
		""",
		"You may first need to open Package Manager (']') to add the package CurveFit",
		x -> x == Main.Polynomial
	),
	Activity(
		"""
		Now we can find the coefficients of our original signal by listing the components
		of fit: fit[:]. Which component of fit contains the vertical intercept of the signal?
		""",
		"Careful: Notice that Polynomials use more indices than you are used to!",
		x -> x==0
	),
	Activity(
		"""
		Right, now let's add some gaussian noise to our signal to create a noisy experiment:

			noisy = signal + 0.01randn(100);
			lines(t,noisy)

		What is your estimate of the gradient of the graph of noisy?
		""",
		"",
		x -> (1.2 < x < 1.4)
	),
	Activity(
		"""
		Now use curve_fit to tell me the approximate gradient of the noisy experiment:
		""",
		"Look back at how you performed curve-fitting earlier for the clean signal",
		x -> (1.27 < x < 1.33)
	),
	Activity(
		"""
		Our last topic in this laboratory is numerical differentiation. Something we often want
		to do is find the minimum value of a curve given just a few data points. For example,
		suppose we know that the following data has a minimum value somewhere:

			data = besselj.(6,9:.5:15);

		As a first approximation, use the findmin function to find the position of the minimum
		value. Which component of the data is the smallest?
		""",
		"",
		x -> x==6
	),
	Activity(
		"""
		But of course, our data has gaps in it between the datapoints. Is there maybe an even
		smaller value in one of these gaps? One way we can track down the minimum is to find the
		gradient of the data, then find out where this gradient is equal to zero. If t is the
		index of the datapoints in the array, we can define the gradient of the data by
		∂data≡Δ(data)/Δt. To calculate this, we subtract from each datapoint the value of the
		PREVIOUS datapoint, and then divide by the gap Δt=1 between consecutive indices of the
		data. We can do this using the circshift function:

			Δt = 1
			Δdata = data - circshift(data,1)
			∂data = Δdata / Δt

		The first components in your ∂data vector should be negative. How many negative components
		can you count, before they become positive?
		""",
		"To make the symbol Δ, type \\Delta then press <Tab>; for ∂, type \\partial, then <Tab>",
		x -> true
	),
	Activity(
		"""
		The circhshift function is an important tool, so let's take a moment to see how it works.
		circshift rotates the elements of an array circularly in the direction given by its second
		argument, for example:

			a = 1:5
			circshift(a,1)
			circhshift(a,-1)

		We can use circshift with arrays of any size. Enter the following lines and in each case,
		follow the shift in position of the number 1 from its position in the original array b.
		Notice how this shift corresponds to the second argument of circhshift:

			b = reshape(1:16,(4,4))
			circshift(b,(-1,0))
			circshift(b,(0,-2))
			circshift(b,(-1,-2))
		""",
		"",
		x -> true
	),
	Activity(
		"""
		OK, so our gradient ∂data looks good, but it has a problem. To calculate it, we subtracted
		the PREVIOUS datapoint from the current datapoint:

			∂data = (data(t)-data(t-Δt)) / Δt

		In the following graphical figure, we use this method to calculate the derivative at x=1:
		
			t = 0.1:.1:2.1
			lines( t, t.^2)
			lines!( t, t, color=:red)
			current_figure()

		At approximately which point between 0 and 1 is the curve parallel to the straight line?
		""",
		"It's about halfway, isn't it?",
		x -> x == 1.5
	),
	Activity(
		"""
		We just saw that subtracting the PREVIOUS datapoint from the current datapoint gives us the
		value of the gradient just BEFORE the current datapoint x. But instead, we could just as
		easily subtract the current datapoint from the NEXT datapoint:

			∂data = (data(t+Δt)-data(t)) / Δt

		The Mean Value Theorem of mathematics tells us that if we use two datapoints at
		t and t+Δt to calculate the gradient, this only gives us an accurate value for the
		gradient at ONE point somewhere between the two datapoints. To see this, enter the
		following lines of code:

			lines!( t, 3t.-2, color=:green)
			current_figure()

		Now, our new calculation gives us the value of the gradient just _______ the current
		datapoint t:
		""",
		"",
		x -> occursin("after", lowercase(x))
	),
	Activity(
		"""
		Look at the figure of your two gradient calculations. In both cases, the gradient is
		parallel to the curve at a point approximately halfway between the two datapoints. This
		idea is the basis of a much more accurate way of approximating the gradient! To calculate
		the approximate gradient AT the datapoint t, we look at the difference between the NEXT
		and the PREVIOUS datapoints:

			∂data(t) ≈ (data(t+Δt)-data(t-Δt)) / 2Δt

		You can see this idea in action by entering the following yellow line in your figure:

			lines!(x,2x,color=:yellow)
			current_figure()

		At approximately which value of t is the curve parallel to the straight line
		between 0 and 2?
		""",
		"It's about halfway, isn't it?",
		x -> x==1
	),
	Activity(
		"""
		Now let's apply our new trick to the data we created earlier. We'll use circshift to
		calculate at each datapoint the value (NEXT-PREVIOUS) / (2Δx):

			Δdata = circshift(data,-1) - circshift(data,1)
			∂data = Δdata / 2Δt

		Look at the component values of the gradient ∂data and make a guess: at what value of t do
		you estimate the minimum value of data to lie?
		""",
		"At approximately which value of t does the gradient ∂data change its sign?",
		x -> (6.4<x<6.6)
	),
	Activity(
		"""
		This idea of calculating gradients more accurately at the centre of an interval is the basis
		of the Runge-Kutta (R-K) method for solving differential equations (DEs) numerically. We
		will now implement the RK-method, and to test our implementation, we will need a simple test
		case. Consider the exponential equation with growth constant and initial value equal to 1:

			dx/dt = x;   x(0) = 1

		What would you expect to be the solution x(1) of this equation when t = 1?
		""",
		"",
		x -> abs(x-exp(1)) < 0.01
	),
	Activity(
		"""
		We can use our knowledge of the value of x(1) to test the goodness of the Euler and RK
		methods for solving a DE. First we set up Euler's method for solving the exponential
		equation over the interval from t=0 to t=1. To stress-test Euler's method, we shall only
		divide this interval into 4 time-steps:

			x = zeros(5); x[1] = 1; Δt = 0.25;
			m(x) = x
			for k in 2:5
				x[k] = x[k-1] + m(x[k-1]) * Δt
			end
			lines(0:4,x)

		You should see the standard increasing slope of the exponential function. What is the final
		value of x (after the four time-steps)?
		""",
		"You can see it by entering x[end]",
		x -> (2.4<x<2.5)
	),
	Activity(
		"""
		We would expect the last value of x to be something like e = 2.718. However, as you may
		already know, Euler's method suffers from the problem that its successive linear
		approximation steps tend to drift outwards away from the solution curve, making our graph
		inaccurate. Would this explanation lead you to expect your value of x[5] to be too low or
		to be too high?
		""",
		"The exponential function curves upwards, and Euler veers to the outside of this curve",
		x -> occursin("low",lowercase(x))
	),
	Activity(
		"""
		The RK-method reduces this problem greatly by first using Euler to approximate the gradient
		of the solution HALFWAY through the timestep, and then using this halfway value to
		approximate the gradient over a COMPLETE timestep. This usually gives excellent results:

			x = zeros(5); x[1] = 1; Δt = 0.25;
			m(x) = x
			for k in 2:5
				xHalfway = x[k-1] + m(x[k-1]) * Δt/2
				x[k] = x[k-1] + m(xHalfway) * Δt
			end
			lines(0:4,x)

		Again, you should see the increasing slope of the exponential function. What is your final
		value of x now (after the four time-steps)?
		""",
		"You can see it by entering x[end]",
		x -> true
	),
	Activity(
		"""
		As you can see, the RK result is much closer to the actual result of exp(1).

		Congratulations! You have implemented your very own DE solver!

		Before moving on to a well-earned rest, and to the next laboratory, take a couple of
		moments to change the number of time-steps in your solver and see how many you need in
		order to achieve the more accurate value of 2.718.
		""",
		"",
		x -> true
	),
]