# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python (base)
#     language: python
#     name: base
# ---

# # Translation of data import from R

# The data we are importing and formatting is for hyperspectral data and spectra (which we also fit). I translate several R functions here. This is preparation for running regression in Python (on a GPU, for faster runtimes than we've been able to get in R).

import tensorflow as tf

import spectral.io.envi as envi

import glob

import pandas as pd

import numpy as np



from rpy2.robjects import numpy2ri

files_list = glob.glob("/scratch2/NSF_GWAS/macroPhor_Array/T16_DEV_genes/EA/wk7/*.hdr")

spectral_library_path = "/scratch2/NSF_GWAS/GMOdetectoR/spectra_library/"

i = 2

file_path = files_list[i]

file_path

img = envi.open(file_path)

img.shape

type(img)

# ##### Get wavelengths

h = envi.read_envi_header(file_path)
wavelengths = h['wavelength']

# ##### Read in spectra

spectral_library = ['DsRed', 'ZsYellow', 'Chl', 'Diffraction']

j = 0

spectra_path = (spectral_library_path + spectral_library[j] + '.csv')

spectra_path

pub_emission_spectrum = pd.read_csv(spectra_path)

pub_emission_spectrum = pub_emission_spectrum.rename(columns={'emission wavelength (nm)': 'wavelength'})
pub_emission_spectrum = pub_emission_spectrum.rename(columns={'Normalized emission': 'intensity'})

#import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
#ts=robjects.r('ts')
stats=importr('stats')

import rpy2.robjects as ro
from rpy2.robjects.conversion import localconverter
from rpy2.robjects import pandas2ri

scales=importr('scales')

import rpy2
print(rpy2.__version__)

# +
# Convert pandas object to something R-friendly https://rpy2.github.io/doc/latest/html/pandas.html
with localconverter(ro.default_converter + ro.pandas2ri.converter):
  r_from_pd_df = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2=2.9.4 may have error here due to no attribute 'py2rpy'

# Use rpy2 so we can use loess from R, as in old version, get same exact results
fit = stats.loess(ro.Formula('intensity~wavelength'), data=r_from_pd_df, span=0.1)

# +
predict_from = {'wavelength':wavelengths}
predict_from = pd.DataFrame(predict_from)

with localconverter(ro.default_converter + ro.pandas2ri.converter):
  r_predict_from = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2=2.9.4 may have error here due to no attribute 'py2rpy'


# -

predictions = stats.predict(fit, r_predict_from)

max(predictions)

predictions = scales.rescale(predictions, center='FALSE') # by default, scales between 0 and 1.

scaled_emission_spectra = {'wavelength': wavelengths, 'intensity': predictions}
scaled_emission_spectra = pd.DataFrame(scaled_emission_spectra)

predictions

# Why are these 0.000201 instead of 0 after scaling?

# +
# print(predictions)
# -

# Because the 0 is in the middle. Will get rid of noise outside spectra with a threshold.

spectra_noise_threshold = float(0.001)

type(spectra_noise_threshold)

type(scaled_emission_spectra)

scaled_emission_spectra.loc[scaled_emission_spectra['intensity'] < spectra_noise_threshold, 'intensity'] = 0

# Did it work?

min(scaled_emission_spectra['intensity'])

# Yes. Min is now 0 instead of 0.000201

scaled_emission_spectra


# ##### Pack up into read_plot_spectra function equivalent to R version (but without integration of plot support yet)

def read_fit_spectra(spectra_path, wavelengths, plot=False, spectra_noise_threshold = 0.001):
    
    numpy2ri.activate()

    pub_emission_spectrum = pd.read_csv(spectra_path)
    pub_emission_spectrum = pub_emission_spectrum.rename(columns={'emission wavelength (nm)': 'wavelength'})
    pub_emission_spectrum = pub_emission_spectrum.rename(columns={'Normalized emission': 'intensity'})

    # Convert pandas object to something R-friendly https://rpy2.github.io/doc/latest/html/pandas.html
    with localconverter(ro.default_converter + ro.pandas2ri.converter):
        r_from_pd_df = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2=2.9.4 may have error here due to no attribute 'py2rpy'

    # Use rpy2 so we can use loess from R, as in old version, get same exact results
    # Fit
    fit = stats.loess(ro.Formula('intensity~wavelength'), data=r_from_pd_df, span=0.1)

    # Prepare to predict
    predict_from = {'wavelength':wavelengths}
    predict_from = pd.DataFrame(predict_from)

    with localconverter(ro.default_converter + ro.pandas2ri.converter):
        r_predict_from = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2>3.0.5 may have error here due to no attribute 'py2rpy'

    # Predict
    predictions = stats.predict(fit, r_predict_from)
    predictions = scales.rescale(predictions, center='FALSE') # by default, scales between 0 and 1.

    # Scale
    scaled_emission_spectra = {'wavelength': wavelengths, 'intensity': predictions}
    scaled_emission_spectra = pd.DataFrame(scaled_emission_spectra)

    # Denoise
    scaled_emission_spectra.loc[scaled_emission_spectra['intensity'] < spectra_noise_threshold, 'intensity'] = 0

    if plot == True:
        print('Plotting not supported here yet.')

    numpy2ri.deactivate()
        
    return(scaled_emission_spectra)


# Test the function

scaled_emission_spectra_test = read_fit_spectra(spectra_path = spectra_path,
                                                wavelengths = wavelengths,
                                                plot = False,
                                                spectra_noise_threshold = 0.001)



scaled_emission_spectra_test

# ##### Now for the whole X, built from intensity columns for each fluorophore constructed as above.

i = 0

import os

os.listdir()

# +
path = (spectral_library_path + str(spectral_library[i]) + ".csv")

spectra = read_fit_spectra(spectra_path = path,
                           wavelengths = wavelengths,
                           plot = False,
                           spectra_noise_threshold = 0.001)
# -

type(wavelengths)

type(0.01)

wavelengths

type(spectra)

spectra

intercept = 1

i

intercept



spectra['intensity']

if(i==0 and intercept==1):
    intercept_vector = [1] * len(wavelengths)
    mm = np.c_[(intercept_vector, np.asarray(spectra['intensity']))] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

intercept = 0

if(i==0 and intercept==0):
    mm = np.array(spectra['intensity'])

mm.shape

i = 1

spectra['intensity']

if(i>0):
    mm = np.c_[mm, np.array(spectra['intensity'])] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

colnames = spectral_library.copy()

colnames

spectral_library

# Did some reading and found out numpy array columns cannot be named unless you do something fancy by creating new class etc

# Moving on, let's subset X to only desired wavelengths.

min_desired_wavelength = 545
max_desired_wavelength = 722

wavelengths = np.asarray(wavelengths)

# https://stackoverflow.com/questions/13869173/numpy-find-index-of-the-elements-within-range
wavelength_indices_desired = np.where(np.logical_and(wavelengths.astype(float)>=min_desired_wavelength,
                                                     wavelengths.astype(float)<=max_desired_wavelength))

mm.shape

mm[wavelength_indices_desired,].shape



mm2 = mm[(mm["wavelength"] >= m 2) | (mm["Pclass"] == 3)]


# ##### Pack it into a function that loops over everything in the fluorophore library

def build_X(fluorophore_ID_vector, intercept, wavelengths,
            min_desired_wavelength, max_desired_wavelength):  
    
    for i in range(0, (len(fluorophore_ID_vector))):
        print(i)
        path = (spectral_library_path + str(fluorophore_ID_vector[i]) + ".csv")
        if(os.path.isfile(path) == False):
            raise NameError("Error: Spectrum for " + fluorophore_ID_vector[i] + " not found in spectra_library folder")

        spectra = read_fit_spectra(spectra_path = path,
                                   wavelengths = wavelengths,
                                   plot = False,
                                   spectra_noise_threshold = 0.001)

        if(i==0 and intercept==1):
            intercept_vector = [1] * len(wavelengths)
            mm = np.c_[(intercept_vector, np.asarray(spectra['intensity']))] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

        if(i==0 and intercept==0):
            mm = np.array(spectra['intensity'])

        if(i>0):
            mm = np.c_[mm, np.array(spectra['intensity'])] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

    wavelengths = np.asarray(wavelengths)
    
    # https://stackoverflow.com/questions/13869173/numpy-find-index-of-the-elements-within-range
    wavelength_indices_desired = np.where(np.logical_and(wavelengths.astype(float)>=min_desired_wavelength,
                                                         wavelengths.astype(float)<=max_desired_wavelength))

    return(mm[wavelength_indices_desired,])


build_X(fluorophore_ID_vector = spectral_library,
        intercept = 1,
        wavelengths = wavelengths,
        min_desired_wavelength = 545,
        max_desired_wavelength = 722)

# ##### Try again outside of the function...

fluorophore_ID_vector = spectral_library

for i in range(0, (len(fluorophore_ID_vector))):
    print(i)
    

i = 0

# +
path = (spectral_library_path + str(fluorophore_ID_vector[i]) + ".csv")
if(os.path.isfile(path) == False):
    raise NameError("Error: Spectrum for " + fluorophore_ID_vector[i] + " not found in spectra_library folder")

spectra = read_fit_spectra(spectra_path = path,
                           wavelengths = wavelengths,
                           plot = False,
                           spectra_noise_threshold = 0.001)

# +
if(i==0 and intercept==1):
    intercept_vector = [1] * len(wavelengths)
    mm = np.c_[(intercept_vector, np.asarray(spectra['intensity']))] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

if(i==0 and intercept==0):
    mm = np.array(spectra['intensity'])

if(i>0):
    mm = np.c_[mm, np.array(spectra['intensity'])] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy
# -

i = 1

# +
path = (spectral_library_path + str(fluorophore_ID_vector[i]) + ".csv")
if(os.path.isfile(path) == False):
    raise NameError("Error: Spectrum for " + fluorophore_ID_vector[i] + " not found in spectra_library folder")

spectra = read_fit_spectra(spectra_path = path,
                           wavelengths = wavelengths,
                           plot = False,
                           spectra_noise_threshold = 0.001)
# -

path = (spectral_library_path + str(fluorophore_ID_vector[i]) + ".csv")

path

if(os.path.isfile(path) == False):
    raise NameError("Error: Spectrum for " + fluorophore_ID_vector[i] + " not found in spectra_library folder")

spectra = read_fit_spectra(spectra_path = path,
                           wavelengths = wavelengths,
                           plot = False,
                           spectra_noise_threshold = 0.001)

# ##### Why doesn't read_fit_spectra work for ZsYellow?

# Work through step-by-step.

# Is the problem that DsRed library file has values for all wavelengths but ZsYellow file doesn't? This is easy to fix if so...

DsRed_spec = pd.read_csv('/scratch2/NSF_GWAS/GMOdetectoR/spectra_library/DsRed.csv')
ZsYellow_spec = pd.read_csv('/scratch2/NSF_GWAS/GMOdetectoR/spectra_library/ZsYellow.csv')

DsRed_spec.shape

ZsYellow_spec.shape

# Yes, that does seem to be the problem. Why is it happening with rpy2 but no tin R version? Work in both RStudio and Jupyter alongside step-by-step to see where it becoems 207 or 318 in each.

# That method of debugging didn't work since the R function is fed a subset of wavelengths, where here we use all wavelengths. Should do some type of merge instead.

# +
pub_emission_spectrum = pd.read_csv('/scratch2/NSF_GWAS/GMOdetectoR/spectra_library/ZsYellow.csv')
pub_emission_spectrum = pub_emission_spectrum.rename(columns={'emission wavelength (nm)': 'wavelength'})
pub_emission_spectrum = pub_emission_spectrum.rename(columns={'Normalized emission': 'intensity'})

# Convert pandas object to something R-friendly https://rpy2.github.io/doc/latest/html/pandas.html
with localconverter(ro.default_converter + ro.pandas2ri.converter):
    r_from_pd_df = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2=2.9.4 may have error here due to no attribute 'py2rpy'

# Use rpy2 so we can use loess from R, as in old version, get same exact results
# Fit
fit = stats.loess(ro.Formula('intensity~wavelength'), data=r_from_pd_df, span=0.1)
# -

# Prepare to predict
predict_from = {'wavelength':wavelengths}
predict_from = pd.DataFrame(predict_from)

predict_from

numpy2ri.activate()

with localconverter(ro.default_converter + ro.pandas2ri.converter):
    r_predict_from = ro.conversion.py2rpy(wavelengths) # If not in rpy2>3.0.5 may have error here due to no attribute 'py2rpy'

# Predict
predictions = stats.predict(fit, wavelengths)

len(predictions)


# ##### Solved. The problem was that I was predicting from original wavelengths instead of desired ones (translation error)

def build_X(fluorophore_ID_vector, intercept, wavelengths,
            min_desired_wavelength, max_desired_wavelength):  
    
    for i in range(0, (len(fluorophore_ID_vector))):
        print(i)
        path = (spectral_library_path + str(fluorophore_ID_vector[i]) + ".csv")
        if(os.path.isfile(path) == False):
            raise NameError("Error: Spectrum for " + fluorophore_ID_vector[i] + " not found in spectra_library folder")

        spectra = read_fit_spectra(spectra_path = path,
                                   wavelengths = wavelengths,
                                   plot = False,
                                   spectra_noise_threshold = 0.01)

        if(i==0 and intercept==1):
            intercept_vector = [1] * len(wavelengths)
            mm = np.c_[(intercept_vector, np.asarray(spectra['intensity']))] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

        if(i==0 and intercept==0):
            mm = np.array(spectra['intensity'])

        if(i>0):
            mm = np.c_[mm, np.array(spectra['intensity'])] # https://stackoverflow.com/questions/43961585/cbind-r-function-equivalent-in-numpy

    wavelengths = np.asarray(wavelengths)
    
    # https://stackoverflow.com/questions/13869173/numpy-find-index-of-the-elements-within-range
    wavelength_indices_desired = np.where(np.logical_and(wavelengths.astype(float)>=min_desired_wavelength,
                                                         wavelengths.astype(float)<=max_desired_wavelength))

    return(mm[wavelength_indices_desired,])


def read_fit_spectra(spectra_path, wavelengths, plot=False, spectra_noise_threshold = 0.01):
    
    numpy2ri.activate()

    pub_emission_spectrum = pd.read_csv(spectra_path)
    pub_emission_spectrum = pub_emission_spectrum.rename(columns={'emission wavelength (nm)': 'wavelength'})
    pub_emission_spectrum = pub_emission_spectrum.rename(columns={'Normalized emission': 'intensity'})

    # Convert pandas object to something R-friendly https://rpy2.github.io/doc/latest/html/pandas.html
    with localconverter(ro.default_converter + ro.pandas2ri.converter):
        r_from_pd_df = ro.conversion.py2rpy(pub_emission_spectrum) # If not in rpy2=2.9.4 may have error here due to no attribute 'py2rpy'

    # Use rpy2 so we can use loess from R, as in old version, get same exact results
    # Fit
    fit = stats.loess(ro.Formula('intensity~wavelength'), data=r_from_pd_df, span=0.1)

    # Predict
    predictions = stats.predict(fit, wavelengths)
    predictions = scales.rescale(predictions, center='FALSE') # by default, scales between 0 and 1.

    # Scale
    scaled_emission_spectra = {'wavelength': wavelengths, 'intensity': predictions}
    scaled_emission_spectra = pd.DataFrame(scaled_emission_spectra)

    # Denoise
    print(scaled_emission_spectra.head())
    print(spectra_noise_threshold)
    scaled_emission_spectra.loc[scaled_emission_spectra['intensity'] < spectra_noise_threshold, 'intensity'] = 0
    print(scaled_emission_spectra.head())
    if plot == True:
        print('Plotting not supported here yet.')

    numpy2ri.deactivate()
        
    return(scaled_emission_spectra)


build_X(fluorophore_ID_vector = spectral_library,
        intercept = 1,
        wavelengths = wavelengths,
        min_desired_wavelength = 445,
        max_desired_wavelength = 725)



StarcoderdataJupyterScriptsDedupFiltered