pyview/view/idl_translation_core/bleach_correction.py

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import logging


# for the linear fit on the logarithm, to estimate initial parameters
def fit_exp_linear_autoOffset(t, y):

    # I assume the decay has reached 90% , i.e. I subtract part of the trace
    percent_left = 0.01

    # to what level do we expect the decay to have gone?
    # e.g. percent_left = 0.1 assumes that 10% of bleachable signal is still there at the end of the measurement
    # 0.001 assumes that at the end there is almost no bleachabel signal any more, i.e. all resting fluorescence is stable

    miny, maxy = np.min(y), np.max(y)
    offset = (1 + percent_left) * miny - percent_left * maxy  # equals: np.min(y) - 0.1*(np.max(y)-np.min(y))
    # the line above fails if the signal is entirely flat - as is the case in some artificial datasets
    if offset == miny:  # i.e., miny == maxy
        return 0, 1, miny
    else:
        # y is guaranteed to be positive after offset removal (all zero case handled above)
        y = y - offset

        y_log = np.log(y)
        K, A_log = np.polyfit(t, y_log, 1)
        A = np.exp(A_log)

        return A, K, offset


# the function to be fitted in the end
def model_func(t, A, K, C):
    return A * np.exp(K * t) + C


# # debug fitlogdecay
# #create fictive time course
# A = 10
# K = -0.01
# C = 3
# length = 100
# t = np.arange(length)
# lineIn = model_func(t,A,K,C)
# lineIn = np.full((length),1000)
# #lineIn[50] = 100
# plt.plot(lineIn)
# weights = np.full((length),1)
# #test fitlogdecay
# (fittedout, opt_parms) = fitlogdecay(lineIn, weights, True, 'test')
# print(opt_parms)

def fitlogdecay(lineIn, weights, showresults=False, measurement_label=""):

    ##python help, look at:
    # https://scipython.com/book/chapter-8-scipy/examples/weighted-and-non-weighted-least-squares-fitting/


    #copy of IDL FitLogDecay
    #but in the essence, not the letter
    #lineIn is an array to fit
    #return: the fitted modelled
    #weights: same length as lineIn. 0: do not consider this value.
    # high value: point is important
    # low value:  point is not important.




#pro FitLogDecay, lineIn, fittedOut, weights, verbose ;, noConverge
#;a fitted output is given
#;Compute the fit to the function we have just defined. First, define the independent and dependent variables:
#
#common data ; necessary in order to write parameters into p structure
    ### in the python translation, parameters are not stored (yet)
#;common CFD
#;common CFDconst
#;set verbose to 1 if not explicitly stated
#IF N_PARAMS() GE 4 THEN begin
#endIF else begin
#	verbose = 1
#endElse
# replaced by keyword argument verbose=1
#Y = reform(lineIn) M reform in IDL removes empty dimensions
    y = lineIn #copy, for now, to keep close to IDL text
#X = FLOAT(INDGEN(N_elements(y)))
    t = np.arange(len(y)) #this assumes that time points are equidistant

##;nov 07, new initial estimates, follow this for python
#    A, K = fit_exp_linear(t, y, offset) #
#    print('ViewCalcData:fitlogdecay: linear fit parameters A,K are: ',A,' ', K, '. Offest is: ',offset)
#    # initial guess done without weights, without excluded values.
## A and K are the estimates

#;A = [1.0D,-0.1D,offset]	;old presets
#A = [logOffset,logdecay,offset]	;Provide an initial guess of the function's parameters.
#
#;;may 07, new approach for initial estimate - removed but kept for record
#;;a0 is offset
#;A(0) = estimate(0)
#;;a1 = (ln(Fx) - ln(a0)) / x
#;;we need the response at timepoint x.
#;fx = yfit(n_elements(y)-1)
#;vx = x(n_elements(x)-1)
#;;A(1) = (alog(fx) - alog(a(0)) ) / vx
#IF verbose THEN print, 'FitLogDecay.pro: Function estimates: ', A

#now look at the weights
#remove where weights are 0
    keep = weights != 0
    t_keep = t[keep]
    y_keep = y[keep]
    weights_keep = weights[keep]
#high weights mean important points, i.e. low variability = sigma
    sigma = 1/weights_keep

#;nov 07, new initial estimates, follow this for python
    A, K, offset = fit_exp_linear_autoOffset(t_keep, y_keep) #
    fittedout_lin = A * np.exp(K*t) + offset
    logging.getLogger("VIEW").debug(
        f'ViewCalcData:fitlogdecay: linear fit parameters A = {A}, K = {K} . Offest={offset}')
    # initial guess done without weights, but excluding excluded values.
# A and K are the estimates

# give feedback on linear fit
#    opt_parms = (A, K, offset)
#    fittedout = A * np.exp(K*t) + offset
#    if showresults: show_fitlogdecay(y, fittedout, t, weights, opt_parms)

#itMax = 1000
#yfit = CURVEFIT(X, Y, weights, A, SIGMA, FUNCTION_NAME='gfunct', ITmax=itMax, ITER=iter, status=status)
    import scipy.optimize as spo
    # Linear fit is better than exponential fit here, as y_keep is flat
    if A == 0:
        opt_parms = (A, K, offset)
    else:
        try:
            # default value of maxfev is 800, not always sufficient
            opt_parms, parm_cov = spo.curve_fit(
                model_func, t_keep, y_keep,  p0=(A,K,offset), sigma = sigma, maxfev=3200,
                # setting bounds on A, to make sure that the fitted curve does not rise/fall faster than y_keep does.
                # it helps to reduce the exponential divergence between fittedout and y for points with 0 weight.
                bounds=[
                    (min(y_keep.min() - offset, A), -np.inf, -np.inf),  # lower bounds
                    (max(y_keep.max() - offset, A), np.inf, np.inf)  # upper bounds
                ]
            )
            logging.getLogger("VIEW").debug("Fitting log function converged! Using Log fit")
        except RuntimeError as rte:
            opt_parms = (A, K, offset)
            logging.getLogger("VIEW").debug("Fitting log function did NOT converge. Using linear parameters.")




#;Compute the parameters.
#;IF ((iter gt itmax) OR (iter eq 1)) THEN noConverge=1 ELSE noConverge=0
#IF (status ge 1) THEN noConverge=1 ELSE noConverge=0
#IF verbose THEN PRINT, 'FitLogDecay.pro: Function parameters (check FitLogDecay.pro for function): ', A	;Print the parameters returned in A.
#IF noConverge THEN print, 'FitLogDecay.pro: Failed to converge, setting parameters to 0; check FitLogDecay.pro'
#
#;copy the parameters into the p structure, so it is available for later retrieval
#p1.BleachPar = A


#;IDL prints:
#;Function parameters:       9.91120    -0.100883      2.07773
#;Thus, the function that best fits the data is:
#;F(x) = 9.91120 * exp(-0.100883x) +  2.07773

    A, K, C = opt_parms
    fittedout = A * np.exp(K*t) + C
### show result
    if showresults: show_fitlogdecay(y, fittedout, fittedout_lin, t, weights, opt_parms, measurement_label)
# end function fitlogdecay
    return (fittedout, opt_parms) #, A, K, C


def show_fitlogdecay(lineIn, fittedout_blue, fittedout_green, t, weights, opt_parms, measurement_label):

    A, K, C = opt_parms #extract function parameters
    fig = plt.figure()
    fig.canvas.set_window_title(measurement_label)
    ax1 = fig.add_subplot(2,1,1) # (2,1,1)
    ax2 = fig.add_subplot(2,1,2)
    ax1.set_title('Fitlogdecay: ' + measurement_label)
    ax2.set_title('weights')
    ax2.plot(t, weights, 'k--',
          label='weights')
    ax1.plot(t, lineIn, 'ro')
    ax1.plot(t, fittedout_blue, 'b-',
          label='$y = %0.2f e^{%0.2f t} + %0.2f$' % (A, K, C))
    ax1.plot(t, fittedout_green, 'g-',
          label='Alternate lin log-fit')
#    ax1.legend(bbox_to_anchor=(1, 1), fancybox=True, shadow=True)
    ax1.legend(fancybox=True, shadow=True)
    plt.show(block=False)
    # move window to the front - no idea if this works


def get_bleach_weights(
        flags, p1_metadata: pd.Series, movie_size: tuple, exclude_stimulus: bool):
    """
    Calculates weights to be used for bleach correction. Weights for the frames 0 to `end_background` are set
    flags["LELog_InitialFactor"]. If `flags["LE_BleachStartFrame"]` is > 0, frames 0 to this value is set to 0. If
    `exclude_stimulus` is True, then the weights for "stimulus frames" are set to 0, where "stimulus frames" are defined
    to be from `end_background` to `flags["LELog_ExcludeSeconds"]` seconds after the onset of first stimulus
    :param FlagsManager flags:
    :param pandas.Series p1_metadata: experimental metadata
    :param tuple movie_size: size of raw data, format XYT
    :param bool exclude_stimulus: when true, "stimulus frames" as defined above are set to have 0 weight
    :return: an array of weights
    :rtype: numpy.ndarray
    """

    n_frames = movie_size[2]
    end_background = p1_metadata.background_frames[1]
    sampling_frequency = p1_metadata.frequency

    # can't exclude stimulus if stimulus info is not provided
    exclude_stimulus = exclude_stimulus and p1_metadata.pulsed_stimuli_handler is not None

    weights = np.ones(n_frames)

    #  change weight to initial portion (with safety net)
    weights[0: end_background] = flags["LELog_InitialFactor"]

    #  exclude first frames
    bleach_start_frame = flags["LE_BleachStartFrame"]
    if bleach_start_frame > 0:
        weights[0:bleach_start_frame] = 0

    #  exclude response interval
    if exclude_stimulus:
        startweight = end_background
        endweight = min([round(end_background + flags["LELog_ExcludeSeconds"] * sampling_frequency), n_frames - 1])
        weights[startweight:endweight] = 0

    return weights