diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..cb5f989
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+__pycache__/
+*.py[cod]
+
+plots/*
+log/*
+
diff --git a/CompressiveSim.py b/CompressiveSim.py
new file mode 100755
index 0000000..75e351d
--- /dev/null
+++ b/CompressiveSim.py
@@ -0,0 +1,82 @@
+#!/usr/bin/python3
+"""
+Ryan A. Colyer
+
+Compressive Fluorescence Lifetime Imaging Microscopy simulation testbed.
+
+This simulation is for working out analysis ideas under a controlled system,
+prior to implementation with the experimental data.
+
+This approach is based on a previous MATLAB simulation done as undergraduate
+research by Sarah Grant.
+
+"""
+
+import library.RunTime as rt
+import library.MakeImages as make_img
+import library.Simulator as cssim
+import library.CompressivePlots as cp
+import numpy as np
+
+run_timer = rt.RunTime()
+
+
+randseed = None
+freq = 2*((32/17)*100e6)/8
+frames_per_s = 60
+sim_cps = 20000
+frame_count = 200
+
+print('randseed', randseed)
+print('freq', freq)
+print('frames_per_s', frames_per_s)
+
+np.random.seed(None)
+img = make_img.MakeCellImage(freq)
+total_brightness = np.sum(img.Intensity())
+
+print('total_brightness', total_brightness)
+print('Input', img.Intensity())
+
+
+def Make_g_array(sim):
+  g_arr = []
+  labels = []
+  index = 'a'
+  for frame_count in range(800,3200+1,800):
+    g_row = []
+    l_row = []
+    for sim_cps in range(80000,320000+1,80000):
+      sim.SimAndOutput(sim_cps, frame_count)
+      g_row.append(sim.g_res)
+      l_row.append('('+index+') '+str(sim_cps)+' cps, '+str(frame_count)+' frames')
+      index = chr(ord(index)+1)
+    g_arr.append(g_row)
+    labels.append(l_row)
+
+  cp.PlotArray(g_arr, labels, 'plots/g_array.png', display=False, zscale=[0,1])
+
+def Make_g_noise_data(sim):
+  frame_count = 2400
+  cps = []
+  sdarr = []
+  print('# cps, stdev')
+  for sim_cps in range(20000,320000+1,20000):
+    gvals = []
+    for i in range(10):
+      sim.SimulateImage(sim_cps, frame_count)
+      gvals.append(sim.g_res[8][8])
+    sd = np.std(gvals, ddof=1)
+    sdarr.append(sd)
+    cps.append(sim_cps)
+    print(str(sim_cps)+', '+str(sd))
+
+
+sim = cssim.CSSimulator(img, total_brightness, freq, frames_per_s, run_timer)
+
+#Make_g_array(sim)
+
+#Make_g_noise_data(sim)
+
+sim.SimAndOutput(sim_cps=320000, frame_count=3200)
+
diff --git a/library/CompressivePlots.py b/library/CompressivePlots.py
new file mode 100644
index 0000000..07ea1b6
--- /dev/null
+++ b/library/CompressivePlots.py
@@ -0,0 +1,169 @@
+"""
+Ryan A. Colyer
+
+Generate comparison plots of Compressive Sensing results
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def Plot2(a, Lab_a, b, Lab_b, filename=None, display=True):
+  plt.figure(figsize=(8, 4))
+  plt.set_cmap('nipy_spectral')
+
+  plt.subplot(1,2,1)
+  plt.imshow(np.flipud(a))
+  plt.title(Lab_a)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(1,2,2)
+  plt.imshow(np.flipud(b))
+  plt.title(Lab_b)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplots_adjust(hspace=0.07, wspace=0.07, top=0.91, bottom=0.02, left=0.005, right=0.99)
+
+  if filename:
+    plt.savefig(filename, dpi=300)
+
+  if display:
+    plt.show()
+
+  plt.close()
+
+
+def Plot3(a, Lab_a, b, Lab_b, c, Lab_c, filename=None, display=True):
+  plt.figure(figsize=(11, 4))
+  plt.set_cmap('nipy_spectral')
+
+  plt.subplot(1,3,1)
+  plt.imshow(np.flipud(a))
+  plt.title(Lab_a)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(1,3,2)
+  plt.imshow(np.flipud(b))
+  plt.title(Lab_b)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(1,3,3)
+  plt.imshow(np.flipud(c))
+  plt.title(Lab_c)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplots_adjust(hspace=0.07, wspace=0.07, top=0.95, bottom=0.02, left=0.005, right=0.99)
+
+  if filename:
+    plt.savefig(filename, dpi=300)
+
+  if display:
+    plt.show()
+
+  plt.close()
+
+
+def Plot4(a, Lab_a, b, Lab_b, c, Lab_c, d, Lab_d, filename=None, display=True, zscale=[None, None, None, None]):
+  plt.figure(figsize=(8, 7))
+  plt.set_cmap('nipy_spectral')
+
+  plt.subplot(2,2,1)
+  if zscale[0]:
+    plt.pcolor(a, vmin=zscale[0][0], vmax=zscale[0][1])
+  else:
+    plt.pcolor(a)
+  plt.title(Lab_a)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(2,2,2)
+  if zscale[1]:
+    plt.pcolor(b, vmin=zscale[1][0], vmax=zscale[1][1])
+  else:
+    plt.pcolor(b)
+  plt.title(Lab_b)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(2,2,3)
+  if zscale[2]:
+    plt.pcolor(c, vmin=zscale[2][0], vmax=zscale[2][1])
+  else:
+    plt.pcolor(c)
+  plt.title(Lab_c)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplot(2,2,4)
+  if zscale[3]:
+    plt.pcolor(d, vmin=zscale[3][0], vmax=zscale[3][1])
+  else:
+    plt.pcolor(d)
+  plt.title(Lab_d)
+  plt.colorbar()
+  plt.axis('off')
+
+  plt.subplots_adjust(hspace=0.07, wspace=0.07, top=0.95, bottom=0.02, left=0.005, right=0.99)
+
+  if filename:
+    plt.savefig(filename, dpi=300)
+
+  if display:
+    plt.show()
+
+  plt.close()
+
+
+def PlotArray(data_arr, label_arr, filename=None, display=True, zscale=None):
+  w = len(data_arr[0])
+  h = len(data_arr)
+  plt.figure(figsize=(3*w+1, 3.5*h))
+  plt.set_cmap('nipy_spectral')
+
+  index = 1
+  for (x,y) in ((x,y) for y in range(h) for x in range(w)):
+    plt.subplot(w, h, index)
+    plt.title(label_arr[y][x])
+    if zscale:
+      plt.pcolor(data_arr[y][x], vmin=zscale[0], vmax=zscale[1])
+    else:
+      plt.pcolor(data_arr[y][x])
+    
+#    if x == w-1:
+#      plt.colorbar()
+
+    plt.axis('off')
+
+    index += 1
+  
+  plt.subplots_adjust(hspace=0.07, wspace=0.07, top=0.95, bottom=0.02, left=0.005, right=0.99)
+
+  if filename:
+    plt.savefig(filename, dpi=300)
+
+  if display:
+    plt.show()
+
+  plt.close()
+
+
+def IntensityComparison(original, result, filename=None, display=True):
+  Plot2(original, '(a) Simulation Input', result, '(b) Compressive Sensing', filename, display)
+
+def IGSComparison(original, intensity, bigG, bigS, filename=None, display=True):
+  Plot4(original, '(a) Simulation Input', intensity, '(b) CS Intensity', bigG, '(c) CS BigG', bigS, '(d) CS BigS', filename, display)
+
+def IgsComparison(original, intensity, g, s, filename=None, display=True):
+  Plot4(original, '(a) Simulation Input', intensity, '(b) CS Intensity', g, '(c) CS g', s, '(d) CS s', filename, display, zscale=[None,None,[0,1],[0,1]])
+
+def ITauPComparison(original, intensity, taup, filename=None, display=True):
+  Plot3(original, '(a) Simulation Input', intensity, '(b) CS Intensity', taup, '(c) CS TauP', filename, display)
+
+def ITauMComparison(original, intensity, taum, filename=None, display=True):
+  Plot3(original, '(a) Simulation Input', intensity, '(b) CS Intensity', taum, '(c) CS TauM', filename, display)
+
diff --git a/library/CompressiveSensing.py b/library/CompressiveSensing.py
new file mode 100644
index 0000000..dea3e83
--- /dev/null
+++ b/library/CompressiveSensing.py
@@ -0,0 +1,65 @@
+"""
+Ryan A. Colyer
+
+Compressive Sensing Analysis tools.
+
+"""
+
+import numpy as np
+import random
+from library.Phasor import Phasor
+from sklearn.linear_model import Lasso
+
+# Stores the data for one compressive FLIM frame, consisting of a mask,
+# the intensity, and the phasor data.
+class Frame:
+  def __init__(self, mask):
+    self.mask = mask
+    self.data = Phasor()
+
+  def Add(self, photon):
+    self.data += Phasor(photon)
+
+
+# A randomly generated mask.
+class Mask:
+  def __init__(self, width, height):
+    self.width = width
+    self.height = height
+    self.data = [[np.random.randint(0,2) for x in range(width)] for y in range(height)]
+
+  def IsOn(self, x, y):
+    return 1 == self.data[y][x]
+
+  def IsOff(self, x, y):
+    return not self.IsOn(x, y)
+
+
+# Reconstruct the base image with L1 penalization using Lasso.
+class Reconstruction:
+  def __init__(self, alpha=0.001):
+    self.model = Lasso(alpha)
+  
+  def ModelFit(self, frames, data):
+    height = frames[0].mask.height
+    width = frames[0].mask.width
+    masks = self.GetMasks(frames)
+    self.model.fit(masks, data)
+    self.result = self.model.coef_.reshape(height, width)
+    return self.result
+
+  def GetMasks(self, frames):
+    return np.array([np.array(f.mask.data).ravel() for f in frames])
+
+  def FitIntensity(self, frames):
+    counts = np.array([f.data.count for f in frames])
+    return self.ModelFit(frames, counts)
+
+  def FitG(self, frames):
+    bigGs = np.array([f.data.big_G for f in frames])
+    return self.ModelFit(frames, bigGs)
+
+  def FitS(self, frames):
+    bigSs = np.array([f.data.big_S for f in frames])
+    return self.ModelFit(frames, bigSs)
+
diff --git a/library/ImageGen.py b/library/ImageGen.py
new file mode 100644
index 0000000..da56d67
--- /dev/null
+++ b/library/ImageGen.py
@@ -0,0 +1,96 @@
+"""
+Ryan A. Colyer
+
+Image simulation tools for compressive sensing with FLIM.
+
+"""
+
+from abc import ABC, abstractmethod
+from math import log
+import random
+from library.Phasor import Photon
+import library.CompressiveSensing as cs
+import copy
+import numpy as np
+
+
+# intensity unit 1 for always on.
+# lifetime unit 1 for laser period.  Special value 0 is uncorrelated light.
+class Pixel:
+  def __init__(self, intensity, lifetime=0):
+    self.intensity = intensity
+    self.lifetime = lifetime
+
+  # Generate random fluorescence lifetime nanotime with periodicity 1.
+  def GetNanotime(self, rng):
+    if self.lifetime == 0:
+      return rng.random()
+    else:
+      return (-self.lifetime*log(rng.random())) % 1
+
+
+
+# Base class for image simulation.
+class ImageSimulator(ABC):
+  def __init__(self, width, height):
+    self.width = width
+    self.height = height
+
+  def GetDim(self):
+    return (self.width, self.height)
+
+  @abstractmethod
+  def GetPixel(self, x, y):
+    return NotImplemented
+
+
+# Provides a static intensity image with uncorrelated (long lifetime) counts.
+class StaticImage(ImageSimulator):
+  def __init__(self, width, height):
+    super().__init__(width, height)
+    self.image_data = [[Pixel(0) for x in range(width)] for y in range(height)]
+
+  def GetPixel(self, x, y):
+    return self.image_data[y][x]
+
+  def SetPixels(self, pixel, xy_list):
+    for xy in xy_list:
+      self.image_data[xy[1]][xy[0]] = copy.deepcopy(pixel)
+
+  def Intensity(self):
+    return [[self.image_data[y][x].intensity for x in range(self.width)] for y in range(self.height)]
+
+
+# Base class for a photon detection sequence.
+class PhotonSource(ABC):
+  @abstractmethod
+  def GetFrame(self):
+    return NotImplemented
+
+
+# Provides a photon detection sequence simulator based on a given image source.
+class PhotonSimSource(PhotonSource):
+  def __init__(self, image_source, sim_intensity):
+    super().__init__()
+    self.image_source = image_source
+    self.sim_intensity = sim_intensity
+
+  def AddPixelPhotons(self, x, y, iframe):
+    pixel = self.image_source.GetPixel(x, y)
+    pix_avgcnt = self.sim_intensity * pixel.intensity;
+    ph_cnt = np.random.poisson(pix_avgcnt)
+    for ph in range(ph_cnt):
+      nanotime = pixel.GetNanotime(np.random)
+      photon = Photon(x, y, nanotime)
+      iframe.Add(photon)
+
+  def GetFrame(self):
+    height = self.image_source.height
+    width = self.image_source.width
+    mask = cs.Mask(width, height)
+    iframe = cs.Frame(mask)
+    for (x,y) in ((x,y) for y in range(height) for x in range(width)):
+      if mask.IsOn(x, y):
+        self.AddPixelPhotons(x, y, iframe)
+    return iframe
+
diff --git a/library/MakeImages.py b/library/MakeImages.py
new file mode 100644
index 0000000..a43c7c3
--- /dev/null
+++ b/library/MakeImages.py
@@ -0,0 +1,49 @@
+"""
+Ryan A. Colyer
+
+Predefined images for simulation purposes.
+"""
+
+import library.ImageGen as ig
+
+
+def MakeSmileImage(freq):
+  dim = [8, 8]
+  img = ig.StaticImage(*dim)
+
+  spec1 = ig.Pixel(1, 4.08e-9*freq) # Rhodamine 6G
+  background = ig.Pixel(0.3, 0*freq) # Uncorrelated background light
+
+  # Draw out a smiley face test image for the theory paper.
+  img.SetPixels(background, [[x,y] for x in range(dim[0]) for y in range(dim[1])])
+  img.SetPixels(spec1, [[2,2],[5,2],[0,5],*[[x,6] for x in range(1,7)],[7,5]])
+
+  return img
+
+
+def MakeCellImage(freq):
+  dim = [16, 16]
+  img = ig.StaticImage(*dim)
+
+  GFP = ig.Pixel(0.8, 3.20e-9*freq) # Membrane bound
+  DAPI = ig.Pixel(1, 2.20e-9*freq) # Nuclear stain
+  autofluor = ig.Pixel(0.4, 100e-9*freq) # Long-lifetime autofluorescence
+  background = ig.Pixel(0.2, 0*freq) # Uncorrelated background light
+
+  img.SetPixels(background, [[x,y] for x in range(dim[0]) for y in range(dim[1])])
+  img.SetPixels(GFP, [[8,2],[9,2],[7,3],[10,3],[11,3],[4,4],[5,4],[6,4],[12,4],[13,4],[3,5],[14,5],[2,6],[14,6],[2,7],[14,7],[1,8],[14,8],[1,9],[13,9],[2,10],[12,10],[3,11],[12,11],[4,12],[11,12],[5,13],[6,13],[9,13],[10,13],[7,14],[8,14]])
+  img.SetPixels(autofluor, [[x,3] for x in range(8,10)])
+  img.SetPixels(autofluor, [[x,4] for x in range(7,12)])
+  img.SetPixels(autofluor, [[x,5] for x in range(4,14)])
+  img.SetPixels(autofluor, [[x,6] for x in range(3,14)])
+  img.SetPixels(autofluor, [[x,7] for x in range(3,14)])
+  img.SetPixels(autofluor, [[x,8] for x in range(2,14)])
+  img.SetPixels(autofluor, [[x,9] for x in range(2,13)])
+  img.SetPixels(autofluor, [[x,10] for x in range(3,12)])
+  img.SetPixels(autofluor, [[x,11] for x in range(4,12)])
+  img.SetPixels(autofluor, [[x,12] for x in range(5,11)])
+  img.SetPixels(autofluor, [[x,13] for x in range(7,9)])
+  img.SetPixels(DAPI, [[7,7],[8,7],[6,8],[7,8],[8,8],[9,8],[6,9],[7,9],[8,9],[9,9],[7,10],[8,10]])
+
+  return img
+
diff --git a/library/Phasor.py b/library/Phasor.py
new file mode 100644
index 0000000..6f6ebba
--- /dev/null
+++ b/library/Phasor.py
@@ -0,0 +1,89 @@
+"""
+Ryan A. Colyer
+
+Store phasor information, corresponding to the linearly accumulated real and
+imaginary components of the first harmonic of the Fourier transform of the
+fluorescence lifetime response for a given photon, pixel, or area.
+
+"""
+
+from math import cos, sin, tan, atan, atan2, pi, sqrt
+from numpy import sign
+
+# Stores photon data with first harmonic phasor information.
+class Photon:
+  def __init__(self, x, y, nanotime):
+    self.x = x
+    self.y = y
+    self.nanotime = nanotime
+    self.g = cos(nanotime*2*pi)
+    self.s = sin(nanotime*2*pi)
+
+# Properly accumulates phasor information for a given point.
+class Phasor:
+  def __init__(self, photon=None):
+    if photon is None:
+      self.big_G = 0
+      self.big_S = 0
+      self.count = 0
+    else:
+      self.big_G = photon.g
+      self.big_S = photon.s
+      self.count = 1
+
+  def __add__(self, phasor):
+    self.big_G += phasor.big_G
+    self.big_S += phasor.big_S
+    self.count += phasor.count
+    return self
+
+  # Return the canonical g coordinate for a phasor.
+  def g(self):
+    return self.big_G / self.count
+
+  # Return the canonical s coordinate for a phasor.
+  def s(self):
+    return self.big_S / self.count
+
+# Return the intensity normalized 2D g and s data.
+def NormalizePhasorData(bigG, bigS, intensity):
+  g = [[G/I for G,I in zip(Grow, Irow)] for Grow,Irow in zip(bigG, intensity)]
+  s = [[S/I for S,I in zip(Srow, Irow)] for Srow,Irow in zip(bigS, intensity)]
+  return (g,s)
+
+# Return a phi image.
+def GetPhiData(g_data, s_data):
+  return [[atan2(s, g) for g,s in zip(grow, srow)] for grow,srow in zip(g_data, s_data)]
+
+# Return a modulation image.
+def GetModData(g_data, s_data):
+  return [[sqrt(g*g+s*s) for g,s in zip(grow, srow)] for grow,srow in zip(g_data, s_data)]
+
+
+phicap_maxv = atan(100)
+def PhiCap(phi):
+  maxv = atan(100)
+  if phi > maxv:
+    return maxv
+  else:
+    return phi
+
+modcap_minv = sqrt(1/(100**2+1))
+def ModCap(m):
+  if m < modcap_minv:
+    return modcap_minv
+  else:
+    return m
+
+# Return a TauP image.
+def GetTauPData(g_data, s_data, freq):
+  phi_data = GetPhiData(g_data, s_data)
+  omega = 2 * pi * freq
+  return [[tan(PhiCap(phi))/omega for phi in phi_row] for phi_row in phi_data]
+
+# Return a TauM image.
+def GetTauMData(g_data, s_data, freq):
+  m_data = GetModData(g_data, s_data)
+  omega = 2 * pi * freq
+  return [[sign(m)*sqrt(abs(1/(ModCap(m)**2)-1))/omega for m in m_row] for m_row in m_data]
+
diff --git a/library/RunTime.py b/library/RunTime.py
new file mode 100644
index 0000000..389e02e
--- /dev/null
+++ b/library/RunTime.py
@@ -0,0 +1,31 @@
+"""
+Ryan A. Colyer
+
+Report process runtime.
+
+"""
+
+from datetime import datetime
+
+
+class RunTime:
+  def __init__(self):
+    self.Start()
+
+  def Start(self):
+    self.start = datetime.now()
+
+  def End(self):
+    self.end = datetime.now()
+
+  def Duration(self):
+    return self.end - self.start
+
+  def DurationStr(self):
+    tlst = str(self.Duration()).split(':')
+    return tlst[0] + 'h ' + tlst[1] + 'm ' + '{0:.2f}s'.format(float(tlst[2]))
+
+  def SinceStart(self):
+    self.End()
+    return self.DurationStr()
+
diff --git a/library/Simulator.py b/library/Simulator.py
new file mode 100644
index 0000000..38c3332
--- /dev/null
+++ b/library/Simulator.py
@@ -0,0 +1,64 @@
+"""
+Ryan A. Colyer
+
+The simulator class for compressive sensing with FLIM.
+"""
+
+import library.ImageGen as ig
+import library.CompressiveSensing as cs
+import library.Phasor as ph
+import library.CompressivePlots as plot
+
+
+class CSSimulator:
+  def __init__(self, img, total_brightness, freq, frames_per_s, run_timer):
+    self.img = img
+    self.total_brightness = total_brightness
+    self.freq = freq
+    self.frames_per_s = frames_per_s
+    self.input_intensity = img.Intensity()
+    self.run_timer = run_timer
+
+  # Acquire the full list of compressive sensing frames.
+  def GetFrames(self, ph_source, frame_count):
+    return [ph_source.GetFrame() for f in range(frame_count)]
+
+  # Run one compressive sensing simulation.
+  def SimulateImage(self, sim_cps, frame_count):
+    self.sim_intensity = (sim_cps/self.frames_per_s)/(self.total_brightness/2)
+
+    ph_source = ig.PhotonSimSource(self.img, self.sim_intensity)
+    self.frames = self.GetFrames(ph_source, frame_count=frame_count)
+
+    l1min = cs.Reconstruction()
+
+    self.intensity_result = l1min.FitIntensity(self.frames)
+    self.G_result = l1min.FitG(self.frames)
+    self.S_result = l1min.FitS(self.frames)
+
+    (self.g_res, self.s_res) = ph.NormalizePhasorData(self.G_result, self.S_result, self.intensity_result)
+    self.taup_res = ph.GetTauPData(self.g_res, self.s_res, self.freq)
+    self.taum_res = ph.GetTauMData(self.g_res, self.s_res, self.freq)
+
+  def SimAndOutput(self, sim_cps, frame_count):
+    self.SimulateImage(sim_cps, frame_count)
+
+    print('sim_cps', sim_cps)
+    print('frame_count', frame_count)
+    print('Intensity', self.intensity_result)
+    print('sim_intensity,', self.sim_intensity)
+    print('frames[0] count', self.frames[0].data.count)
+    print("g", self.g_res)
+    print("s", self.s_res)
+    print('taup', self.taup_res)
+    print('taum', self.taum_res)
+    print('Simulation run time: ', self.run_timer.SinceStart())
+
+    pltdir = 'plots/'
+    suffix = '_' + str(sim_cps) + '_' + str(frame_count)
+    plot.IntensityComparison(self.input_intensity, self.intensity_result, pltdir+'intensity_comparison'+suffix+'.png', display=False)
+    plot.ITauPComparison(self.input_intensity, self.intensity_result, self.taup_res, pltdir+'taup_comparison'+suffix+'.png', display=False)
+    plot.ITauMComparison(self.input_intensity, self.intensity_result, self.taum_res, pltdir+'taum_comparison'+suffix+'.png', display=False)
+    plot.IgsComparison(self.input_intensity, self.intensity_result, self.g_res, self.s_res, pltdir+'gs_comparison'+suffix+'.png', display=False)
+
+
diff --git a/run b/run
new file mode 100755
index 0000000..a1befbe
--- /dev/null
+++ b/run
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+DATE="$(date +%Y-%m-%d_%H-%M-%S)"
+LOGFILE="log/simulation_$DATE.log"
+
+python3 CompressiveSim.py >"$LOGFILE"
+