123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240 |
- 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
|