treetime/node_interpolator.py
def convolve(cls, node_interp, branch_interp, max_or_integral='integral',
                 n_grid_points = ttconf.NODE_GRID_SIZE, n_integral=ttconf.N_INTEGRAL,
                 inverse_time=True, rel_tol=0.05, yc=10):

        '''
        calculate H(t) = \int_tau f(t-tau)g(tau) if inverse_time=True
                  H(t) = \int_tau f(t+tau)g(tau) if inverse_time=False

        This function determines the time points of the grid of the result to
        ensure an accurate approximation.
        '''

        if max_or_integral not in ['max', 'integral']:
            raise Exception("Max_or_integral expected to be 'max' or 'integral', got "
                            + str(max_or_integral)  + " instead.")

        def conv_in_point(time_point):

            if max_or_integral == 'integral': # compute integral of the convolution
                return _evaluate_convolution(time_point, node_interp, branch_interp,
                                               n_integral=n_integral, return_log=True,
                                               inverse_time = inverse_time)

            else: # compute max of the convolution
                return _max_of_integrand(time_point, node_interp, branch_interp,
                                               return_log=True, inverse_time = inverse_time)

        # estimate peak and width
        joint_fwhm  = (node_interp.fwhm + branch_interp.fwhm)
        min_fwhm  = min(node_interp.fwhm, branch_interp.fwhm)
        # determine support of the resulting convolution
        # in order to be positive, the flipped support of f, shifted by t and g need to overlap
        if inverse_time:
            new_peak_pos = node_interp.peak_pos + branch_interp.peak_pos
            tmin = node_interp.xmin+branch_interp.xmin
            tmax = node_interp.xmax+branch_interp.xmax
        else:
            new_peak_pos = node_interp.peak_pos - branch_interp.peak_pos
            tmin = node_interp.xmin - branch_interp.xmax
            tmax = node_interp.xmax - branch_interp.xmin

        # make initial node grid consisting of linearly spaced points around
        # the center and quadratically spaced points at either end
        n = n_grid_points/3
        center_width = 3*joint_fwhm
        grid_center = new_peak_pos + np.linspace(-1, 1, n)*center_width

        # add the right and left grid if it is needed
        right_range = (tmax - grid_center[-1])
        if right_range>4*center_width:
            grid_right = grid_center[-1] + right_range*(np.linspace(0, 1, n)**2.0)
        elif right_range>0: # use linear grid the right_range is comparable to center_width
            grid_right = grid_center[-1] + right_range*np.linspace(0,1, int(min(n,1+0.5*n*right_range/center_width)))
        else:
            grid_right =[]

        left_range = grid_center[0]-tmin
        if left_range>4*center_width:
            grid_left = tmin + left_range*(np.linspace(0, 1, n)**2.0)
        elif left_range>0:
            grid_left = tmin + left_range*np.linspace(0,1, int(min(n,1+0.5*n*left_range/center_width)))
        else:
            grid_left =[]


        if tmin>-1:
            grid_zero_left = tmin + (tmax-tmin)*np.linspace(0,0.01,11)**2
        else:
            grid_zero_left = [tmin]
        if tmax<1:
            grid_zero_right = tmax - (tmax-tmin)*np.linspace(0,0.01,11)**2
        else:
            grid_zero_right = [tmax]

        # make grid and calculate convolution
        t_grid_0 = np.unique(np.concatenate([grid_zero_left, grid_left[:-1], grid_center, grid_right[1:], grid_zero_right]))
        t_grid_0 = t_grid_0[(t_grid_0 > tmin-ttconf.TINY_NUMBER) & (t_grid_0 < tmax+ttconf.TINY_NUMBER)]

        # res0 - the values of the convolution (integral or max)
        # t_0  - the value, at which the res0 achieves maximum
        #        (when determining the maximum of the integrand, otherwise meaningless)
        res_0, t_0 = np.array([conv_in_point(t_val) for t_val in t_grid_0]).T

        # refine grid as necessary and add new points
        # calculate interpolation error at all internal points [2:-2] bc end points are sometime off scale
        interp_error = np.abs(res_0[3:-1]+res_0[1:-3]-2*res_0[2:-2])
        # determine the number of extra points needed, criterion depends on distance from peak dy
        dy = (res_0[2:-2]-res_0.min())
        dx = np.diff(t_grid_0)
        refine_factor = np.minimum(np.minimum(np.array(np.floor(np.sqrt(interp_error/(rel_tol*(1+(dy/yc)**4)))), dtype=int),
                                   np.array(100*(dx[1:-2]+dx[2:-1])/min_fwhm, dtype=int)), 10)

        insert_point_idx = np.zeros(interp_error.shape[0]+1, dtype=int)
        insert_point_idx[1:] = refine_factor
        insert_point_idx[:-1] += refine_factor
        # add additional points if there are any to add
        if np.sum(insert_point_idx):
            add_x = np.concatenate([np.linspace(t1,t2,n+2)[1:-1] for t1,t2,n in
                               zip(t_grid_0[1:-2], t_grid_0[2:-1], insert_point_idx) if n>0])
            # calculate convolution at these points
            add_y, add_t = np.array([conv_in_point(t_val) for t_val in add_x]).T

            t_grid_0 = np.concatenate((t_grid_0, add_x))
            res_0 = np.concatenate ((res_0, add_y))
            t_0 = np.concatenate ((t_0, add_t))

        # instantiate the new interpolation object and return
        res_y = cls(t_grid_0, res_0, is_log=True, kind='linear')

        # the interpolation object, which is used to store the value of the
        # grid, which maximizes the convolution (for 'max' option),
        # or flat -1 distribution (for 'integral' option)
        # this grid is the optimal branch length
        res_t = Distribution(t_grid_0, t_0, is_log=True,
                             min_width=node_interp.min_width, kind='linear')

        return res_y, res_t
def convolve(cls, node_interp, branch_interp, max_or_integral='integral',
                 n_grid_points = ttconf.NODE_GRID_SIZE, n_integral=ttconf.N_INTEGRAL,
                 inverse_time=True, rel_tol=0.05, yc=10):

        '''
        calculate H(t) = \int_tau f(t-tau)g(tau) if inverse_time=True
                  H(t) = \int_tau f(t+tau)g(tau) if inverse_time=False

        This function determines the time points of the grid of the result to
        ensure an accurate approximation.
        '''

        if max_or_integral not in ['max', 'integral']:
            raise Exception("Max_or_integral expected to be 'max' or 'integral', got "
                            + str(max_or_integral)  + " instead.")

        def conv_in_point(time_point):

            if max_or_integral == 'integral': # compute integral of the convolution
                return _evaluate_convolution(time_point, node_interp, branch_interp,
                                               n_integral=n_integral, return_log=True,
                                               inverse_time = inverse_time)

            else: # compute max of the convolution
                return _max_of_integrand(time_point, node_interp, branch_interp,
                                               return_log=True, inverse_time = inverse_time)

        # estimate peak and width
        joint_fwhm  = (node_interp.fwhm + branch_interp.fwhm)
        min_fwhm  = min(node_interp.fwhm, branch_interp.fwhm)
        # determine support of the resulting convolution
        # in order to be positive, the flipped support of f, shifted by t and g need to overlap
        if inverse_time:
            new_peak_pos = node_interp.peak_pos + branch_interp.peak_pos
            tmin = node_interp.xmin+branch_interp.xmin
            tmax = node_interp.xmax+branch_interp.xmax
        else:
            new_peak_pos = node_interp.peak_pos - branch_interp.peak_pos
            tmin = node_interp.xmin - branch_interp.xmax
            tmax = node_interp.xmax - branch_interp.xmin

        # make initial node grid consisting of linearly spaced points around
        # the center and quadratically spaced points at either end
        n = n_grid_points/3
        center_width = 3*joint_fwhm
        grid_center = new_peak_pos + np.linspace(-1, 1, n)*center_width

        # add the right and left grid if it is needed
        right_range = (tmax - grid_center[-1])
        if right_range>4*center_width:
            grid_right = grid_center[-1] + right_range*(np.linspace(0, 1, n)**2.0)
        elif right_range>0: # use linear grid the right_range is comparable to center_width
            grid_right = grid_center[-1] + right_range*np.linspace(0,1, int(min(n,1+0.5*n*right_range/center_width)))
        else:
            grid_right =[]

        left_range = grid_center[0]-tmin
        if left_range>4*center_width:
            grid_left = tmin + left_range*(np.linspace(0, 1, n)**2.0)
        elif left_range>0:
            grid_left = tmin + left_range*np.linspace(0,1, int(min(n,1+0.5*n*left_range/center_width)))
        else:
            grid_left =[]


        if tmin>-1:
            grid_zero_left = tmin + (tmax-tmin)*np.linspace(0,0.01,11)**2
        else:
            grid_zero_left = [tmin]
        if tmax<1:
            grid_zero_right = tmax - (tmax-tmin)*np.linspace(0,0.01,11)**2
        else:
            grid_zero_right = [tmax]

        # make grid and calculate convolution
        t_grid_0 = np.unique(np.concatenate([grid_zero_left, grid_left[:-1], grid_center, grid_right[1:], grid_zero_right]))
        t_grid_0 = t_grid_0[(t_grid_0 > tmin-ttconf.TINY_NUMBER) & (t_grid_0 < tmax+ttconf.TINY_NUMBER)]

        # res0 - the values of the convolution (integral or max)
        # t_0  - the value, at which the res0 achieves maximum
        #        (when determining the maximum of the integrand, otherwise meaningless)
        res_0, t_0 = np.array([conv_in_point(t_val) for t_val in t_grid_0]).T

        # refine grid as necessary and add new points
        # calculate interpolation error at all internal points [2:-2] bc end points are sometime off scale
        interp_error = np.abs(res_0[3:-1]+res_0[1:-3]-2*res_0[2:-2])
        # determine the number of extra points needed, criterion depends on distance from peak dy
        dy = (res_0[2:-2]-res_0.min())
        dx = np.diff(t_grid_0)
        refine_factor = np.minimum(np.minimum(np.array(np.floor(np.sqrt(interp_error/(rel_tol*(1+(dy/yc)**4)))), dtype=int),
                                   np.array(100*(dx[1:-2]+dx[2:-1])/min_fwhm, dtype=int)), 10)

        insert_point_idx = np.zeros(interp_error.shape[0]+1, dtype=int)
        insert_point_idx[1:] = refine_factor
        insert_point_idx[:-1] += refine_factor
        # add additional points if there are any to add
        if np.sum(insert_point_idx):
            add_x = np.concatenate([np.linspace(t1,t2,n+2)[1:-1] for t1,t2,n in
                               zip(t_grid_0[1:-2], t_grid_0[2:-1], insert_point_idx) if n>0])
            # calculate convolution at these points
            add_y, add_t = np.array([conv_in_point(t_val) for t_val in add_x]).T

            t_grid_0 = np.concatenate((t_grid_0, add_x))
            res_0 = np.concatenate ((res_0, add_y))
            t_0 = np.concatenate ((t_0, add_t))

        # instantiate the new interpolation object and return
        res_y = cls(t_grid_0, res_0, is_log=True, kind='linear')

        # the interpolation object, which is used to store the value of the
        # grid, which maximizes the convolution (for 'max' option),
        # or flat -1 distribution (for 'integral' option)
        # this grid is the optimal branch length
        res_t = Distribution(t_grid_0, t_0, is_log=True,
                             min_width=node_interp.min_width, kind='linear')

        return res_y, res_t
def convolve ( cls , node_interp , branch_interp , max_or_integral = 'integral' , n_grid_points = ttconf . NODE_GRID_SIZE , n_integral = ttconf . N_INTEGRAL , inverse_time = True , rel_tol = 0.05 , yc = 10 ) : if max_or_integral not in [ 'max' , 'integral' ] : raise Exception ( "Max_or_integral expected to be 'max' or 'integral', got " + str ( max_or_integral ) + " instead." ) def conv_in_point ( time_point ) : if max_or_integral == 'integral' : # compute integral of the convolution return _evaluate_convolution ( time_point , node_interp , branch_interp , n_integral = n_integral , return_log = True , inverse_time = inverse_time ) else : # compute max of the convolution return _max_of_integrand ( time_point , node_interp , branch_interp , return_log = True , inverse_time = inverse_time ) # estimate peak and width joint_fwhm = ( node_interp . fwhm + branch_interp . fwhm ) min_fwhm = min ( node_interp . fwhm , branch_interp . fwhm ) # determine support of the resulting convolution # in order to be positive, the flipped support of f, shifted by t and g need to overlap if inverse_time : new_peak_pos = node_interp . peak_pos + branch_interp . peak_pos tmin = node_interp . xmin + branch_interp . xmin tmax = node_interp . xmax + branch_interp . xmax else : new_peak_pos = node_interp . peak_pos - branch_interp . peak_pos tmin = node_interp . xmin - branch_interp . xmax tmax = node_interp . xmax - branch_interp . xmin # make initial node grid consisting of linearly spaced points around # the center and quadratically spaced points at either end n = n_grid_points / 3 center_width = 3 * joint_fwhm grid_center = new_peak_pos + np . linspace ( - 1 , 1 , n ) * center_width # add the right and left grid if it is needed right_range = ( tmax - grid_center [ - 1 ] ) if right_range > 4 * center_width : grid_right = grid_center [ - 1 ] + right_range * ( np . linspace ( 0 , 1 , n ) ** 2.0 ) elif right_range > 0 : # use linear grid the right_range is comparable to center_width grid_right = grid_center [ - 1 ] + right_range * np . linspace ( 0 , 1 , int ( min ( n , 1 + 0.5 * n * right_range / center_width ) ) ) else : grid_right = [ ] left_range = grid_center [ 0 ] - tmin if left_range > 4 * center_width : grid_left = tmin + left_range * ( np . linspace ( 0 , 1 , n ) ** 2.0 ) elif left_range > 0 : grid_left = tmin + left_range * np . linspace ( 0 , 1 , int ( min ( n , 1 + 0.5 * n * left_range / center_width ) ) ) else : grid_left = [ ] if tmin > - 1 : grid_zero_left = tmin + ( tmax - tmin ) * np . linspace ( 0 , 0.01 , 11 ) ** 2 else : grid_zero_left = [ tmin ] if tmax < 1 : grid_zero_right = tmax - ( tmax - tmin ) * np . linspace ( 0 , 0.01 , 11 ) ** 2 else : grid_zero_right = [ tmax ] # make grid and calculate convolution t_grid_0 = np . unique ( np . concatenate ( [ grid_zero_left , grid_left [ : - 1 ] , grid_center , grid_right [ 1 : ] , grid_zero_right ] ) ) t_grid_0 = t_grid_0 [ ( t_grid_0 > tmin - ttconf . TINY_NUMBER ) & ( t_grid_0 < tmax + ttconf . TINY_NUMBER ) ] # res0 - the values of the convolution (integral or max) # t_0  - the value, at which the res0 achieves maximum #        (when determining the maximum of the integrand, otherwise meaningless) res_0 , t_0 = np . array ( [ conv_in_point ( t_val ) for t_val in t_grid_0 ] ) . T # refine grid as necessary and add new points # calculate interpolation error at all internal points [2:-2] bc end points are sometime off scale interp_error = np . abs ( res_0 [ 3 : - 1 ] + res_0 [ 1 : - 3 ] - 2 * res_0 [ 2 : - 2 ] ) # determine the number of extra points needed, criterion depends on distance from peak dy dy = ( res_0 [ 2 : - 2 ] - res_0 . min ( ) ) dx = np . diff ( t_grid_0 ) refine_factor = np . minimum ( np . minimum ( np . array ( np . floor ( np . sqrt ( interp_error / ( rel_tol * ( 1 + ( dy / yc ) ** 4 ) ) ) ) , dtype = int ) , np . array ( 100 * ( dx [ 1 : - 2 ] + dx [ 2 : - 1 ] ) / min_fwhm , dtype = int ) ) , 10 ) insert_point_idx = np . zeros ( interp_error . shape [ 0 ] + 1 , dtype = int ) insert_point_idx [ 1 : ] = refine_factor insert_point_idx [ : - 1 ] += refine_factor # add additional points if there are any to add if np . sum ( insert_point_idx ) : add_x = np . concatenate ( [ np . linspace ( t1 , t2 , n + 2 ) [ 1 : - 1 ] for t1 , t2 , n in zip ( t_grid_0 [ 1 : - 2 ] , t_grid_0 [ 2 : - 1 ] , insert_point_idx ) if n > 0 ] ) # calculate convolution at these points add_y , add_t = np . array ( [ conv_in_point ( t_val ) for t_val in add_x ] ) . T t_grid_0 = np . concatenate ( ( t_grid_0 , add_x ) ) res_0 = np . concatenate ( ( res_0 , add_y ) ) t_0 = np . concatenate ( ( t_0 , add_t ) ) # instantiate the new interpolation object and return res_y = cls ( t_grid_0 , res_0 , is_log = True , kind = 'linear' ) # the interpolation object, which is used to store the value of the # grid, which maximizes the convolution (for 'max' option), # or flat -1 distribution (for 'integral' option) # this grid is the optimal branch length res_t = Distribution ( t_grid_0 , t_0 , is_log = True , min_width = node_interp . min_width , kind = 'linear' ) return res_y , res_t
calculate H(t) = \int_tau f(t-tau)g(tau) if inverse_time=True
                  H(t) = \int_tau f(t+tau)g(tau) if inverse_time=False

        This function determines the time points of the grid of the result to
        ensure an accurate approximation.
calculate H ( t ) = \ int_tau f ( t - tau ) g ( tau ) if inverse_time = True H ( t ) = \ int_tau f ( t + tau ) g ( tau ) if inverse_time = False