cochlea-positional-information_fork-test/notebooks/sandbox-pluto-nb backup 4.jl

### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 485d8da0-be94-11eb-183b-e5e105e25b55
begin
	using DrWatson
	using PlutoUI
	using JuMP
	using Optim
	using Random
	import Ipopt
	import Juniper
	import Alpine
	# using NLopt
	using Downloads
	using ZipFile
	using DataFrames
	using MAT
	using MATLAB
	using Plots
	using Statistics
	using NaNStatistics
	using Images
	using DelimitedFiles
	using CSV
	using JLD
	using Interpolations
	using Convex, SCS # other solvers: GLPK, Mosek, Gurobi, ECOS
	using Distributions
	using StatsBase
	using LinearAlgebra
	plotly()
	PlutoUI.TableOfContents()
end

# ╔═╡ 8d62b88b-f77f-4d8d-8c6a-ada1d1626b8c
let
	using Colors, ColorSchemes
using CuArrays, GPUArrays

# get the number of steps
function get_steps(c::Complex, max_steps)
    z = Complex(0.0, 0.0) # 0 + 0im
    for i=1:max_steps
        z = z*z+c
        if abs2(z) >= 4
            return i
        end
    end
    return max_steps+1
end

function get_color(colorscheme, step, max_steps)
    if step == max_steps+1
        return [0.0, 0.0, 0.0]
    end
    color = get(colorscheme, step, (1, max_steps))
    return [color.r, color.g, color.b]
end

function get_cmap(colorscheme, max_steps)
    colors = zeros(Float64, (3, max_steps+1))
    for i=1:max_steps 
        colors[:,i] = get_color(colorscheme, i, max_steps)
    end
    colors[:,max_steps+1] = [0.0, 0.0, 0.0]
    return colors
end

function mandelbrot_plot()
    width = 1000
    height = 600

    max_steps = 3500
    steps = zeros(Int, (height, width))
   
    # range for real values
    cr_min = -0.7491597623
    cr_max = -0.7491597623+0.0000000004

    # range for imaginary values
    ci_min = 0.1005089256
    
    range = cr_max - cr_min
    dot_size = range/width
    ci_max = ci_min + height*dot_size

    # println("cr: $cr_min - $cr_max")
    # println("ci: $ci_min - $ci_max")

    image = zeros(Float64, (3, height, width))
    complexes = zeros(ComplexF64, (height, width))
    steps = zeros(Int, (height, width))
    cu_steps = CuArray(zeros(Int, (height, width)))

    colorscheme = ColorSchemes.inferno
    colorscheme_sized = get_cmap(colorscheme, max_steps)
    x, y = 1,1
    for ci=ci_min:dot_size:ci_max-dot_size
        x = 1
        for cr=cr_min:dot_size:cr_max-dot_size
            complexes[y,x] = Complex(cr, ci)
            x += 1
        end
        y += 1
    end
    cu_complexes = CuArray(complexes)
    cu_steps .= get_steps.(cu_complexes, max_steps)
    GPUArrays.synchronize(cu_steps)
    steps = Array(cu_steps)
    
    image = colorscheme_sized[:, steps]
    
    save("images/test.bmp", colorview(RGB, image))
end

@time mandelbrot_plot();
end

# ╔═╡ a4dd3c3d-26a8-4c55-8dde-3c7f6301f01e
# begin
# 	using Pkg
# 	Pkg.activate("..")
# 	Pkg.instantiate()
# end

# ╔═╡ ea9907cf-6f0e-495e-a172-bc6fd0cf2863
# @quickactivate

# ╔═╡ 4f275ac2-e986-4eec-bcef-fd7178e96d3d
#= html"""
<style>
  main {
    max-width: 1200px;
  }
</style>
""" =#

# ╔═╡ f0638606-63c0-4541-bd76-2887002557f0
md"""
# Interactive Stuff
In this section, we'll try making interactive features.
"""

# ╔═╡ 8b0d2a11-ca25-4d95-bb74-33a5a8f5b5e5
md"""
## Interactive Plotting
"""

# ╔═╡ 69573234-68a8-4e41-bd34-7c5bb013e7f9
md"""
Set the range of x and the amplitude for the signal (A) and noise (N).
"""

# ╔═╡ ed292325-e42e-43af-936f-6a98b1c86e3f
md"""
xₘₐₓ = $(@bind xₘₐₓ Slider(0:0.01:62.8, default=7.0, show_value=true))
|
A = $(@bind A Slider(0.01:0.01:10, default=1.0, show_value=true))
|
N = $(@bind N Slider(0.0:0.01:5.0, default=0.1, show_value=true))


Replicates = $(@bind n_signal_replicates Slider(1:1000, default=10, show_value=true))
"""

# ╔═╡ d4e4a684-1b0f-45d5-a377-35f005333bcc
x=(0:0.01:xₘₐₓ)

# ╔═╡ b40521ca-9982-4c61-954d-920a766380d2
begin
	x_mat = repeat(x,1,n_signal_replicates);
	noisy_signals = A.*sin.(x_mat) + N.*randn(Float64,size(x_mat))
	μ_signal = mean(noisy_signals,dims=2)
	σ_signal = std(noisy_signals,dims=2)
	plot(x,μ_signal,
	color="black", linewidth=2, legend=false, ribbon=σ_signal, ylims = (-1.5*A, 1.5*A))
	title!("Carrier Signal in Noise")

end

# ╔═╡ 097cad13-82b7-4859-adc0-4a0bba8bc632
md"""
# Drosophila Data
"""

# ╔═╡ ab4a23bb-dd9e-4455-ba6a-08e007f7092e
md"""
OK, well that was fun. Now, let's try to get the Drosophila data from the Cell paper. The data was downloaded from "https://ars.els-cdn.com/content/image/1-s2.0-S0092867419300406-mmc4.zip", then extracted into `cochlea-positional-information > data > exp_raw > Data_doi.org10.1016j.cell.2019.01.007`
"""

# ╔═╡ 51e30cad-30b7-44ff-adf3-d06018541e16
md"""
## ✅ Access drosophila data.
"""

# ╔═╡ 6d6bfc47-b00d-4575-a62f-ac99fd37e641
vars = matread(joinpath("..","data","exp_raw","doi.org10.1016j.cell.2019.01.007","Gap","gap_data_raw_dorsal_wt.mat"))

# ╔═╡ b165c538-a5f5-4ea6-8ce5-6d60943b3cfa
n_replicates = size(vars["data"]["index"])[2]

# ╔═╡ 7456dc7a-fed8-4134-ba4c-f05ad99665f4
gene_set = ["Hb","Gt","Kr","Kni"]

# ╔═╡ 19fa8dc3-fc71-4cfa-87c7-961d94fc2ade
md"""
Extract the data from the MATLAB struct into a data type native to Julia.
Let's first see how we can index into the struct array pulled in from the .mat file and plot one profile at a time for each gene.
"""

# ╔═╡ ca5b9a0e-59b4-478d-9347-f500a364a18d
md"""
Now, let's get the data out of the struct array and put it into a dictionary.
"""

# ╔═╡ c8a39986-acb9-41ca-bb0b-ad5d005d4c08
begin
	# Initialize an empty dictionary with gene names as keys and empty matrices as values to populate in the loop.
	gene_data = Dict(name => Matrix{Float64}(undef,1000,102) for name in gene_set)
	
	for gene in gene_set
			gene_data 
		for i_replicate = 1:102
			gene_data[gene][:,i_replicate] = vars["data"][gene][i_replicate]'
		end
	end
end

# ╔═╡ 335bae5d-d0b9-4620-9ec9-71e427bf5dea
begin
	gene_data_array = Array{Float64}(undef,1000,102,4)
	for i_gene in 1:4
		gene_data_array[:,:,i_gene] = gene_data[gene_set[i_gene]]
	end
	n_points, ~, n_targets = size(gene_data_array)
end

# ╔═╡ 0b49770b-a6cb-44bf-a1e1-f9db3e361cb5
size(vars["data"]["Hb"][1])

# ╔═╡ e64c4b85-3550-43d2-bdf6-8979ff5fcb27
size(gene_data["Hb"][:,1])

# ╔═╡ 8cf2b1e5-a403-4e69-aabd-d712b3b746ac
md"""
Move the slider below to select an embryo from the dataset.
"""

# ╔═╡ 1fb671db-4ee7-4756-87ab-7970a701e7bb
@bind i_replicate Slider(1:1:n_replicates, show_value=true)

# ╔═╡ e56f2074-5f94-43ba-9b10-662f65723033
begin
	x_dros = collect(range(0.001,1,length=1000))
	p = plot()
	idx=0
	for gene in gene_set
		
		μ_dros = nanmean(gene_data[gene],dims=2)		
		σ_dros = std(gene_data[gene],dims=2)
		
		μ_dros_v = vec(μ_dros)
		σ_dros_v = vec(σ_dros)
		
		n_pts = length(σ_dros)
		
		nan_idcs = findall(isnan,σ_dros_v)
		
		x_nan = x_dros
	
		x_nan[nan_idcs] .= NaN
		
		bnds = zeros(2)
		for idx_σ = 2:n_pts-1
			if !isnan(σ_dros[idx_σ]) & isnan(σ_dros[idx_σ-1])
				bnds[1] = idx_σ
			end
			if !isnan(σ_dros[idx_σ]) & isnan(σ_dros[idx_σ+1])
				bnds[2] = idx_σ
			end
		end
		
		idx = idx+1
# <<<<<<< HEAD
# 		# plot!(p,x_dros,gene_data[gene],
# 		# 	color=idx, legend=false, ylims=[0,2100],alpha=0.05)
# =======
# 		# All individual profiles
# 		plot!(p,x_dros,gene_data[gene],
# 		color=idx, legend=false, ylims=[0,2100],alpha=0.05)
# >>>>>>> dev-matt
		
		# Highlight one individual profile
		plot!(p,x_dros,gene_data[gene][:,i_replicate],
			color=idx, ylims=[0,2100],linewidth=1)
		
		σ_dros = nanstd(gene_data[gene],dims=2)
		
		plot!(p,x_dros,nanmean(gene_data[gene],dims=2),
			color=idx, ylims=[0,2100],linewidth=4,ribbon=σ_dros)
			# color=idx, ylims=[0,2100],linewidth=1.5)
# >>>>>>> dev-matt
		
		# Mean and standard deviation
		plot!(p,x_nan,μ_dros,
			color=idx, ylims=[0,2100],linewidth=4,ribbon=σ_dros)
	end
	current()
end

# ╔═╡ 7123c8bc-7e48-4d7a-91dc-d35885407ac7
md"""
We have the data in a reasonable format now. So let's process it.
"""

# ╔═╡ 05d085bb-878d-4e93-9985-373e9955892a
md"""
## ✅ Anchor mean from 0 to 1.
"""

# ╔═╡ ef7a4517-be6e-46ab-b85f-f3254e61905a
md"""
## ✅ Align Y with JuMP.jl
"""

# ╔═╡ 7170e92a-ab07-4ee9-aad3-dc30783f016e
# Plot drosophila alignment results.

# ╔═╡ b9c644e8-86c4-4935-8117-bfde3173b0f7
md"""
____________________________________________________________________________________
"""

# ╔═╡ 6cb4b869-dcd0-414b-88c5-67bf69d084a2
md"""
# Cochlea Data
"""

# ╔═╡ 3e9036f5-5d11-4950-bba8-a2acfb84bf5e
md"""
## ✅ Access cochlea images and data.
"""

# ╔═╡ 9fabf548-d740-4989-a0e6-4700f8436cfe
### Set up path string lists.

begin
	e12_img_dir = joinpath("..","data","exp_pro","sw","wt","e12.5","psmad-sox2-topro3","tif-orient","sep-ch","mip","_data_")
	e13_img_dir = joinpath("..","data","exp_pro","sw","wt","e13.5","psmad-sox2-topro3","tif-orient","sep-ch","mip","_data_")
	
	e12_img_paths = filter(endswith(".tif"),readdir(e12_img_dir,join=true))
	e13_img_paths = filter(endswith(".tif"),readdir(e13_img_dir,join=true))
	
	img_paths = Dict([("e12",e12_img_paths),("e13",e13_img_paths)])

	e12_prof_paths = filter(endswith(".csv"),readdir(normpath(joinpath(e12_img_dir,"..","quant-prof")),join=true))
	e13_prof_paths = filter(endswith(".csv"),readdir(normpath(joinpath(e13_img_dir,"..","quant-prof")),join=true))
	
	prof_paths = Dict([("e12",e12_prof_paths),("e13",e13_prof_paths)])
end

# ╔═╡ b695d442-89e2-4f21-aca6-52d26ca305c2
md"""
## ✅ Get the data into a dataframe.
"""

# ╔═╡ 4c86e3ce-cc11-4644-88d9-c52cb41d09e3
md"""
## ✅ Resample raw data uniformly.
"""

# ╔═╡ ef16f350-ab20-4aa1-ae2f-d68f0d6f8f00
md"""
## 🚧 Exclude outer 10%.
"""

# ╔═╡ a80ca2b0-c4bf-4d9d-9ead-9fe5fff2ae58
md"""
## ✅ Regress nuclear noise.
"""

# ╔═╡ 6895787b-70f3-4cd1-9e9f-243c22662f99
md"""
## ✅ Anchor mean from 0 to 1.
"""

# ╔═╡ 8e6682fa-6790-4275-9ad2-3aa9c1c72c1c
md"""
## 🚧 Smooth data.
"""

# ╔═╡ c3c4517b-aca3-437f-a3ca-8653b91b1e8a


# ╔═╡ 1689ecbf-0750-46c1-82f2-d0e598d15448
md"""
## ✅ Align Y with JuMP.jl.

Using MATLAB's `fmincon` is a straightforward implementation of a $\chi$``^2`` alignment.

With Julia, we can use JuMP.
"""

# ╔═╡ f69ce4e9-3f7f-4928-ad8f-b8d98384565f
md"""
### ⚠ Mean should be same before and after y-align.
"""

# ╔═╡ 3bf6395b-acf2-4396-8b1a-f65291df41ed
md"""
## 🚧 Align X and Y with JuMP.jl.
"""

# ╔═╡ 8f6133fd-2508-4e06-97a9-65985b779a1c
md"""
### 🚧 begin construction 🚧
"""

# ╔═╡ 92bceb54-b2a0-4d82-a2d4-f3770c141834
ones(2,2)

# ╔═╡ bc001ccd-321b-4a9c-80b9-9bf36e79335a
# # ✅ Initialize a 3D matrix (n_points x n_replicates x n_targets)
# begin
# 	Y_0 = zeros(length(df[1,:topro3_yalign]),nrow(df),2)
# end

# ╔═╡ 1c0ec19a-4229-43e3-83fd-2c07c03b8e44
# nan_pad_percent = 0.1

# ╔═╡ a1965e64-9449-47b2-a706-ebe36080f2a3
# n_nans = convert(Int,round(nan_pad_percent*n_p))

# ╔═╡ c5cdaccd-00b6-4fa8-b2ad-294feead5c90
# nan_pad = replace!(zeros(n_nans,n_r,n_t),0=>NaN)

# ╔═╡ 1934c23e-478f-4404-ad4c-b37c33f526a2
# begin
# 	shift_array = round.(rand(Uniform(-n_nans,n_nans),1,11))
# 	Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
# 	Y_0_shift = Y_0_padded
# 	for i_r = 1:n_r
# 		Y_0_shift[:,i_r,:] = circshift(Y_0_padded[:,i_r,:],shift_array[i_r])
# 	end
# 	ybar = reshape(nanmean(Y_0_shift,dims=2),n_p+2*n_nans,2)
# end

# ╔═╡ d9edc6c1-07cd-417e-9f81-cb5db71c5185
model_xy = Model(Ipopt.Optimizer)

# ╔═╡ ed85ac59-4332-4c76-8e35-a1e193d5d185
# function xy_objective(α,β,γ,nan_pad_percent,Y_0)
	
# 	n_p, n_r, n_t = size(Y_0)
# 	n_nans = convert(Int,round(nan_pad_percent*n_p))
	
# 	# 1. circshift by γ
	
# 	shift_array = round.(rand(Uniform(-n_nans,n_nans),1,11))
# 	nan_pad = replace!(zeros(n_nans,n_r,n_t),0=>NaN)
# 	Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
# 	Y_0_shift = Y_0
# 	for i_r = 1:n_r
# 		Y_0_shift[:,i_r,:] = circshift(Y_0_padded[:,i_r,:],γ[i_r])
# 	end
	
# 	circshift(Y,(γ,0))
# 	# 2. take nanmean of each target to reduce variance against
# 	# ybar = reshape(nanmean(Y_0_shift,dims=2),1200,2) # for reference

# 	ȳ = reshape(nanmean(Y_0_shift,dims=2),n_p + 2*n_nans, n_t)
# 	# 3. calculate cost
	
# 	χ² = sum((Y_0_shift - (reshape(α,(1,n_replicates)) .+ ȳ*reshape(β,(1,n_replicates)))).^2)


# 	return χ²
# end

# ╔═╡ 7420d1ae-dd36-478e-96c8-f3c65a3746ce


# ╔═╡ 45cf206d-7405-464c-a86d-a4cfa8af3119
md"""
### 🚧 end construction zone 🚧
"""

# ╔═╡ 9f64c871-04a4-4d33-ac5d-1bc49b4d0485
# xy_align script implementation

# ╔═╡ 64b5deca-ab94-4f74-b960-711d2df20894
md"""
## 🚧 Positional error.
"""

# ╔═╡ 16a12f25-b06a-496c-bd66-c57ff3e15783
md"""
## 👀 🚧 Direct estimate of I({g};x)
"""

# ╔═╡ 0bf8238e-4217-45cb-bb78-28d6ba420a8e
begin
	# Using drosophila data
	# size(gene_data_array)
	# 	n_points
	# 	n_replicates
	# 	n_targets
	xx = collect(1/n_points:1/n_points:1)[100:900]
	
end

# ╔═╡ 158ba835-4537-485f-b623-5564d1c57338
floor.(xx)

# ╔═╡ 66da211b-0e53-41c1-9ce4-75fb95a35991
begin
	n_trials = 1
	n_boots = 1e1
	bin_counts = collect( 10:2:50)
	subsamps = [0.50, 0.75, 0.80, 0.85, 0.90, 0.95]
end

# ╔═╡ ba630ae6-8ab2-4ff8-bfd4-cab5f0cf0d3f
mm = convert.(Int,round.(subsamps*n_replicates))

# ╔═╡ 4d9dded4-923c-41f9-86c5-9a8982e3a5de
n_subsamp_sets = length(mm)

# ╔═╡ e7c1b228-7945-4be6-8c48-891301693b5e
replicate_idcs = 1:n_replicates

# ╔═╡ f3e394f0-401f-41d9-a667-445a3c3d2329
n_bin_sizes = length(bin_counts)

# ╔═╡ 22d52f0f-64cf-48fb-8895-41f504ef2da3
# perm_order = ... # For use with n-D hist

# ╔═╡ 9c2130eb-031c-4fb4-b983-075f6af2ca7b
md"""
✅ P_g(g)
"""

# ╔═╡ 0eb04b09-02d4-4151-9f05-a306eac00719
g = gene_data_array[100:900,:,1][:];

# ╔═╡ 9d1678e1-079e-469a-b9f7-656ae688a628
h_hb = normalize(fit(Histogram,g,nbins=40),mode=:probability)

# ╔═╡ a915181f-4030-4a86-88ed-df804d970a02
length(range(minimum(g),stop=maximum(g),length=29))-1

# ╔═╡ 89ab03c0-1906-4ca2-a7f2-ac3c0d0ccea1
plot(h_hb)

# ╔═╡ d56fd8e0-74b0-40e5-b0ef-b8a7ed54432e
begin
	xxx = [1, 1, 1, 1, 1, 1]
	h = fit(Histogram, x, range(minimum(xxx), stop=maximum(xxx), length=4))
end

# ╔═╡ e8fc5beb-989f-4486-a798-45b34b4a0753
sum(h_hb.weights)

# ╔═╡ 5d01e544-05ce-4df7-982b-265952bc9bdd
size(h_hb.weights)

# ╔═╡ fd7d9cbe-6b30-4187-9831-b55053600b69
md"""
✅ Joint P(g,x)
"""

# ╔═╡ d96c3c1e-7023-4ae5-8c56-1a4c63ef4263
g1 = gene_data_array[100:900,:,1]

# ╔═╡ f1f5ae30-bbfb-4922-9754-a2f0f907d2dd
G = gene_data_array[100:900,:,:]

# ╔═╡ ca0dec61-20e8-457f-81fd-e4d7df52e6c2
x_column = repeat(xx,n_replicates)

# ╔═╡ e11a583d-d580-42ce-9266-300464761324
size(x_column)

# ╔═╡ 86f2c5a7-8930-4c03-af65-4c56670c1354
begin
	hist_data_joint = (x_column,g1[:])
	pyx_joint = normalize(fit(Histogram,hist_data_joint, nbins = 30),mode=:probability)
end

# ╔═╡ 064e1fa8-d714-4372-82c4-a5197e11675c
plot(pyx_joint)

# ╔═╡ 113c7c76-b732-45e6-a45c-c2156b5c80fe
md"""
✅ Random subsample of replicates
"""

# ╔═╡ ff14e4b5-9527-4132-8bb4-e6ae76536b00
k = mm[2]

# ╔═╡ 92a21016-f530-4987-beb8-b8ac9cad91a6
# xy-align with test data
begin
	
	# 👇 This is where I don't know what I'm doing 👇
	
	# I don't know how to pass the variables and parameters,
	# and I don't know how to access them from within the function.
	function xy_objective(α...)
	
	# 👆 This is where I don't know what I'm doing 👆
		

		n_p, n_r, n_t = size(Y_0)
		n_nans = convert(Int,round(nan_pad_fraction*n_p))

		# 1. circshift by γ
		nan_pad = replace!(zeros(n_nans,n_r,n_t),0=>NaN)
		Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
		Y_0_shift = Y_0

		for i_r = 1:n_r
			Y_0_shift[:,i_r,:] = circshift(Y_0_padded[:,i_r,:],γ[i_r])
		end

		# 2. take nanmean of each target to reduce variance against target mean
		ȳ = reshape(nanmean(Y_0_shift,dims=2),n_p + 2*n_nans, n_t)
		
		# 3. calculate cost
		χ²_t = zeros(n_t)
		for i_t = 1:n_t
			χ²_t[i_t] = 
			sum(Y_0_shift[:,:,i_t] - 
				(reshape(α,(1,n_r)) .+ ȳ[:,i_t]*reshape(β,(1,n_r))).^2)
		end		

		χ² = sum(χ²_t)

		return χ²
	end
	
	# Sample data
	n_pts = 1000 # number of data sample points per experiment
	n_reps = 10 # number of experimental replicates
	n_targ = 2 # number of experimental targets (simultaneously acquired in each rep)
	
	Y_0 = rand(n_pts, n_reps, n_targ) # data matrix
	
	# Goal: reduce total variance by:
	# 	adding α to each column
	# 	scaling scaling by β for each column
	# 	shifting by γ for each target
	
	# Pad edges with NaN values so data can shift 
	np_frac = 0.1 
	# e.g. a vector with 100 elements along the first dimension will have 100*np_frac NaNs added to beginning and end to become 100*(1+np_frac) elements long in the first dimension
	
	# Build model
	m = Model(Ipopt.Optimizer)
	
	register(m, 
		:xy_objective, 
		n_reps*(2*n_targ + 1), 
		xy_objective; autodiff=true)

	α_test, β_test, γ_test = nothing, nothing, nothing
	nan_pad_fraction, Y_0_prob = nothing, nothing
		
	# α: additive variable, 1 for each replicate for each target
	@variable(m, 
		α_test[i=1:n_reps,j=1:n_targ])
	set_start_value(α_test, zeros(n_reps,n_targ))
	# β: scaling variable, 1 for each replicate for each target
	@variable(m, 
		β_test[i=1:n_reps,j=1:n_targ])
	# γ: shifting variable, 1 for each replicate (shift all targets by same amount)
	@variable(m, 
		-convert(Int,floor(np_frac*n_pts)) <= 
		γ_test[i=1:n_reps] <= 
		convert(Int,floor(np_frac*n_pts)) ) # also need integer constraint
	
	# Fraction of data length to allow for shifting
	@NLparameter(m, 
		nan_pad_fraction == np_frac)
	# Starting point for experimental data
	@NLparameter(m, 
		Y_0_prob[i=1:n_pts,j=1:n_reps,k=1:n_targ] == Y_0[i,j,k])
	
	# 👇 This is where I don't know what I'm doing 👇

	# I don't know how to pass the variables and parameters,
	# and I don't know how to access them from within the function.
	@NLobjective(m, 
		Min, 
		xy_objective(α_test..., β_test..., γ_test..., Y_0_prob..., nan_pad_fraction))
	
	# 👆 This is where I don't know what I'm doing 👆

	optimize!(m)
	
	# α_test_opt = resha
	
	
	
# 		optimize!(model)

# 	α_opt = reshape(value.(α),(1,n_replicates))
# 	β_opt = reshape(value.(β),(1,n_replicates))

# 	Y_align = reshape(Y,(n_points,n_replicates))

# 	Y_align = (Y_align .- α_opt)./β_opt

# 	for i=1:3 println(".") end

# 	println(solution_summary(model))

# 	return Y_align, α_opt, β_opt, model
	
end

# ╔═╡ fc4483ec-d716-4279-9ca2-d827214e64b8
n_p, n_r, n_t = size(Y_0)

# ╔═╡ 0c3adb9e-8bf1-4f74-a72a-7d89f7129b13
begin
	# Align X,Y by doing align_y() for every possible shift in x.
	# nan_pad_fraction = 0.1
	
	
	# n_p, n_r, n_t = size(Y)
	# n_nans = convert(Int,round(nan_pad_fraction*n_p_gap))

	# 1. circshift by γ
	# nan_pad = replace!(zeros(n_nans,n_r,n_t),0=>NaN)
	nan_pad = fill(NaN, n_nans, n_r, n_t)
	Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
	Y_0_shift = Y_0

	for i_r = 1:n_r
		Y_0_shift[:,i_r,:] = circshift(Y_0_padded[:,i_r,:],γ[i_r])
	end

	# 2. take nanmean of each target to reduce variance against target mean
	ȳ = reshape(nanmean(Y_0_shift,dims=2),n_p + 2*n_nans, n_t)
end

# ╔═╡ 33517f40-2661-4a1b-ba7b-6f8f4d802b4b
function align_xy(Y_0)
	
	
	shift_array = round.(rand(Uniform(-n_nans,n_nans),1,11))
	Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
	Y_0_shift = Y_0_padded
	for i_r = 1:n_r
		Y_0_shift[:,i_r,:] = circshift(Y_0_padded[:,i_r,:],shift_array[i_r])
	end
	
# REFERENCE 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 
# 	model = Model(Ipopt.Optimizer)

# 	n_points, n_replicates = size(Y)

# 	α, β = nothing, nothing

# 	@variable(model, α[i=1:n_replicates])
# 	@variable(model, β[i=1:n_replicates])

# 	@objective(model, Min, sum((Y - (reshape(α,(1,n_replicates)) .+ ȳ*reshape(β,(1,n_replicates)))).^2))

# 	optimize!(model)

# 	α_opt = reshape(value.(α),(1,n_replicates))
# 	β_opt = reshape(value.(β),(1,n_replicates))

# 	Y_align = reshape(Y,(n_points,n_replicates))

# 	Y_align = (Y_align .- α_opt)./β_opt

# 	for i=1:3 println(".") end

# 	println(solution_summary(model))

# 	return Y_align, α_opt, β_opt, model
# REFERENCE 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 
	
	
end

# ╔═╡ 158a1b8b-c2e5-43fe-a05e-b49d823be9eb
subsample_idcs = sample(replicate_idcs, k, replace=false)

# ╔═╡ 74ce856f-5f00-49f3-a228-9da74e83f359
g1_subsample = g1[:,subsample_idcs]

# ╔═╡ 7189f9ac-b083-425a-b883-7c4fb350dbfe
G_subsample_temp = G[:,subsample_idcs,:] # see line 63 in midirectestimate.m

# ╔═╡ b7067b2a-ee9e-4ddf-bf22-57405585c204
# Reshape G_subsample 
# first argument: number of entries in G_subsample[:,:,1]
# second argument: n_targets

# final dims will be ~1bajillion x 4

# ╔═╡ 042c9c93-b7fd-4783-88a1-4106522d27f1
G_subsample = reshape(G_subsample_temp, length(G_subsample_temp[:,:,1]), n_targets)

# ╔═╡ bd47febd-c099-40a2-919f-52e948da0f8b
# Make an x column to correspond to the data columns (subsampled n_replicates)
x_column_sub = repeat(xx,k)

# ╔═╡ c07fec15-c32b-4ae1-857f-87892dc6eecf
size(x_column_sub)

# ╔═╡ d0a77da5-0d0b-4997-be39-82fabc24938d
begin
	# Get the joint distribution of the subsampled data and x
	
	hist_data_sub_joint = tuple(x_column_sub,[G_subsample[:,c] for c in 1:size(G_subsample,2)]...) 
	
	pyx_joint_sub = normalize(fit(Histogram,hist_data_sub_joint, nbins = (3,3,3,3,3)),mode=:probability)

end

# ╔═╡ f0537a73-c9c7-421a-82d8-32fa22f8ab72
pyx_joint_sub.weights = pyx_joint_sub.weights .+ eps()

# ╔═╡ 619a2c99-b8a9-4b00-b81f-889cdbf726fc
# Permute the order of the n-D histogram to put x at the end ... not sure why this is needed, but it's required to work in MATLAB

perm_order = vcat(collect(2:1:n_targets+1),1)

# ╔═╡ 9c369db4-ff2a-4af6-ad21-3879f7d81e3a
begin
	g2 = G[:,:,2]
	
	g1_1 = g1[:,1]/maximum(g1[:,1])
	g2_1 = g2[:,1]/maximum(g1[:,1])
end

# ╔═╡ 5d265b5e-3ed7-4ad1-b370-54ca4396fc32
begin
	using Trapz
	u1 = g1_1
	u2 = g2_1 
	trapz(u1)
end

# ╔═╡ d9399fcf-dfe5-4ca9-acbc-251a2aaf3916
sum(u1.*u2)

# ╔═╡ 69b1aa9e-8bc3-440f-a040-eb09883bb54b


# ╔═╡ f14f5199-87ee-414f-8a02-bc5f7ac58a55


# ╔═╡ 81041290-7b64-46ae-b315-1c3489c52499


# ╔═╡ 80518999-65c1-4a26-8a6c-89319430ed08


# ╔═╡ dd7b721d-fab0-4849-b305-dab72df2c73f


# ╔═╡ db4639b0-17c4-4efc-a42b-d3e51166fd26


# ╔═╡ bdb73836-cd26-4869-b8d3-bc0dac2305ed
md"""
## 🚧 Extrapolate I({g};x).
"""

# ╔═╡ f31e63ad-81c5-48f6-a572-e925141a1b97
md"""
## 🚧 Decoding maps.
"""

# ╔═╡ 90205ce4-dc4b-4b27-9e12-5f9a04e655b4
md"""
# 🚧 Functions
"""

# ╔═╡ cba983b5-1813-4c4a-b0e5-bb72f463ce39
md"""
## Main Analysis Functions
"""

# ╔═╡ 4a442a3c-df99-4149-8289-0ca3146cd1ec
md"""
### ✔ raw\_data\_frame()
"""

# ╔═╡ 5f271db7-1dca-4249-aa09-c4260c9da2e3
md"""
### ✔ uniform\_linear\_sample()
"""

# ╔═╡ 4618f68a-ba5a-48a2-ad0e-5595d71368a4
function uniform_linear_sample(data_in, n_samps)

	nodes_in = range(1, size(data_in)[1], length=length(data_in))
	nodes_eval = range(1, nodes_in[end], length = n_samps)

	int = interpolate((nodes_in,), data_in, Gridded(Linear()))

	data_out = zeros(length(nodes_eval))
	idx = trues(length(nodes_eval))

	data_out[idx] = int(nodes_eval[idx])

	return data_out

end

# ╔═╡ 9a2e237d-d913-4974-b278-968ca0f9453e
md"""
### ✔ regress\_noise()
"""

# ╔═╡ 9502d39e-a04d-4a44-8415-9e68f05031be
function regress_noise(n_ref,x)
	# Correct a signal measurement by accounting for variations in nuclear density.
	
	# Inputs
	# > n_ref: measured noise reference (nuclear channel)
	# > x: measured signal + noise
	
	# Outputs
	# s̃: approximate pure signal (mean ~ 0)
	
	β = (transpose(x)*n_ref)/(transpose(n_ref)*n_ref)
	
	s̃ = x - β*n_ref
end

# ╔═╡ 45e668ef-f791-40d1-b244-524d18f5eb56
md"""
### ✔ anchor\_mean\_0\_to\_1()
"""

# ╔═╡ 2e7252a1-5a1b-49e5-854a-ae75f082eb6d
function anchor_mean_0_to_1(data)
	# Inputs
	# > data: n_points x n_replicates
	#
	# Outputs
	# > data_anchor: array of same size as data where the minimum and maximum values 		for the mean profile is 0 and 1, respectively.
	# > data_anchor_mean: the mean of the anchored data

	n_points, n_reps = size(data)

	data_μ= nanmean(data, dims=2)
	data_μ_reshape = reshape(data_μ,n_points)

	data_μ_mins = reshape(nanminimum(data_μ_reshape, dims=1),1,1)
	data_μ_maxes = reshape(nanmaximum(data_μ_reshape, dims=1),1,1)

	data_⚓ = (data .- data_μ_mins)./
		(data_μ_maxes .- data_μ_mins)

	data_⚓_μ = nanmean(data_⚓, dims=2)
	data_⚓_μ_reshape = reshape(data_⚓_μ, n_points)

	return data_⚓, data_⚓_μ_reshape

end

# ╔═╡ d6a76d2b-f924-42e4-8a17-931baefb51c1
begin
	data_anchor = zeros(size(gene_data_array))
	data_anchor_mean = zeros(1000,102)
	for i_t = 1:4
		gene_array = reshape(gene_data_array[:,:,i_t],(1000,102))
		
		data_anchor[:,:,i_t], data_anchor_mean[:,i_t] = anchor_mean_0_to_1(gene_array)
	end
end

# ╔═╡ cbf17e2e-47a5-4390-99e6-c3f5882c3a10
md"""
### ✔ anchor\_mean\_to\_1()
"""

# ╔═╡ beb19c28-f5dd-4f36-9e55-85c6b9508024
	function anchor_mean_to_1(data)
		# Inputs
		# > data: n_points x n_replicates
		#
		# Outputs
		# > data_anchor: array of same size as data where the mean value of all profiles is anchored to 1 (for nuclear channel)
		# > data_anchor_mean: the mean of the anchored data
		
		n_points, n_reps = size(data)

		data_μ= nanmean(data, dims=2)
		data_μ_reshape = reshape(data_μ,n_points)

		data_μ_mins = reshape(nanminimum(data_μ_reshape, dims=1),1,1)
		data_μ_maxes = reshape(nanmaximum(data_μ_reshape, dims=1),1,1)

		data_⚓ = (data)./(nanmean(data_μ))

		data_⚓_μ = nanmean(data_⚓, dims=2)
		data_⚓_μ_reshape = reshape(data_⚓_μ, n_points)
		
		
		return data_⚓, data_⚓_μ_reshape
	end

# ╔═╡ 327d1ba2-226c-4fb7-9a78-6acca773322c
md"""
### smooth\_data()
"""

# ╔═╡ ff920b1e-a2b3-4c5b-8cc4-80aebbffe4b5
function smooth_data()
	
end

# ╔═╡ 4957bf80-0f7f-4ca2-997b-7c6ff6972004
md"""
### align\_y()
"""

# ╔═╡ 833c1369-b584-4161-b92a-bb9f99280af3
function align_y(Y,ȳ)
	# Using JuMP, perform a χ² variance minimization alignment using a linear model to bring individual replicate profiles closer to the sample mean.

	# Inputs
	# > Y: n_points x n_replicates 2D array
	# > ȳ: n_points x 1 2D array

	# Outpus
	# > Y_align: n_points x n_replicates optimally aligned 2D array
	# > α_opt: n_replicates additive parameter 1D array
	# > β_opt: n_replicates scaling parameter 1D array
	# > model: the full model for minimization

	model = Model(Ipopt.Optimizer)

	n_points, n_replicates = size(Y)

	α, β = nothing, nothing

	@variable(model, α[i=1:n_replicates])
	@variable(model, β[i=1:n_replicates])

	@objective(model, 
		Min, 
		sum((Y - (reshape(α,(1,n_replicates)) .+ ȳ*reshape(β,(1,n_replicates)))).^2))

	optimize!(model)

	α_opt = reshape(value.(α),(1,n_replicates))
	β_opt = reshape(value.(β),(1,n_replicates))

	Y_align = reshape(Y,(n_points,n_replicates))

	Y_align = (Y_align .- α_opt)./β_opt

	for i=1:3 println(".") end

	println(solution_summary(model))

	return Y_align, α_opt, β_opt, model

end

# ╔═╡ d381176e-7653-448a-80d2-d862031743d5
begin
	n_p_gap, n_r_gap, n_t_gap = size(data_anchor[100:900,:,:])
	Y_yalign_gap = zeros(n_p_gap, n_r_gap, n_t_gap)
	for i_t = 1:4
		Y_temp = data_anchor[100:900,:,i_t]
		ȳ_temp = data_anchor_mean[100:900,i_t]

		Y_yalign_gap[:,:,i_t], α_opt_gap, β_opt_gap, model_gap = align_y(Y_temp, ȳ_temp)
	end
end

# ╔═╡ 1daee640-7efc-448f-9477-5cac9c5ae9d0
md"""
### 👀 align\_xy()
"""

# ╔═╡ ff9bb147-a271-435b-9cbd-f914b3315864
# function align_xy(Y_0)
# 	# Input
# 	# > Y_0: data to align n_points x n_replicates x n_targets
	
# 	# Output
# 	# > Y_xyalign: aligned data
# 	# > α_opt: additive parameter (1 per profile per replicate)
# 	# > β_opt: scaling parameter (1 per profile per replicate)
# 	# > γ_opt: left-right shift parameter (1 per replicate, all targets within a replicate are shifted by the same amount)
	
# 	n_points, n_replicates, n_targets = size(Y_0)
	
# 	# Pad data with NaN values for shifting indices.
# 	nan_pad_percent = 0.1

# 	n_nans = convert(Int64,round(nan_pad_percent*n_points))	
# 	nan_pad = replace!(zeros(n_nans,n_replicates), 0=>NaN)
# 	Y_0_padded = vcat(nan_pad,Y_0,nan_pad)
	
# 	model = Model(Ipopt.Optimizer)
	
# 	α, β, γ = nothing, nothing, nothing
	
# 	@variable(model, α[i=1:n_replicates,j=1:n_targets]) # n_targest x n_replicates
# 	@variable(model, β[i=1:n_replicates,j=1:n_targets]) # n_targets x n_replicates
# 	@variable(model, γ[i=1:n_replicates],Int) # n_replicates - integer array
	
	
# 	@objective(model, Min, xy_objective(α,β,γ,Y_0))
	
# 	function xy_objective(α,β,γ,Y_0)
		
# 		# 1. circshift by γ
# 		circshift(Y,(γ,0))
# 		# 2. take nanmean of each target to reduce variance against
		
		
		
# 		return χ²
# 	end
	
# # REFERENCE 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 👇 
# # 	model = Model(Ipopt.Optimizer)

# # 	n_points, n_replicates = size(Y)

# # 	α, β = nothing, nothing

# # 	@variable(model, α[i=1:n_replicates])
# # 	@variable(model, β[i=1:n_replicates])

# # 	@objective(model, Min, sum((Y - (reshape(α,(1,n_replicates)) .+ ȳ*reshape(β,(1,n_replicates)))).^2))

# # 	optimize!(model)

# # 	α_opt = reshape(value.(α),(1,n_replicates))
# # 	β_opt = reshape(value.(β),(1,n_replicates))

# # 	Y_align = reshape(Y,(n_points,n_replicates))

# # 	Y_align = (Y_align .- α_opt)./β_opt

# # 	for i=1:3 println(".") end

# # 	println(solution_summary(model))

# # 	return Y_align, α_opt, β_opt, model
# # REFERENCE 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 👆 
	
	
# end

# ╔═╡ 950d747a-d37a-40f2-b872-624f54419de5
md"""
## Helper Functions
"""

# ╔═╡ 64104b9b-cdec-4f8f-a114-c9759d5e5fa3
md"""
### ✔ path\_matches\_id()
"""

# ╔═╡ 393f89d0-2a08-4f34-918a-f6a80b5a278f
function path_matches_id(id, channel, test_string)
# Check if a path contains a specific cochlea ID.

	pattern_string = "MAX_($id)_($channel)"
	match_obj = match(Regex(pattern_string), test_string)

	if match_obj == nothing
		return false
	else
		return true
	end
end

# ╔═╡ fef63029-1df9-45cd-81a0-4a511948c358
function raw_data_frame(prof_paths,img_paths)
### Get the data into data frames.

	df_data = DataFrame(id=[],
		x=[],
		topro3=[], sox2=[], psmad=[])

	df_img = DataFrame(topro3_image_data=[], sox2_image_data=[], psmad_image_data=[])

	for prof_path in prof_paths
		df_data_csv = DataFrame(CSV.File(prof_path))

		df_img_temp = DataFrame()

		id_obj = match(r"\d{2}\_.+(?=.csv)", prof_path)
		id = id_obj.match

		for img_path in img_paths

			# Get topro3 image for current cochlea
			tf_topro3 = path_matches_id(id, "topro3", img_path)
			tf_sox2 = path_matches_id(id, "sox2", img_path)
			tf_psmad = path_matches_id(id, "psmad", img_path)

			if tf_topro3
				insertcols!(df_img_temp, 1, :topro3_image_data => [load(img_path)])
			elseif tf_sox2
				insertcols!(df_img_temp, 1, :sox2_image_data => [load(img_path)])
			elseif tf_psmad
				insertcols!(df_img_temp, 1, :psmad_image_data => [load(img_path)])
			end

		end

		push!(df_data,[id, df_data_csv.x, df_data_csv.topro3, df_data_csv.sox2,df_data_csv.psmad])

		push!(df_img,
			[df_img_temp.topro3_image_data, 
			 df_img_temp.sox2_image_data,
			 df_img_temp.psmad_image_data])

		# WIP: add ROI column

		global df = hcat(df_data,df_img)

	end
	return df
end

# ╔═╡ 46fbc4fe-d1dd-4447-917b-031bdb129612
### Get the data into dataframs.
begin
	df_e12 = raw_data_frame(e12_prof_paths, e12_img_paths);
	df_e13 = raw_data_frame(e13_prof_paths, e13_img_paths);
end

# ╔═╡ c0c0e2e0-5bc9-435e-a73c-d5a0e69c9640
df_e12

# ╔═╡ 1087a093-3966-450c-b8fe-1c28039c4769
begin
	n_e12 = nrow(df_e12)
	n_e13 = nrow(df_e13)
	
	n_samples = 1000
	
	raw_data = ["x","topro3","sox2","psmad"]
	n_raw_data = length(raw_data)
	
	# can this be condensed with broadcasting???
	for i_t = 1:n_raw_data		
		raw_col = string(raw_data[i_t])
		new_col = string(raw_data[i_t],"_usamp")
		
		df_e12[:,new_col] = uniform_linear_sample.(df_e12[:,raw_col], n_samples)
		
		df_e13[:,new_col] = uniform_linear_sample.(df_e13[:,raw_col], n_samples)
	end
end

# ╔═╡ c4933e11-8c4c-4e5d-b9f5-975266a0c229
df_e12

# ╔═╡ 00008a70-45ca-42e7-9472-d1353f473b26
begin
	noise_col = string("topro3_usamp")
	signal_targets = ["topro3","sox2","psmad"]
	
	for t in signal_targets
		signal_col = string(t,"_usamp")
		est_col = string(t,"_regress")
		
		df_e12[:,est_col] = regress_noise.(df_e12[:,noise_col],df_e12[:,signal_col])
		df_e13[:,est_col] = regress_noise.(df_e13[:,noise_col],df_e13[:,signal_col])
	end
end

# ╔═╡ 823d5226-79e6-40c9-94db-8df0d75a662e
begin
	row_idx = 1
	plot(df_e12[row_idx,:topro3_usamp])
	plot!(df_e12[row_idx,:topro3_regress])
end
			

# ╔═╡ 770a6c9d-2ece-4272-8d7f-df8f003ccbe3
df_e12

# ╔═╡ 74a9c8ec-0d6d-4d13-b010-d0a97676bc1a
plot(df_e12[:,:sox2_⚓])

# ╔═╡ 47fd720b-2374-46a9-9208-f9bb9bc3a15e
begin
	targets = ["topro3","sox2","psmad"]
	n_targets_cochlea = length(targets)
	
	df_e12_⚓ = df_e12
	df_e13_⚓ = df_e13
	
	for t in targets
		usamp_col = string(t,"_usamp")

		d_12 = df_e12[:,usamp_col]
		d_13 = df_e13[:,usamp_col]
		
		data_12 = hcat(d_12...)
		data_13 = hcat(d_13...)
		
		if t == "topro3"
			data_⚓_12, data_⚓_μ_12 = anchor_mean_to_1(data_12)
			data_⚓_13, data_⚓_μ_13 = anchor_mean_to_1(data_13)
		else
			data_⚓_12, data_⚓_μ_12 = anchor_mean_0_to_1(data_12)
			data_⚓_13, data_⚓_μ_13 = anchor_mean_0_to_1(data_13)
		end
		
		# Next, get each column of data in data_⚓ into a row of df_e12
		n_pts_12, n_profs_12 = size(data_⚓_12)
		n_pts_13, n_profs_13 = size(data_⚓_13)

		new_col = string(t,"_usamp_⚓")
		
		df_temp_12 = DataFrame(new_col => [])
		df_temp_13 = DataFrame(new_col => [])

		for i_p = 1:n_profs_12
			push!(df_temp_12,[data_⚓_12[:,i_p]])
		end	
		for i_p = 1:n_profs_13
			push!(df_temp_13,[data_⚓_13[:,i_p]])
		end	
		
		# Can't get this to add to the global df_e12 ... don't know why
		df_e12_⚓ = hcat(df_e12_⚓, df_temp_12)
		df_e13_⚓ = hcat(df_e13_⚓, df_temp_13)
		
	end
end

# ╔═╡ 6700db87-4976-4b71-961a-530e6dfea489
begin
	for t in targets
		
		if t === "topro3"
			col = string(t,"_usamp")
			
			# Convert from a df column of arrays to a matrix first.
			d_12 = hcat(df_e12[:,col]...)
			d_⚓_12,_ = anchor_mean_to_1(d_12)
	
			d_13 = hcat(df_e13[:,col]...)
			d_⚓_13,_ = anchor_mean_to_1(d_13)
		else
			col = string(t,"_regress")
			
			d_12 = hcat(df_e12[:,col]...)
			d_⚓_12,_ = anchor_mean_0_to_1(d_12)
			
			d_13 = hcat(df_e13[:,col]...)
			d_⚓_13,_ = anchor_mean_0_to_1(d_13)
		end
		_, n_profs_12 = size(d_⚓_12)
		_, n_profs_13 = size(d_⚓_13)
		
		# Convert the matrix into a temp df to add back to the original df.
		new_col = string(t,"_⚓")
		df_temp_12 = DataFrame(new_col => [])
		df_temp_13 = DataFrame(new_col => [])
		
		for i_p = 1:n_profs_12
			push!(df_temp_12,[d_⚓_12[:,i_p]])
		end	
		for i_p = 1:n_profs_13
			push!(df_temp_13,[d_⚓_13[:,i_p]])
		end	
		
		df_e12[:,new_col] = df_temp_12[:,new_col]
		df_e13[:,new_col] = df_temp_13[:,new_col]
		
	end
end

# ╔═╡ bbbb002a-2098-43ce-8f3e-ebced0ef48f7
df_e12_⚓

# ╔═╡ 58fef615-a305-48a2-9daa-e203a139ea6a
df_e13_⚓

# ╔═╡ 9fe1f2a9-4d87-42b6-ac3d-1d7282c836eb
targets

# ╔═╡ 67f0b3b7-4201-42fd-86ca-16c49d9567fd
begin
	df_e12_align = df_e12_⚓
	df_e13_align = df_e13_⚓
	
	for t in targets		
		⚓_col = string(t,"_usamp_⚓")

		d_12 = hcat(df_e12_⚓[:,⚓_col]...)			
		d_13 = hcat(df_e13_⚓[:,⚓_col]...)
		
		loc = joinpath("..","data","exp_pro")
		
		if t == "psmad"
			matwrite(joinpath(loc,"julia-data-psmad.mat"),Dict("p12" => d_12, "p13" => d_13))
		elseif t == "sox2"
			matwrite(joinpath(loc,"julia-data-sox2.mat"),Dict("s12" => d_12, "s13" => d_13))
		elseif t == "topro3"
			matwrite(joinpath(loc,"julia-data-topro3.mat"),Dict("t12" => d_12, "t13" => d_13))
		end


		d_mean_12 = mean(d_12,dims=2)
		d_mean_13 = mean(d_13,dims=2)

		d_yalign_12,_,_,_ = align_y(d_12, d_mean_12)
		d_yalign_13,_,_, _ = align_y(d_13, d_mean_13)
		
		if t == "psmad"
			matwrite(joinpath(loc,"julia-data-aligned-psmad.mat"),Dict("p12_al" => d_yalign_12, "p13_al" => d_yalign_13))
		elseif t == "sox2"
			matwrite(joinpath(loc,"julia-data-aligned-sox2.mat"),Dict("s12_al" => d_yalign_12, "s13_al" => d_yalign_13))
		elseif t == "topro3"
			matwrite(joinpath(loc,"julia-data-aligned-topro3.mat"),Dict("t12_al" => d_yalign_12, "t13_al" => d_yalign_13))
		end

		# Convert the matrix into a temp df to add back to the original df.
		_, n_profs_12 = size(d_yalign_12)
		_, n_profs_13 = size(d_yalign_13)

		new_col = string(t,"_yalign")

		df_temp_12 = DataFrame(new_col => [])
		df_temp_13 = DataFrame(new_col => [])

		for i_p = 1:n_profs_12
			push!(df_temp_12,[d_yalign_12[:,i_p]])
		end
		for i_p = 1:n_profs_13
			push!(df_temp_13,[d_yalign_13[:,i_p]])
		end
		
		df_e12[:,new_col] = df_temp_12[:,new_col]
		df_e13[:,new_col] = df_temp_13[:,new_col]

	end
end

# ╔═╡ a20f2a88-9b12-4b66-a6aa-480e9bc038ac
df_e12

# ╔═╡ 7f8a3758-5d41-4d1d-a03e-70a0ab9dd0a4
df_e13

# ╔═╡ 3ac6aef2-ff95-4828-835e-cb5e0ae8f15a
begin	
	# Plot results to check.
	
	plot_dict_12 = Dict("topro3"=>plot(title="E12.5 Nuclear Signal"),"sox2"=>plot(title="Sox2"),"psmad"=>plot(title="P-Smad"))
	plot_dict_13 = Dict("topro3"=>plot(title="E13.5 Nuclear Signal"),"sox2"=>plot(title="Sox2"),"psmad"=>plot(title="P-Smad"))

	
	for t in targets
		raw_col_12 = string(t,"_⚓")
		aln_col_12 = string(t,"_yalign")
		
		raw_data_12 = hcat(df_e12[:,raw_col_12]...)
		aln_data_12 = hcat(df_e12[:,aln_col_12]...)
		
		mean_raw_12 = mean(raw_data_12,dims=2)
		std_raw_12 = std(raw_data_12,dims=2)
		
		mean_aln_12 = mean(aln_data_12,dims=2)
		std_aln_12 = std(aln_data_12,dims=2)
		
		plot!(plot_dict_12[t], 
		mean_raw_12, linewidth=2, ribbon= std_raw_12,
		color=1)
		plot!(plot_dict_12[t], 
		mean_aln_12, linewidth=2, ribbon= std_aln_12,
		color=2)
		
		
		raw_col_13 = string(t,"_⚓")
		aln_col_13 = string(t,"_yalign")
		
		raw_data_13 = hcat(df_e13[:,raw_col_13]...)
		aln_data_13 = hcat(df_e13[:,aln_col_13]...)
		
		mean_raw_13 = mean(raw_data_13,dims=2)
		std_raw_13 = std(raw_data_13,dims=2)
		
		mean_aln_13 = mean(aln_data_13,dims=2)
		std_aln_13 = std(aln_data_13,dims=2)
		
		plot!(plot_dict_13[t], 
		mean_raw_13, linewidth=2, ribbon= std_raw_13,
		color=3)
		plot!(plot_dict_13[t], 
		mean_aln_13, linewidth=2, ribbon= std_aln_13,
		color=4)
		
	end
end

# ╔═╡ 05373280-3e42-4eb3-847c-e9a931909db4
plot(plot_dict_12["topro3"],
	plot_dict_12["sox2"],
	plot_dict_12["psmad"],
	layout=(1,3),legend=false,ylims=(-0.1,1.3))

# ╔═╡ 144a99c9-e3cb-44c2-8b77-6d685f0773d7
begin
	plot(plot_dict_12["sox2"],
		plot_dict_12["psmad"], 
		layout=(1,2),legend=false,ylims=(-0.1,1.3))
end

# ╔═╡ 80c62ecc-f612-4975-8392-0fb5c6a7c18e
plot(plot_dict_13["sox2"],
	plot_dict_13["psmad"],
	layout=(1,2),legend=false,ylims=(-0.1,1.3))

# ╔═╡ b066c1cb-5121-4a3e-bf51-f31844b7dda3
# ✅ Build the 3D matrix for alignment
begin	
	i_t = 0
	for t in targets[2:3] # Exclude topro3 ... only shift it after.
		i_t = i_t + 1
		yalign_col = string(t,"_yalign")
		d_12 = hcat(df_e12[:,yalign_col]...)
		Y_0[:,:,i_t] = d_12 # working
	end
end

# ╔═╡ be2d2b73-c63c-43e0-89d2-ebef2cd0e6c9


# ╔═╡ ce817927-c6ec-4d9b-84da-11b34eb08022


# ╔═╡ Cell order:
# ╠═a4dd3c3d-26a8-4c55-8dde-3c7f6301f01e
# ╠═485d8da0-be94-11eb-183b-e5e105e25b55
# ╠═ea9907cf-6f0e-495e-a172-bc6fd0cf2863
# ╠═4f275ac2-e986-4eec-bcef-fd7178e96d3d
# ╟─f0638606-63c0-4541-bd76-2887002557f0
# ╟─8b0d2a11-ca25-4d95-bb74-33a5a8f5b5e5
# ╟─69573234-68a8-4e41-bd34-7c5bb013e7f9
# ╟─ed292325-e42e-43af-936f-6a98b1c86e3f
# ╟─d4e4a684-1b0f-45d5-a377-35f005333bcc
# ╟─b40521ca-9982-4c61-954d-920a766380d2
# ╟─097cad13-82b7-4859-adc0-4a0bba8bc632
# ╟─ab4a23bb-dd9e-4455-ba6a-08e007f7092e
# ╟─51e30cad-30b7-44ff-adf3-d06018541e16
# ╠═6d6bfc47-b00d-4575-a62f-ac99fd37e641
# ╠═b165c538-a5f5-4ea6-8ce5-6d60943b3cfa
# ╠═7456dc7a-fed8-4134-ba4c-f05ad99665f4
# ╟─19fa8dc3-fc71-4cfa-87c7-961d94fc2ade
# ╟─ca5b9a0e-59b4-478d-9347-f500a364a18d
# ╠═c8a39986-acb9-41ca-bb0b-ad5d005d4c08
# ╠═335bae5d-d0b9-4620-9ec9-71e427bf5dea
# ╠═0b49770b-a6cb-44bf-a1e1-f9db3e361cb5
# ╠═e64c4b85-3550-43d2-bdf6-8979ff5fcb27
# ╟─8cf2b1e5-a403-4e69-aabd-d712b3b746ac
# ╠═1fb671db-4ee7-4756-87ab-7970a701e7bb
# ╠═e56f2074-5f94-43ba-9b10-662f65723033
# ╟─7123c8bc-7e48-4d7a-91dc-d35885407ac7
# ╟─05d085bb-878d-4e93-9985-373e9955892a
# ╠═d6a76d2b-f924-42e4-8a17-931baefb51c1
# ╟─ef7a4517-be6e-46ab-b85f-f3254e61905a
# ╠═d381176e-7653-448a-80d2-d862031743d5
# ╠═0c3adb9e-8bf1-4f74-a72a-7d89f7129b13
# ╠═7170e92a-ab07-4ee9-aad3-dc30783f016e
# ╟─b9c644e8-86c4-4935-8117-bfde3173b0f7
# ╟─6cb4b869-dcd0-414b-88c5-67bf69d084a2
# ╟─3e9036f5-5d11-4950-bba8-a2acfb84bf5e
# ╠═9fabf548-d740-4989-a0e6-4700f8436cfe
# ╟─b695d442-89e2-4f21-aca6-52d26ca305c2
# ╠═46fbc4fe-d1dd-4447-917b-031bdb129612
# ╠═c0c0e2e0-5bc9-435e-a73c-d5a0e69c9640
# ╟─4c86e3ce-cc11-4644-88d9-c52cb41d09e3
# ╠═1087a093-3966-450c-b8fe-1c28039c4769
# ╠═c4933e11-8c4c-4e5d-b9f5-975266a0c229
# ╠═ef16f350-ab20-4aa1-ae2f-d68f0d6f8f00
# ╟─a80ca2b0-c4bf-4d9d-9ead-9fe5fff2ae58
# ╠═00008a70-45ca-42e7-9472-d1353f473b26
# ╠═823d5226-79e6-40c9-94db-8df0d75a662e
# ╠═770a6c9d-2ece-4272-8d7f-df8f003ccbe3
# ╟─6895787b-70f3-4cd1-9e9f-243c22662f99
# ╠═6700db87-4976-4b71-961a-530e6dfea489
# ╠═74a9c8ec-0d6d-4d13-b010-d0a97676bc1a
# ╠═47fd720b-2374-46a9-9208-f9bb9bc3a15e
# ╠═bbbb002a-2098-43ce-8f3e-ebced0ef48f7
# ╠═58fef615-a305-48a2-9daa-e203a139ea6a
# ╟─8e6682fa-6790-4275-9ad2-3aa9c1c72c1c
# ╠═c3c4517b-aca3-437f-a3ca-8653b91b1e8a
# ╟─1689ecbf-0750-46c1-82f2-d0e598d15448
# ╠═9fe1f2a9-4d87-42b6-ac3d-1d7282c836eb
# ╠═67f0b3b7-4201-42fd-86ca-16c49d9567fd
# ╠═a20f2a88-9b12-4b66-a6aa-480e9bc038ac
# ╠═7f8a3758-5d41-4d1d-a03e-70a0ab9dd0a4
# ╟─3ac6aef2-ff95-4828-835e-cb5e0ae8f15a
# ╠═05373280-3e42-4eb3-847c-e9a931909db4
# ╠═144a99c9-e3cb-44c2-8b77-6d685f0773d7
# ╠═80c62ecc-f612-4975-8392-0fb5c6a7c18e
# ╟─f69ce4e9-3f7f-4928-ad8f-b8d98384565f
# ╟─3bf6395b-acf2-4396-8b1a-f65291df41ed
# ╟─8f6133fd-2508-4e06-97a9-65985b779a1c
# ╠═92a21016-f530-4987-beb8-b8ac9cad91a6
# ╠═92bceb54-b2a0-4d82-a2d4-f3770c141834
# ╠═158ba835-4537-485f-b623-5564d1c57338
# ╠═bc001ccd-321b-4a9c-80b9-9bf36e79335a
# ╠═b066c1cb-5121-4a3e-bf51-f31844b7dda3
# ╠═fc4483ec-d716-4279-9ca2-d827214e64b8
# ╠═1c0ec19a-4229-43e3-83fd-2c07c03b8e44
# ╠═a1965e64-9449-47b2-a706-ebe36080f2a3
# ╠═c5cdaccd-00b6-4fa8-b2ad-294feead5c90
# ╠═1934c23e-478f-4404-ad4c-b37c33f526a2
# ╠═d9edc6c1-07cd-417e-9f81-cb5db71c5185
# ╠═ed85ac59-4332-4c76-8e35-a1e193d5d185
# ╠═7420d1ae-dd36-478e-96c8-f3c65a3746ce
# ╠═33517f40-2661-4a1b-ba7b-6f8f4d802b4b
# ╟─45cf206d-7405-464c-a86d-a4cfa8af3119
# ╠═9f64c871-04a4-4d33-ac5d-1bc49b4d0485
# ╠═64b5deca-ab94-4f74-b960-711d2df20894
# ╟─16a12f25-b06a-496c-bd66-c57ff3e15783
# ╠═0bf8238e-4217-45cb-bb78-28d6ba420a8e
# ╠═66da211b-0e53-41c1-9ce4-75fb95a35991
# ╠═ba630ae6-8ab2-4ff8-bfd4-cab5f0cf0d3f
# ╠═4d9dded4-923c-41f9-86c5-9a8982e3a5de
# ╠═e7c1b228-7945-4be6-8c48-891301693b5e
# ╠═f3e394f0-401f-41d9-a667-445a3c3d2329
# ╠═22d52f0f-64cf-48fb-8895-41f504ef2da3
# ╟─9c2130eb-031c-4fb4-b983-075f6af2ca7b
# ╠═0eb04b09-02d4-4151-9f05-a306eac00719
# ╠═9d1678e1-079e-469a-b9f7-656ae688a628
# ╠═a915181f-4030-4a86-88ed-df804d970a02
# ╠═89ab03c0-1906-4ca2-a7f2-ac3c0d0ccea1
# ╠═d56fd8e0-74b0-40e5-b0ef-b8a7ed54432e
# ╠═e8fc5beb-989f-4486-a798-45b34b4a0753
# ╠═5d01e544-05ce-4df7-982b-265952bc9bdd
# ╟─fd7d9cbe-6b30-4187-9831-b55053600b69
# ╠═d96c3c1e-7023-4ae5-8c56-1a4c63ef4263
# ╠═f1f5ae30-bbfb-4922-9754-a2f0f907d2dd
# ╠═ca0dec61-20e8-457f-81fd-e4d7df52e6c2
# ╠═e11a583d-d580-42ce-9266-300464761324
# ╠═86f2c5a7-8930-4c03-af65-4c56670c1354
# ╠═064e1fa8-d714-4372-82c4-a5197e11675c
# ╟─113c7c76-b732-45e6-a45c-c2156b5c80fe
# ╠═ff14e4b5-9527-4132-8bb4-e6ae76536b00
# ╠═158a1b8b-c2e5-43fe-a05e-b49d823be9eb
# ╠═74ce856f-5f00-49f3-a228-9da74e83f359
# ╠═7189f9ac-b083-425a-b883-7c4fb350dbfe
# ╠═b7067b2a-ee9e-4ddf-bf22-57405585c204
# ╠═042c9c93-b7fd-4783-88a1-4106522d27f1
# ╠═bd47febd-c099-40a2-919f-52e948da0f8b
# ╠═c07fec15-c32b-4ae1-857f-87892dc6eecf
# ╠═d0a77da5-0d0b-4997-be39-82fabc24938d
# ╠═f0537a73-c9c7-421a-82d8-32fa22f8ab72
# ╠═619a2c99-b8a9-4b00-b81f-889cdbf726fc
# ╠═9c369db4-ff2a-4af6-ad21-3879f7d81e3a
# ╠═5d265b5e-3ed7-4ad1-b370-54ca4396fc32
# ╠═d9399fcf-dfe5-4ca9-acbc-251a2aaf3916
# ╠═8d62b88b-f77f-4d8d-8c6a-ada1d1626b8c
# ╠═69b1aa9e-8bc3-440f-a040-eb09883bb54b
# ╠═f14f5199-87ee-414f-8a02-bc5f7ac58a55
# ╠═81041290-7b64-46ae-b315-1c3489c52499
# ╠═80518999-65c1-4a26-8a6c-89319430ed08
# ╠═dd7b721d-fab0-4849-b305-dab72df2c73f
# ╠═db4639b0-17c4-4efc-a42b-d3e51166fd26
# ╠═bdb73836-cd26-4869-b8d3-bc0dac2305ed
# ╠═f31e63ad-81c5-48f6-a572-e925141a1b97
# ╟─90205ce4-dc4b-4b27-9e12-5f9a04e655b4
# ╟─cba983b5-1813-4c4a-b0e5-bb72f463ce39
# ╟─4a442a3c-df99-4149-8289-0ca3146cd1ec
# ╠═fef63029-1df9-45cd-81a0-4a511948c358
# ╟─5f271db7-1dca-4249-aa09-c4260c9da2e3
# ╠═4618f68a-ba5a-48a2-ad0e-5595d71368a4
# ╟─9a2e237d-d913-4974-b278-968ca0f9453e
# ╠═9502d39e-a04d-4a44-8415-9e68f05031be
# ╟─45e668ef-f791-40d1-b244-524d18f5eb56
# ╠═2e7252a1-5a1b-49e5-854a-ae75f082eb6d
# ╟─cbf17e2e-47a5-4390-99e6-c3f5882c3a10
# ╠═beb19c28-f5dd-4f36-9e55-85c6b9508024
# ╟─327d1ba2-226c-4fb7-9a78-6acca773322c
# ╠═ff920b1e-a2b3-4c5b-8cc4-80aebbffe4b5
# ╟─4957bf80-0f7f-4ca2-997b-7c6ff6972004
# ╠═833c1369-b584-4161-b92a-bb9f99280af3
# ╟─1daee640-7efc-448f-9477-5cac9c5ae9d0
# ╠═ff9bb147-a271-435b-9cbd-f914b3315864
# ╟─950d747a-d37a-40f2-b872-624f54419de5
# ╟─64104b9b-cdec-4f8f-a114-c9759d5e5fa3
# ╠═393f89d0-2a08-4f34-918a-f6a80b5a278f
# ╠═be2d2b73-c63c-43e0-89d2-ebef2cd0e6c9
# ╠═ce817927-c6ec-4d9b-84da-11b34eb08022