<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def haarpsi_similarity_map(img1, img2, axis, c, a):
 r"""Local similarity map for directional features along an axis. Parameters img1, img2 : array-like The images to compare. They must have equal shape. axis : {0, 1} Direction in which to look for edge similarities. c : positive float Constant determining the score of maximally dissimilar values. Smaller constant means higher penalty for dissimilarity. See Notes for details. a : positive float Parameter in the logistic function. Larger value leads to a steeper curve, thus lowering the threshold for an input to be mapped to an output close to 1. See Notes for details. Returns ------- local_sim : `numpy.ndarray` Pointwise similarity of directional edge features of ``img1`` and ``img2``, measured using two Haar wavelet detail levels. Notes ----- For input images :math:`f_1, f_2` this function is defined as .. math:: \mathrm{HS}_{f_1, f_2}^{(k)}(x) = l_a \left( \frac{1}{2} \sum_{j=1}^2 S\left(\left|g_j^{(k)} \ast f_1 \right|(x), \left|g_j^{(k)} \ast f_2 \right|(x), c\right) \right), see `[Rei+2016] <https://arxiv.org/abs/1607.06140>`_ equation (10). Here, the superscript :math:`(k)` refers to the axis (0 or 1) in which edge features are compared, :math:`l_a` is the logistic function :math:`l_a(x) = (1 + \mathrm{e}^{-a x})^{-1}`, and :math:`S` is the pointwise similarity score .. math:: S(x, y, c) = \frac{2xy + c^2}{x^2 + y^2 + c^2}, Hence, :math:`c` is the :math:`y`-value at which the score drops to :math:`1 / 2` for :math:`x = 0`. In other words, the smaller :math:`c` is chosen, the more dissimilarity is penalized. The filters :math:`g_j^{(k)}` are high-pass Haar wavelet filters in the axis :math:`k` and low-pass Haar wavelet filters in the other axes. The index :math:`j` refers to the scaling level of the wavelet. In code, these filters can be computed as :: f_lo_level1 = [np.sqrt(2), np.sqrt(2)] # low-pass Haar filter f_hi_level1 = [-np.sqrt(2), np.sqrt(2)] # high-pass Haar filter f_lo_level2 = np.repeat(f_lo_level1, 2) f_hi_level2 = np.repeat(f_hi_level1, 2) f_lo_level3 = np.repeat(f_lo_level2, 2) f_hi_level3 = np.repeat(f_hi_level2, 2) The logistic function :math:`l_a` transforms values in :math:`[0, \infty)` to :math:`[1/2, 1)`, where the parameter :math:`a` determines how fast the curve attains values close to 1. Larger :math:`a` means that smaller :math:`x` will yield a value :math:`l_a(x)` close to 1 (and thus result in a higher score). In other words, the larger :math:`a`, the more forgiving the similarity measure. References [Rei+2016] Reisenhofer, R, Bosse, S, Kutyniok, G, and Wiegand, T. *A Haar Wavelet-Based Perceptual Similarity Index for Image Quality Assessment*. arXiv:1607.06140 [cs], Jul. 2016. """

    # TODO: generalize for nD
    import scipy.special
    impl = 'pyfftw' if PYFFTW_AVAILABLE else 'numpy'

    # Haar wavelet filters for levels 1 and 2
    dec_lo_lvl1 = np.array([np.sqrt(2), np.sqrt(2)])
    dec_lo_lvl2 = np.repeat(dec_lo_lvl1, 2)
    dec_hi_lvl1 = np.array([-np.sqrt(2), np.sqrt(2)])
    dec_hi_lvl2 = np.repeat(dec_hi_lvl1, 2)

    if axis == 0:
        # High-pass in axis 0, low-pass in axis 1
        fh_lvl1 = dec_hi_lvl1
        fv_lvl1 = dec_lo_lvl1
        fh_lvl2 = dec_hi_lvl2
        fv_lvl2 = dec_lo_lvl2
    elif axis == 1:
        # Low-pass in axis 0, high-pass in axis 1
        fh_lvl1 = dec_lo_lvl1
        fv_lvl1 = dec_hi_lvl1
        fh_lvl2 = dec_lo_lvl2
        fv_lvl2 = dec_hi_lvl2
    else:
        raise ValueError('`axis` out of the valid range 0 -> 1')

    # Filter images with level 1 and 2 filters
    img1_lvl1 = filter_image_sep2d(img1, fh_lvl1, fv_lvl1, impl=impl)
    img1_lvl2 = filter_image_sep2d(img1, fh_lvl2, fv_lvl2, impl=impl)

    img2_lvl1 = filter_image_sep2d(img2, fh_lvl1, fv_lvl1, impl=impl)
    img2_lvl2 = filter_image_sep2d(img2, fh_lvl2, fv_lvl2, impl=impl)

    c = float(c)

    def S(x, y):
        """Return ``(2 * x * y + c ** 2) / (x ** 2 + y ** 2 + c ** 2)``."""
        num = 2 * x
        num *= y
        num += c ** 2
        denom = x ** 2
        denom += y ** 2
        denom += c ** 2
        frac = num
        frac /= denom
        return frac

    # Compute similarity scores for both levels
    np.abs(img1_lvl1, out=img1_lvl1)
    np.abs(img2_lvl1, out=img2_lvl1)
    np.abs(img1_lvl2, out=img1_lvl2)
    np.abs(img2_lvl2, out=img2_lvl2)

    sim_lvl1 = S(img1_lvl1, img2_lvl1)
    sim_lvl2 = S(img1_lvl2, img2_lvl2)

    # Return logistic of the mean value
    sim = sim_lvl1
    sim += sim_lvl2
    sim /= 2
    sim *= a
    return scipy.special.expit(sim)