pyart中关于线性规划部分代码的疑问

# 第一部分
# parse the field parameters
if refl_field is None:
    refl_field = get_field_name("reflectivity")
if ncp_field is None:
    ncp_field = get_field_name("normalized_coherent_power")
if rhv_field is None:
    rhv_field = get_field_name("cross_correlation_ratio")
if phidp_field is None:
    phidp_field = get_field_name("differential_phase")
if kdp_field is None:
    kdp_field = get_field_name("specific_differential_phase")
if unf_field is None:
    unf_field = get_field_name("unfolded_differential_phase")

# 第二部分
# prepare reflectivity field
refl = copy.deepcopy(radar.fields[refl_field]["data"]) + offset
is_low_z = (refl) < low_z
is_high_z = (refl) > high_z
refl[np.where(is_high_z)] = high_z
refl[np.where(is_low_z)] = low_z
z_mod = refl

# 第三部分
# unfold Phi_DP
if debug:
    print("Unfolding")
my_unf = get_phidp_unf(
    radar,
    ncp_lev=min_ncp,
    rhohv_lev=min_rhv,
    debug=debug,
    ncpts=2,
    doc=None,
    sys_phase=sys_phase,
    nowrap=nowrap,
    overide_sys_phase=overide_sys_phase,
    refl_field=refl_field,
    ncp_field=ncp_field,
    rhv_field=rhv_field,
    phidp_field=phidp_field,
)
my_new_ph = copy.deepcopy(radar.fields[phidp_field])
my_unf[:, -1] = my_unf[:, -2] 
my_new_ph["data"] = my_unf
radar.fields.update({unf_field: my_new_ph})

# 小回波的剔除
x_ma = ma.masked_where(notmeteo, my_phidp[radial, :])
try:
    ma.notmasked_contiguous(x_ma)
    for slc in ma.notmasked_contiguous(x_ma):
        # so trying to get rid of clutter and small things that
        # should not add to phidp anyway
        if slc.stop - slc.start < ncpts or slc.start < ncpts:
            x_ma.mask[slc.start - 1 : slc.stop + 1] = True 
    c = 0
except TypeError:  # non sequence, no valid regions
    c = 1  # ie do nothing
    x_ma.mask = True 
except AttributeError:
    # sys.stderr.write('No Valid Regions, ATTERR \n ')
    # sys.stderr.write(myfile.times['time_end'].strftime('%Y-%m-%dT%H:%M:%SZ') + '\n')
    # print x_ma
    # print x_ma.mask
    c = 1  # also do nothing
    x_ma.mask = True

# 第四部分
for sweep in range(len(radar.sweep_start_ray_index["data"])):
    if debug:
        print("Doing ", sweep)
    end_gate, start_ray, end_ray = det_process_range(radar, sweep, fzl, doc=15)
    start_gate = 0 

    A_Matrix = construct_A_matrix(
        len(radar.range["data"][start_gate:end_gate]), St_Gorlv_differential_5pts
    )   

    B_vectors = construct_B_vectors(
        phidp_mod[start_ray:end_ray, start_gate:end_gate],
        z_mod[start_ray:end_ray, start_gate:end_gate],
        St_Gorlv_differential_5pts,
        dweight=self_const,
        coef=coef,
    )   

    weights = np.ones(phidp_mod[start_ray:end_ray, start_gate:end_gate].shape)

    nw = np.bmat([weights, np.zeros(weights.shape)])

    if LP_solver == "pyglpk":
        mysoln = LP_solver_pyglpk(
            A_Matrix, B_vectors, nw, really_verbose=really_verbose
        )   
    elif LP_solver == "cvxopt":
        mysoln = LP_solver_cvxopt(A_Matrix, B_vectors, nw) 
    elif LP_solver == "cylp":
        mysoln = LP_solver_cylp(
            A_Matrix, B_vectors, nw, really_verbose=really_verbose
        )   
    elif LP_solver == "cylp_mp":
        mysoln = LP_solver_cylp_mp(
            A_Matrix, B_vectors, nw, really_verbose=really_verbose, proc=proc
        )   
    else:
        raise ValueError("unknown LP_solver:" + LP_solver)

    proc_ph["data"][start_ray:end_ray, start_gate:end_gate] = mysoln

# 第五部分
# prepare output
sobel = 2.0 * np.arange(window_len) / (window_len - 1.0) - 1.0 
sobel = sobel / (abs(sobel).sum())
sobel = sobel[::-1]
gate_spacing = (radar.range["data"][1] - radar.range["data"][0]) / 1000.0
kdp = scipy.ndimage.convolve1d(proc_ph["data"], sobel, axis=1) / ( 
    (window_len / 3.0) * 2.0 * gate_spacing
)