decide whether to use ipmag.smooth() within rockmag.py vs new functions

[Swanson-Hysell]

Within ipmag.py there is a function smooth(). We could use this function rather than a new moving average within rockmag.py. The default for smoothing is bartlett which seems sensible rather than a flat window. Note that there is preprocessing of the data within ipmag.curie() to make it so that the input array is evenly spaced.

as a note about code functionality: if the window parameter is 'flat', a simple moving average is created by generating an array of ones using np.ones(window_len, 'd'), where window_len is the length of the window and 'd' specifies the dtype as double precision floating point.
For other window types ('hanning', 'hamming', 'bartlett', 'blackman'), the code dynamically constructs a string that corresponds to the command for creating the specified window type in NumPy. For example, if window is 'hanning', the string np.hanning(window_len) is constructed and then evaluated using eval. This results in a call to np.hanning(window_len), which generates a Hanning window of length window_len. After constructing that window, it is then processed through np.convolve()

def smooth(x, window_len, window='bartlett'):
    """
    Smooth the data using a sliding window with requested size - meant to be
    used with the ipmag function curie().
    This method is based on the convolution of a scaled window with the signal.
    The signal is prepared by padding the beginning and the end of the signal
    with average of the first (last) ten values of the signal, to avoid jumps
    at the beginning/end. Output is an array of the smoothed signal.

    Required Parameters
    ----------
    x : the input signal, equally spaced!
    window_len : the dimension of the smoothing window

    Optional Parameters (defaults are used if not specified)
    ----------
    window : type of window from numpy library ['flat','hanning','hamming','bartlett','blackman']
        (default is Bartlett)
        -flat window will produce a moving average smoothing.
        -Bartlett window is very similar to triangular window,
            but always ends with zeros at points 1 and n.
        -hanning,hamming,blackman are used for smoothing the Fourier transform
    """
    if x.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")

    if x.size < window_len:
        raise ValueError("Input vector needs to be bigger than window size.")

    if window_len < 3:
        return x

    # numpy available windows
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError(
            "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")

    # padding the beginning and the end of the signal with an average value to
    # evoid edge effect
    start = [np.average(x[0:10])] * window_len
    end = [np.average(x[-10:])] * window_len
    s = start + list(x) + end

    # s=numpy.r_[2*x[0]-x[window_len:1:-1],x,2*x[-1]-x[-1:-window_len:-1]]
    if window == 'flat':  # moving average
        w = ones(window_len, 'd')
    else:
        w = eval('np.' + window + '(window_len)')
    y = np.convolve(w/w.sum(), s, mode='same')
    return np.array(y[window_len:-window_len])

The code that makes the input be evenly spaced can be seen in ipmag.curie() in the block of code before M_smooth = smooth(M, window_len).

def curie(path_to_file='.', file_name='', magic=False,
          window_length=3, save=False, save_folder='.', fmt='svg', t_begin="", t_end=""):
    """
    Plots and interprets curie temperature data.
    The 1st derivative is calculated from smoothed M-T curve (convolution
    with triangular window with width= <-w> degrees)
    The 2nd derivative is calculated from smoothed 1st derivative curve
    (using the same sliding window width)
    The estimated curie temperation is the maximum of the 2nd derivative.
    Temperature steps should be in multiples of 1.0 degrees.

    Parameters:
        file_name : name of file to be opened
        path_to_file : path to directory that contains file (default is current directory, '.')
        window_length : dimension of smoothing window (input to smooth() function)
        save : boolean argument to save plots (default is False)
        save_folder : relative directory where plots will be saved (default is current directory, '.')
        fmt : format of saved figures (default is svg)
        t_begin: start of truncated window for search (default is beginning of data)
        t_end: end of truncated window for search (default is end of data)
        magic : True if MagIC formatted measurements.txt file
        
    Returns:
        A plot is shown and saved if save=True.
    """
    plot = 0
    window_len = window_length

    # read data from file
    complete_path = os.path.join(path_to_file, file_name)
    if magic:
        data_df = pd.read_csv(complete_path, sep='\t', header=1)
        T = data_df['meas_temp'].values-273
        magn_key = cb.get_intensity_col(data_df)
        M = data_df[magn_key].values
    else:
        Data = np.loadtxt(complete_path, dtype=np.float)
        T = Data.transpose()[0]
        M = Data.transpose()[1]
    T = list(T)
    M = list(M)
    # cut the data if -t is one of the flags
    if t_begin != "":
        while T[0] < t_begin:
            M.pop(0)
            T.pop(0)
        while T[-1] > t_end:
            M.pop(-1)
            T.pop(-1)

    # prepare the signal:
    # from M(T) array with unequal deltaT
    # to M(T) array with deltaT=(1 degree).
    # if delataT is larger, then points are added using linear fit between
    # consecutive data points.
    # exit if deltaT is not integer
    i = 0
    while i < (len(T) - 1):
        if (T[i + 1] - T[i]) % 1 > 0.001:
            print("delta T should be integer, this program will not work!")
            print("temperature range:", T[i], T[i + 1])
            sys.exit()
        if (T[i + 1] - T[i]) == 0.:
            M[i] = np.average([M[i], M[i + 1]])
            M.pop(i + 1)
            T.pop(i + 1)
        elif (T[i + 1] - T[i]) < 0.:
            M.pop(i + 1)
            T.pop(i + 1)
            print("check data in T=%.0f ,M[T] is ignored" % (T[i]))
        elif (T[i + 1] - T[i]) > 1.:
            slope, b = np.polyfit([T[i], T[i + 1]], [M[i], M[i + 1]], 1)
            for j in range(int(T[i + 1]) - int(T[i]) - 1):
                M.insert(i + 1, slope * (T[i] + 1.) + b)
                T.insert(i + 1, (T[i] + 1.))
                i = i + 1
        i = i + 1

    # calculate the smoothed signal
    M = np.array(M, 'f')
    T = np.array(T, 'f')
    M_smooth = []
    M_smooth = smooth(M, window_len)

    # plot the original data and the smooth data
    PLT = {'M_T': 1, 'der1': 2, 'der2': 3, 'Curie': 4}
    plt.figure(num=PLT['M_T'], figsize=(5, 5))
    string = 'M-T (sliding window=%i)' % int(window_len)
    pmagplotlib.plot_xy(PLT['M_T'], T, M_smooth, sym='-')
    pmagplotlib.plot_xy(PLT['M_T'], T, M, sym='--',
                        xlab='Temperature C', ylab='Magnetization', title=string)

    # calculate first derivative
    d1, T_d1 = [], []
    for i in range(len(M_smooth) - 1):
        Dy = M_smooth[i - 1] - M_smooth[i + 1]
        Dx = T[i - 1] - T[i + 1]
        d1.append(Dy/Dx)
    T_d1 = T[1:len(T - 1)]
    d1 = np.array(d1, 'f')
    d1_smooth = smooth(d1, window_len)

    # plot the first derivative
    plt.figure(num=PLT['der1'], figsize=(5, 5))
    string = '1st derivative (sliding window=%i)' % int(window_len)
    pmagplotlib.plot_xy(PLT['der1'], T_d1, d1_smooth,
                        sym='-', xlab='Temperature C', title=string)
    pmagplotlib.plot_xy(PLT['der1'], T_d1, d1, sym='b--')

    # calculate second derivative
    d2, T_d2 = [], []
    for i in range(len(d1_smooth) - 1):
        Dy = d1_smooth[i - 1] - d1_smooth[i + 1]
        Dx = T[i - 1] - T[i + 1]
        # print Dy/Dx
        d2.append(Dy/Dx)
    T_d2 = T[2:len(T - 2)]
    d2 = np.array(d2, 'f')
    d2_smooth = smooth(d2, window_len)

    # plot the second derivative
    plt.figure(num=PLT['der2'], figsize=(5, 5))
    string = '2nd derivative (sliding window=%i)' % int(window_len)
    pmagplotlib.plot_xy(PLT['der2'], T_d2, d2, sym='-',
                        xlab='Temperature C', title=string)
    d2 = list(d2)
    print('second derivative maximum is at T=%i' %
          int(T_d2[d2.index(max(d2))]))

    # calculate Curie temperature for different width of sliding windows
    curie, curie_1 = [], []
    wn = list(range(5, 50, 1))
    for win in wn:
        # calculate the smoothed signal
        M_smooth = []
        M_smooth = smooth(M, win)
        # calculate first derivative
        d1, T_d1 = [], []
        for i in range(len(M_smooth) - 1):
            Dy = M_smooth[i - 1] - M_smooth[i + 1]
            Dx = T[i - 1] - T[i + 1]
            d1.append(Dy/Dx)
        T_d1 = T[1:len(T - 1)]
        d1 = np.array(d1, 'f')
        d1_smooth = smooth(d1, win)
        # calculate second derivative
        d2, T_d2 = [], []
        for i in range(len(d1_smooth) - 1):
            Dy = d1_smooth[i - 1] - d1_smooth[i + 1]
            Dx = T[i - 1] - T[i + 1]
            d2.append(Dy/Dx)
        T_d2 = T[2:len(T - 2)]
        d2 = np.array(d2, 'f')
        d2_smooth = smooth(d2, win)
        d2 = list(d2)
        d2_smooth = list(d2_smooth)
        curie.append(T_d2[d2.index(max(d2))])
        curie_1.append(T_d2[d2_smooth.index(max(d2_smooth))])

    # plot Curie temp for different sliding window length
    plt.figure(num=PLT['Curie'], figsize=(5, 5))
    pmagplotlib.plot_xy(PLT['Curie'], wn, curie, sym='.',
                        xlab='Sliding Window Width (degrees)', ylab='Curie Temp', title='Curie Statistics')
    files = {}
    for key in list(PLT.keys()):
        files[key] = str(key) + '.' + fmt
    if save == True:
        for key in list(PLT.keys()):
            try:
                plt.figure(num=PLT[key])
                plt.savefig(save_folder + '/' + files[key].replace('/', '-'))
            except:
                print('could not save: ', PLT[key], files[key])
                print("output file format not supported ")
    plt.show()

[Swanson-Hysell]

One thing to be aware of with this code is that the preprocessing function within ipmag.curie() requires that deltaT be an integer. One approach could be to develop another resampling algorithm that relaxes this requirement while outputting a resampled evenly spaced array that can then be fed into ipmag.smooth()

[josh-feinberg]

I like your last suggestion of creating a modified algorithm that relaxes the integer requirement for deltaT, while still outputting an evenly spaced array.