text:
"""
This module is the computational part of the geometrical module of ToFu
"""

# Built-in
import os
import warnings
from xml.dom import minidom

# Common
import numpy as np
import scipy.interpolate as scpinterp
import scipy.integrate as scpintg
import matplotlib.colors as mplcol
from inspect import signature as insp

# ToFu-specific
try:
    import tofu.geom._def as _def
    import tofu.geom._GG as _GG
except Exception:
    from . import _def as _def
    from . import _GG as _GG

_LTYPES = [int, float, np.int_, np.float_]
_RES = 0.1

###############################################################################
#                            Default parameters
###############################################################################

_SAMPLE_RES = {
    'edge': 0.02,
    'cross': 0.1,
    'surface': 0.1,
    'volume': 0.1,
}
_SAMPLE_RESMODE = {
    'edge': 'abs',
    'cross': 'abs',
    'surface': 'abs',
    'volume': 'abs',
}

def _check_float(var=None, varname=None, vardef=None):
    if var is None:
        var = vardef
    if not type(var) in _LTYPES:
        msg = (
            "Arg {} must be a float!\n".format(varname)
            + "Provided: {}".format(type(var))
        )
        raise Exception(msg)
    return var

###############################################################################
#                            Ves functions
###############################################################################

# ==============================================================================
# Interfacing functions
# ==============================================================================

def _get_pts_from_path_svg(
    path_str=None,
    res=None,
    k0=None,
):

    # Check inputs
    res = _check_float(var=res, varname='res', vardef=_RES)

    # try loading
    try:
        from svg.path import parse_path
    except Exception as err:
        msg = (
            str(err)
            + "\n\nYou do not seem to have svg.path installed\n"
            + "It is an optional dependency only used for this method\n"
            + "To use from_svg(), please install svg.path using:\n"
            + "\tpip install svg.path"
        )
        raise Exception(msg)

    lpath = parse_path(path_str)

    lpath._calc_lengths()
    fract = lpath._fractions

    pos = []
    for ii, pat in enumerate(lpath):
        if pat.__class__.__name__ == 'Line':
            pos.append(np.r_[fract[ii]])
        elif pat.__class__.__name__ == 'Move':
            pos.append(np.r_[fract[ii]])
        elif pat.__class__.__name__ == 'Close':
            pos.append(np.r_[fract[ii]])
        else:
            npts = int(np.ceil(pat.length() / res))
            pos.append(
                np.linspace(fract[ii], fract[ii+1], npts, endpoint=False)
            )

    pos = np.unique(np.concatenate(pos))
    ind1 = np.abs(pos-1.) < 1e-14
    if np.sum(ind1) == 1:
        pos[ind1] = 1.
    elif np.sum(ind1) > 1:
        msg = "Several 1!"
        raise Exception(msg)
    pts = np.array([lpath.point(po) for po in pos])
    pts = np.array([pts.real, pts.imag])

    # Check for reference line
    isref = False
    if 'z' not in path_str.lower():
        if pts.shape[1] == 2:
            isref = True
        elif np.allclose(pts[:, 0], pts[:, -1]):
            pass
        else:
            pts = np.concatenate((pts, pts[:, 0:1]), axis=1)
            msg = (
                f"Non-conform path '{k0}' identified!\n"
                "All path must be either:\n"
                "\t- closed\n"
                "\t- or a unique straight line with 2 points\n"
                "Provided:\n"
                f"path_str: {path_str}\n"
                f"  => closed automatically"
            )
            warnings.warn(msg)

    return pts, isref

def get_paths_from_svg(
    pfe=None,
    res=None,
    r0=None,
    z0=None,
    point_ref1=None,
    point_ref2=None,
    length_ref=None,
    scale=None,
    verb=None,
):

    # check input
    c0 = isinstance(pfe, str) and os.path.isfile(pfe) and pfe.endswith('.svg')
    if not c0:
        msg = (
            "Arg pfe should be a path to a valid .svg file!\n"
            + "Provided:\n\t{}".format(pfe)
        )
        raise Exception(msg)
    pfe = os.path.abspath(pfe)

    # r0, z0, scale
    z0 = _check_float(var=z0, varname='z0', vardef=0.)
    r0 = _check_float(var=r0, varname='r0', vardef=0.)
    scale = _check_float(var=scale, varname='scale', vardef=1.)

    # verb
    if verb is None:
        verb = True
    if not isinstance(verb, bool):
        msg = (
            "Arg verb must be a bool!\n"
            + "Provided:\n\t{}".format(verb)
        )
        raise Exception(msg)

    # Predefine useful var
    doc = minidom.parse(pfe)

    # Try extract raw data
    try:
        dpath = {
            path.getAttribute('id').replace('\n', '').replace('""', ''): {
                'poly': path.getAttribute('d'),
                'color': path.getAttribute('style')
            }
            for path in doc.getElementsByTagName('path')
        }
    except Exception as err:
        msg = (
            "Could not extract path coordinates from {}".format(pfe)
        )
        raise Exception(msg)

    # Derive usable data
    kstr = 'fill:'
    dpath = {k0: v0 for k0, v0 in dpath.items() if kstr in v0['color']}
    lk = list(dpath.keys())
    ref = None
    for ii, k0 in enumerate(lk):

        v0 = dpath[k0]
        poly, isref = _get_pts_from_path_svg(v0['poly'], res=res, k0=k0)
        if isref is True:
            ref = poly
            del dpath[k0]
            continue
        dpath[k0]['poly'] = poly

        # class and color
        color = v0['color'][v0['color'].index(kstr) + len(kstr):].split(';')[0]
        if color == 'none':
            dpath[k0]['cls'] = 'Ves'
            color = None
        else:
            if mplcol.to_rgb(color) == (1., 0., 0.):
                dpath[k0]['cls'] = 'CoilPF'
            else:
                dpath[k0]['cls'] = 'PFC'
        dpath[k0]['color'] = color

    # Check for negative r
    lkneg = [k0 for k0, v0 in dpath.items() if np.any(v0['poly'][0, :] <= 0.)]
    if len(lkneg) > 0.:
        lstr = ['\t- {}'.format(k0) for k0 in lkneg]
        msg = (
            "With the chosen r0 ({}) some structure have negative r values\n"
            + "This is impossible in a toroidal coordinate system\n"
            + "  => the following structures are removed:\n"
            + "\n".join(lstr)
        )
        if len(lkneg) == len(dpath):
            raise Exception(msg)
        else:
            warnings.warn(msg)
        dpath = {k0: dpath[k0] for k0 in dpath.keys() if k0 not in lkneg}

    # Set origin and rescale
    if ref is not None:
        lc = [
            point_ref1 is not None and point_ref2 is not None,
            point_ref1 is not None and length_ref is not None,
        ]
        if not any(lc):
            msg = (
                "Arg reference line for scaling has been detected!\n"
                + "But it cannot be used without providing:\n"
                + "\t- point_ref1 + point_ref2: iterables of len() = 2\n"
                + "\t- point_Ref1 + length_ref: iterable len() = 2 + scalar\n"
            )
            warnings.warn(msg)
        else:
            unit = np.diff(ref, axis=1)
            unit = unit / np.linalg.norm(unit)
            unit = np.array([[unit[0, 0]], [-unit[1, 0]]])
            if not lc[0]:
                point_ref2 = np.array(point_ref1)[:, None] + length_ref*unit

            # if horizontal (resp. vertical line) => coef = inf
            # => assume equal scale for r and z instead to avoid inf
            eps = 1.e-8
            if np.abs(unit[0, 0]) > eps:
                r_coef = (
                    (point_ref2[0]-point_ref1[0]) / (ref[0, 1] - ref[0, 0])
                )
            if np.abs(unit[1, 0]) > eps:
                z_coef = (
                    (point_ref2[1]-point_ref1[1]) / (ref[1, 1] - ref[1, 0])
                )
            if np.abs(unit[0, 0]) < eps:
                # vertical line => assume rscale  = zscale
                r_coef = -z_coef
            if np.abs(unit[1, 0]) < eps:
                # horizontal line => assume rscale  = zscale
                z_coef = -r_coef
            r_offset = point_ref1[0] - r_coef*ref[0, 0]
            z_offset = point_ref1[1] - z_coef*ref[1, 0]

            for k0 in dpath.keys():
                dpath[k0]['poly'] = np.array([
                    r_coef*dpath[k0]['poly'][0, :] + r_offset,
                    z_coef*dpath[k0]['poly'][1, :] + z_offset,
                ])
    else:
        for k0 in dpath.keys():
            dpath[k0]['poly'] = np.array([
                scale*(dpath[k0]['poly'][0, :] - r0),
                scale*(-dpath[k0]['poly'][1, :] - z0),
            ])

    # verb
    if verb is True:
        lVes = sorted([k0 for k0, v0 in dpath.items() if v0['cls'] == 'Ves'])
        lPFC = sorted([k0 for k0, v0 in dpath.items() if v0['cls'] == 'PFC'])
        lCoilPF = sorted([
            k0 for k0, v0 in dpath.items() if v0['cls'] == 'CoilPF'
        ])
        lobj = [
            '\t- {}: {} ({} pts, {})'.format(
                dpath[k0]['cls'], k0,
                dpath[k0]['poly'].shape[1], dpath[k0]['color'],
            )
            for k0 in lVes + lPFC + lCoilPF
        ]
        msg = (
            "The following structures were loaded:\n".format(pfe)
            + "\n".join(lobj)
            + "\nfrom {}".format(pfe)
        )
        print(msg)

    return dpath

# ==============================================================================
# = Ves sub-functions
# ==============================================================================

def _Struct_set_Poly(
    Poly, pos=None, extent=None, arrayorder="C", Type="Tor", Clock=False
):
    """ Compute geometrical attributes of a Struct object """

    # Make Poly closed, counter-clockwise, with '(cc,N)' layout and arrayorder
    try:
        Poly = _GG.format_poly(Poly, order="C", Clock=False, close=True,
                               Test=True)
    except Exception as excp:
        print(excp)
    assert Poly.shape[0] == 2, "Arg Poly must be a 2D polygon !"
    fPfmt = np.ascontiguousarray if arrayorder == "C" else np.asfortranarray

    # Get all remarkable points and moments
    NP = Poly.shape[1] - 1
    P1Max = Poly[:, np.argmax(Poly[0, :])]
    P1Min = Poly[:, np.argmin(Poly[0, :])]
    P2Max = Poly[:, np.argmax(Poly[1, :])]
    P2Min = Poly[:, np.argmin(Poly[1, :])]
    BaryP = np.sum(Poly[:, :-1], axis=1, keepdims=False) / (Poly.shape[1] - 1)
    BaryL = np.array(
        [(P1Max[0] + P1Min[0]) / 2.0, (P2Max[1] + P2Min[1]) / 2.0]
    )
    BaryS, Surf = _GG.poly_area_and_barycenter(Poly, NP)

    # Get lim-related indicators
    noccur = int(pos.size)
    Multi = noccur > 1

    # Get Tor-related quantities
    if Type.lower() == "lin":
        Vol, BaryV = None, None
    else:
        Vol, BaryV = _GG.Poly_VolAngTor(Poly)
        if Vol <= 0.0:
            msg = ("Pb. with volume computation for Struct of type 'Tor' !\n"
                   + "\t- Vol = {}\n".format(Vol)
                   + "\t- Poly = {}\n\n".format(str(Poly))
                   + "  => Probably corrupted polygon\n"
                   + "  => Please check polygon is not self-intersecting")
            raise Exception(msg)

    # Compute the non-normalized vector of each side of the Poly
    Vect = np.diff(Poly, n=1, axis=1)
    Vect = fPfmt(Vect)

    # Compute the normalised vectors directed inwards
    Vin = np.array([Vect[1, :], -Vect[0, :]])
    Vin = -Vin  # Poly is Counter Clock-wise as defined above
    Vin = Vin / np.hypot(Vin[0, :], Vin[1, :])[np.newaxis, :]
    Vin = fPfmt(Vin)

    poly = _GG.format_poly(
        Poly,
        order=arrayorder,
        Clock=Clock,
        close=False,
        Test=True,
    )

    # Get bounding circle
    circC = BaryS
    r = np.sqrt(np.sum((poly - circC[:, np.newaxis]) ** 2, axis=0))
    circr = np.max(r)

    dout = {
        "Poly": poly,
        "pos": pos,
        "extent": extent,
        "noccur": noccur,
        "Multi": Multi,
        "nP": NP,
        "P1Max": P1Max,
        "P1Min": P1Min,
        "P2Max": P2Max,
        "P2Min": P2Min,
        "BaryP": BaryP,
        "BaryL": BaryL,
        "BaryS": BaryS,
        "BaryV": BaryV,
        "Surf": Surf,
        "VolAng": Vol,
        "Vect": Vect,
        "VIn": Vin,
        "circ-C": circC,
        "circ-r": circr,
        "Clock": Clock,
    }
    return dout

def _Ves_get_InsideConvexPoly(
    Poly,
    P2Min,
    P2Max,
    BaryS,
    RelOff=_def.TorRelOff,
    ZLim="Def",
    Spline=True,
    Splprms=_def.TorSplprms,
    NP=_def.TorInsideNP,
    Plot=False,
    Test=True,
):
    if Test:
        assert type(RelOff) is float, "Arg RelOff must be a float"
        assert (
            ZLim is None or ZLim == "Def" or type(ZLim) in [tuple, list]
        ), "Arg ZLim must be a tuple (ZlimMin, ZLimMax)"
        assert type(Spline) is bool, "Arg Spline must be a bool !"
    if ZLim is not None:
        if ZLim == "Def":
            ZLim = (
                P2Min[1] + 0.1 * (P2Max[1] - P2Min[1]),
                P2Max[1] - 0.05 * (P2Max[1] - P2Min[1]),
            )
        indZLim = (Poly[1, :] < ZLim[0]) | (Poly[1, :] > ZLim[1])
        if Poly.shape[1] - indZLim.sum() < 10:
            msg = "Poly seems to be Convex and simple enough !"
            msg += "\n  Poly.shape[1] - indZLim.sum() < 10"
            warnings.warn(msg)
            return Poly
        Poly = np.delete(Poly, indZLim.nonzero()[0], axis=1)
    if np.all(Poly[:, 0] == Poly[:, -1]):
        Poly = Poly[:, :-1]
    Np = Poly.shape[1]
    if Spline:
        BarySbis = np.tile(BaryS, (Np, 1)).T
        Ptemp = (1.0 - RelOff) * (Poly - BarySbis)
        # Poly = BarySbis + Ptemp
        Ang = np.arctan2(Ptemp[1, :], Ptemp[0, :])
        Ang, ind = np.unique(Ang, return_index=True)
        Ptemp = Ptemp[:, ind]
        # spline parameters
        ww = Splprms[0] * np.ones((Np + 1,))
        ss = Splprms[1] * (Np + 1)  # smoothness parameter
        kk = Splprms[2]  # spline order
        nest = int(
            (Np + 1) / 2.0
        )  # estimate of number of knots needed (-1 = maximal)
        # Find the knot points

        # TODO @DV : we can probably get rid of this
        # tckp,uu = scpinterp.splprep([np.append(Ptemp[0,:],Ptemp[0,0]),
        # np.append(Ptemp[1,:],Ptemp[1,0]),np.append(Ang,Ang[0]+2.*np.pi)],
        # w=ww, s=ss, k=kk, nest=nest)
        tckp, uu = scpinterp.splprep(
            [
                np.append(Ptemp[0, :], Ptemp[0, 0]),
                np.append(Ptemp[1, :], Ptemp[1, 0]),
            ],
            u=np.append(Ang, Ang[0] + 2.0 * np.pi),
            w=ww,
            s=ss,
            k=kk,
            nest=nest,
            full_output=0,
        )
        xnew, ynew = scpinterp.splev(np.linspace(-np.pi, np.pi, NP), tckp)
        Poly = np.array([xnew + BaryS[0], ynew + BaryS[1]])
        Poly = np.concatenate((Poly, Poly[:, 0:1]), axis=1)
    if Plot:
        import matplotlib.pyplot as plt

        f = plt.figure(facecolor="w", figsize=(8, 10))
        ax = f.add_axes([0.1, 0.1, 0.8, 0.8])
        ax.plot(Poly[0, :], Poly[1, :], "-k", Poly[0, :], Poly[1, :], "-r")
        ax.set_aspect(aspect="equal", adjustable="datalim"), ax.set_xlabel(
            r"R (m)"
        ), ax.set_ylabel(r"Z (m)")
        f.canvas.draw()
    return Poly

# ==============================================================================
# = Ves sampling functions
# ==============================================================================

def _Ves_get_sample_checkinputs(
    res=None,
    domain=None,
    resMode=None,
    ind=None,
    which='volume'
):
    """ Check inputs for all sampling routines """

    # res
    dres = {'edge': 1, 'cross': 2, 'surface': 2, 'volume': 3}
    if res is None:
        res = _SAMPLE_RES[which]
    ltypes = [int, float, np.int_, np.float_]
    c0 = (type(res) in ltypes
          or (hasattr(res, "__iter__")
              and len(res) == dres[which]
              and all([type(ds) in ltypes for ds in res])))
    if not c0:
        msg = ("Arg res must be either:\n"
               + "\t- float: unique resolution for all directions\n"
               + "\t- iterable of {} floats\n".format(dres[which])
               + "  You provided:\n{}".format(res))
        raise Exception(msg)
    if type(res) in ltypes:
        if which != 'edge':
            res = [float(res) for ii in range(dres[which])]
    else:
        if which == 'edge':
            msg = ("For edge, res cannot be an iterable!\n"
                   + "\t- res: {}".format(res))
            raise Exception(msg)
        res = [float(res[ii]) for ii in range(dres[which])]

    # domain (i.e.: sub-domain to be sampled, defined by its limits)
    ddomain = {'edge': 2, 'cross': 2, 'surface': 3, 'volume': 3}
    if domain is None:
        domain = [None for ii in range(ddomain[which])]
    c0 = (hasattr(domain, "__iter__")
          and len(domain) == ddomain[which]
          and all([dd is None
                   or (hasattr(dd, "__iter__")
                       and len(dd) == 2
                       and all([ss is None
                                or type(ss) in ltypes
                                for ss in dd]))
                   for dd in domain]))
    if not c0:
        msg = ("Arg domain must be a len()={} iterable".format(ddomain[which])
               + " where each element can be:\n"
               + "\t- an iterable 2 floats: [lower, upper] bounds\n"
               + "\t- None: no bounds\n"
               + "  You provided:\n{}".format(domain))
        raise Exception(msg)
    for ii in range(len(domain)):
        if domain[ii] is not None:
            domain[ii] = [float(domain[ii][0])
                          if domain[ii][0] is not None else None,
                          float(domain[ii][1])
                          if domain[ii][1] is not None else None]

    # resMode
    if resMode is None:
        resMode = _SAMPLE_RESMODE[which]
    c0 = isinstance(resMode, str) and resMode.lower() in ["abs", "rel"]
    if not c0:
        msg = ("Arg resMode must be in ['abs','rel']!\n"
               + "  You provided:\n{}".format(resMode))
        raise Exception(msg)
    resMode = resMode.lower()

    # ind (indices of points to be recovered)
    c0 = (ind is None
          or (isinstance(ind, np.ndarray)
              and ind.ndim == 1
              and ind.dtype == np.int_
              and np.all(ind >= 0))
          or (which == 'surface'
              and isinstance(ind, list)
              and all([isinstance(indi, np.ndarray)
                       and indi.ndim == 1
                       and indi.dtype == np.int_
                       and np.all(indi >= 0) for indi in ind])))
    if not c0:
        msg = ("Arg ind must be either:\n"
               + "\t- None: domain is used instead\n"
               + "\t- 1d np.ndarray of positive int: indices\n")
        if which == 'surface':
            msg += "\t- list of 1d np.ndarray of positive indices\n"
        msg += "  You provided:\n{}".format(ind)
        raise Exception(msg)

    return res, domain, resMode, ind

def _Ves_get_sampleEdge(
    VPoly,
    res=None,
    domain=None,
    resMode=None,
    offsetIn=0.0,
    VIn=None,
    margin=1.0e-9
):

    # -------------
    #  Check inputs

    # standard
    res, domain, resMode, ind = _Ves_get_sample_checkinputs(
        res=res,
        domain=domain,
        resMode=resMode,
        ind=None,
        which='edge',
    )

    # specific
    ltypes = [int, float, np.int_, np.float_]
    assert type(offsetIn) in ltypes

    # -------------
    #  Compute

    Pts, reseff, ind, N, Rref, VPolybis = _GG.discretize_vpoly(
        VPoly,
        float(res),
        mode=resMode,
        D1=domain[0],
        D2=domain[1],
        margin=margin,
        DIn=float(offsetIn),
        VIn=VIn,
    )
    return Pts, reseff, ind

def _Ves_get_sampleCross(
    VPoly,
    Min1,
    Max1,
    Min2,
    Max2,
    res=None,
    domain=None,
    resMode=None,
    ind=None,
    margin=1.0e-9,
    mode="flat",
):

    # -------------
    #  Check inputs

    # standard
    res, domain, resMode, ind = _Ves_get_sample_checkinputs(
        res=res,
        domain=domain,
        resMode=resMode,
        ind=ind,
        which='cross',
    )

    # specific
    assert mode in ["flat", "imshow"]

    # -------------
    #  Compute

    MinMax1 = np.array([Min1, Max1])
    MinMax2 = np.array([Min2, Max2])
    if ind is None:
        if mode == "flat":
            Pts, dS, ind, d1r, d2r = _GG.discretize_segment2d(
                MinMax1,
                MinMax2,
                res[0],
                res[1],
                D1=domain[0],
                D2=domain[1],
                mode=resMode,
                VPoly=VPoly,
                margin=margin,
            )
            out = (Pts, dS, ind, (d1r, d2r))
        else:
            x1, d1r, ind1, N1 = _GG._Ves_mesh_dlfromL_cython(
                MinMax1, res[0], domain[0], Lim=True,
                dLMode=resMode, margin=margin
            )
            x2, d2r, ind2, N2 = _GG._Ves_mesh_dlfromL_cython(
                MinMax2, res[1], domain[1], Lim=True,
                dLMode=resMode, margin=margin
            )
            xx1, xx2 = np.meshgrid(x1, x2)
            pts = np.squeeze([xx1, xx2])
            extent = (
                x1[0] - d1r / 2.0,
                x1[-1] + d1r / 2.0,
                x2[0] - d2r / 2.0,
                x2[-1] + d2r / 2.0,
            )
            out = (pts, x1, x2, extent)
    else:
        assert mode == "flat"
        c0 = type(ind) is np.ndarray and ind.ndim == 1
        c0 = c0 and ind.dtype in ["int32", "int64"] and np.all(ind >= 0)
        assert c0, "Arg ind must be a np.ndarray of int !"
        Pts, dS, d1r, d2r = _GG._Ves_meshCross_FromInd(
            MinMax1,
            MinMax2,
            res[0],
            res[1],
            ind,
            dSMode=resMode,
            margin=margin,
        )
        out = (Pts, dS, ind, (d1r, d2r))
    return out

def _Ves_get_sampleS(
    VPoly,
    res=None,
    domain=None,
    resMode="abs",
    ind=None,
    offsetIn=0.0,
    VIn=None,
    VType="Tor",
    VLim=None,
    nVLim=None,
    returnas="(X,Y,Z)",
    margin=1.0e-9,
    Multi=False,
    Ind=None,
):
    """ Sample the surface """

    # -------------
    #  Check inputs

    # standard
    res, domain, resMode, ind = _Ves_get_sample_checkinputs(
        res=res,
        domain=domain,
        resMode=resMode,
        ind=ind,
        which='surface',
    )

    # nVLim and VLim
    if not (type(nVLim) in [int, np.int_] and nVLim >= 0):
        msg = ("Arg nVLim must be a positive int\\n"
               + "  You provided:\n{} ({})".format(nVLim, type(nVLim)))
        raise Exception(msg)
    VLim = None if (VLim is None or nVLim == 0) else np.array(VLim)

    if not isinstance(Multi, bool):
        msg = ("Arg Multi must be a bool!\n"
               + "  You provided:\n{}".format(Multi))
        raise Exception(msg)

    # Check if Multi
    if nVLim > 1:
        assert VLim is not None, "For multiple Struct, Lim cannot be None !"
        assert all([hasattr(ll, "__iter__") and len(ll) == 2 for ll in VLim])
        if Ind is None:
            Ind = np.arange(0, nVLim)
        else:
            Ind = [Ind] if not hasattr(Ind, "__iter__") else Ind
            Ind = np.asarray(Ind).astype(int)
        if ind is not None:
            if isinstance(ind, np.ndarray):
                ind = [ind for ii in range(len(Ind))]
            elif not (isinstance(ind, list) and len(ind) == len(Ind)):
                msg = ("Arg ind must be a list of same len() as Ind!\n"
                       + "\t- provided: {}".format(ind))
                raise Exception(msg)

    else:
        VLim = [None] if VLim is None else [VLim.ravel()]
        if ind is not None:
            if not isinstance(ind, np.ndarray):
                msg = ("ind must be a np.ndarray if nVLim == 1\n"
                       + "\t- provided: {}".format(ind))
                raise Exception(msg)
        Ind = [0]

    if ind is None:
        pts, dS, ind, reseff = (
            [0 for ii in Ind],
            [0 for ii in Ind],
            [0 for ii in Ind],
            [[0, 0] for ii in Ind],
        )
        if VType.lower() == "tor":
            for ii in range(0, len(Ind)):
                if VLim[Ind[ii]] is None:
                    (pts[ii],
                     dS[ii],
                     ind[ii],
                     NL,
                     reseff[ii][0],
                     Rref,
                     reseff[ii][1],
                     nRPhi0,
                     VPbis) = _GG._Ves_Smesh_Tor_SubFromD_cython(
                         res[0],
                         res[1],
                         VPoly,
                         DR=domain[0],
                         DZ=domain[1],
                         DPhi=domain[2],
                         DIn=offsetIn,
                         VIn=VIn,
                         PhiMinMax=None,
                         Out=returnas,
                         margin=margin,
                     )
                else:
                    (pts[ii],
                     dS[ii],
                     ind[ii],
                     NL,
                     reseff[ii][0],
                     Rref,
                     dR0r,
                     dZ0r,
                     reseff[ii][1],
                     VPbis) = _GG._Ves_Smesh_TorStruct_SubFromD_cython(
                         VLim[Ind[ii]],
                         res[0],
                         res[1],
                         VPoly,
                         DR=domain[0],
                         DZ=domain[1],
                         DPhi=domain[2],
                         DIn=offsetIn,
                         VIn=VIn,
                         Out=returnas,
                         margin=margin,
                     )
                    reseff[ii] += [dR0r, dZ0r]
        else:
            for ii in range(0, len(Ind)):
                (pts[ii],
                 dS[ii],
                 ind[ii],
                 NL,
                 reseff[ii][0],
                 Rref,
                 reseff[ii][1],
                 dY0r,
                 dZ0r,
                 VPbis) = _GG._Ves_Smesh_Lin_SubFromD_cython(
                     VLim[Ind[ii]],
                     res[0],
                     res[1],
                     VPoly,
                     DX=domain[0],
                     DY=domain[1],
                     DZ=domain[2],
                     DIn=offsetIn,
                     VIn=VIn,
                     margin=margin,
                 )
                reseff[ii] += [dY0r, dZ0r]
    else:
        ind = ind if Multi else [ind]
        pts, dS, reseff = (
            [np.ones((3, 0)) for ii in Ind],
            [0 for ii in Ind],
            [[0, 0] for ii in Ind],
        )
        if VType.lower() == "tor":
            for ii in range(0, len(Ind)):
                if ind[Ind[ii]].size > 0:
                    if VLim[Ind[ii]] is None:
                        out_loc = _GG._Ves_Smesh_Tor_SubFromInd_cython(
                            res[0],
                            res[1],
                            VPoly,
                            ind[Ind[ii]],
                            DIn=offsetIn,
                            VIn=VIn,
                            PhiMinMax=None,
                            Out=returnas,
                            margin=margin,
                        )
                        pts[ii], dS[ii], NL, reseff[ii][0], Rref = out_loc[:5]
                        reseff[ii][1], nRPhi0, VPbis = out_loc[5:]
                    else:
                        out_loc = _GG._Ves_Smesh_TorStruct_SubFromInd_cython(
                            VLim[Ind[ii]],
                            res[0],
                            res[1],
                            VPoly,
                            ind[Ind[ii]],
                            DIn=offsetIn,
                            VIn=VIn,
                            Out=returnas,
                            margin=margin,
                        )
                        pts[ii], dS[ii], NL, reseff[ii][0], Rref = out_loc[:5]
                        dR0r, dZ0r, reseff[ii][1], VPbis = out_loc[5:]
                        reseff[ii] += [dR0r, dZ0r]
        else:
            for ii in range(0, len(Ind)):
                if ind[Ind[ii]].size > 0:
                    out_loc = _GG._Ves_Smesh_Lin_SubFromInd_cython(
                        VLim[Ind[ii]],
                        res[0],
                        res[1],
                        VPoly,
                        ind[Ind[ii]],
                        DIn=offsetIn,
                        VIn=VIn,
                        margin=margin,
                    )
                    pts[ii], dS[ii], NL, reseff[ii][0], Rref = out_loc[:5]
                    reseff[ii][1], dY0r, dZ0r, VPbis = out_loc[5:]
                    reseff[ii] += [dY0r, dZ0r]

    if len(VLim) == 1:
        pts, dS, ind, reseff = pts[0], dS[0], ind[0], reseff[0]
    return pts, dS, ind, reseff

def _Ves_get_sampleV(
    VPoly,
    Min1,
    Max1,
    Min2,
    Max2,
    res=None,
    domain=None,
    resMode=None,
    ind=None,
    VType="Tor",
    VLim=None,
    returnas="(X,Y,Z)",
    margin=1.0e-9,
    algo="new",
    num_threads=48,
):
    """ Sample the volume """

    # -------------
    #  Check inputs
    res, domain, resMode, ind = _Ves_get_sample_checkinputs(
        res=res,
        domain=domain,
        resMode=resMode,
        ind=ind,
        which='volume',
    )

    # ------------
    # Computation

    MinMax1 = np.array([Min1, Max1])
    MinMax2 = np.array([Min2, Max2])
    VLim = None if VType.lower() == "tor" else np.array(VLim).ravel()
    reseff = [None, None, None]
    if ind is None:
        if VType.lower() == "tor":
            if algo.lower() == "new":
                (pts, dV, ind,
                 reseff[0],
                 reseff[1],
                 reseff[2],
                 sz_r, sz_z,
                 ) = _GG._Ves_Vmesh_Tor_SubFromD_cython(
                    res[0],
                    res[1],
                    res[2],
                    MinMax1,
                    MinMax2,
                    DR=domain[0],
                    DZ=domain[1],
                    DPhi=domain[2],
                    limit_vpoly=VPoly,
                    out_format=returnas,
                    margin=margin,
                    num_threads=num_threads,
                )
            else:
                (pts, dV, ind,
                 reseff[0],
                 reseff[1],
                 reseff[2]) = _GG._Ves_Vmesh_Tor_SubFromD_cython_old(
                    res[0],
                    res[1],
                    res[2],
                    MinMax1,
                    MinMax2,
                    DR=domain[0],
                    DZ=domain[1],
                    DPhi=domain[2],
                    VPoly=VPoly,
                    Out=returnas,
                    margin=margin,
                )

        else:
            (pts, dV, ind,
             reseff[0],
             reseff[1],
             reseff[2]) = _GG._Ves_Vmesh_Lin_SubFromD_cython(
                res[0],
                res[1],
                res[2],
                VLim,
                MinMax1,
                MinMax2,
                DX=domain[0],
                DY=domain[1],
                DZ=domain[2],
                limit_vpoly=VPoly,
                margin=margin,
            )
    else:
        if VType.lower() == "tor":
            if algo.lower() == "new":
                (pts, dV,
                 reseff[0],
                 reseff[1],
                 reseff[2]) = _GG._Ves_Vmesh_Tor_SubFromInd_cython(
                    res[0],
                    res[1],
                    res[2],
                    MinMax1,
                    MinMax2,
                    ind,
                    Out=returnas,
                    margin=margin,
                    num_threads=num_threads,
                )
            else:
                (pts, dV,
                 reseff[0],
                 reseff[1],
                 reseff[2]) = _GG._Ves_Vmesh_Tor_SubFromInd_cython_old(
                    res[0],
                    res[1],
                    res[2],
                    MinMax1,
                    MinMax2,
                    ind,
                    Out=returnas,
                    margin=margin,
                )
        else:
            (pts, dV,
             reseff[0],
             reseff[1],
             reseff[2]) = _GG._Ves_Vmesh_Lin_SubFromInd_cython(
                 res[0],
                 res[1],
                 res[2],
                 VLim,
                 MinMax1,
                 MinMax2,
                 ind,
                 margin=margin
             )
    return pts, dV, ind, reseff

# ==============================================================================
# =  phi / theta projections for magfieldlines
# ==============================================================================

def _Struct_get_phithetaproj(ax=None, poly_closed=None, lim=None, noccur=0):

    # phi = toroidal angle
    if noccur == 0:
        Dphi = np.array([[-np.pi, np.pi]])
        nphi = np.r_[1]
    else:
        assert lim.ndim == 2, str(lim)
        nphi = np.ones((noccur,), dtype=int)
        ind = (lim[:, 0] > lim[:, 1]).nonzero()[0]
        Dphi = np.concatenate((lim, np.full((noccur, 2), np.nan)), axis=1)
        if ind.size > 0:
            for ii in ind:
                Dphi[ii, :] = [lim[ii, 0], np.pi, -np.pi, lim[ii, 1]]
                nphi[ii] = 2

    # theta = poloidal angle
    Dtheta = np.arctan2(poly_closed[1, :] - ax[1], poly_closed[0, :] - ax[0])
    Dtheta = np.r_[np.min(Dtheta), np.max(Dtheta)]
    if Dtheta[0] > Dtheta[1]:
        ntheta = 2
        Dtheta = [Dtheta[0], np.pi, -np.pi, Dtheta[1]]
    else:
        ntheta = 1

    return nphi, Dphi, ntheta, Dtheta

def _get_phithetaproj_dist(
    poly_closed,
    ax,
    Dtheta,
    nDtheta,
    Dphi,
    nDphi,
    theta,
    phi,
    ntheta,
    nphi,
    noccur,
):

    if nDtheta == 1:
        ind = (theta >= Dtheta[0]) & (theta <= Dtheta[1])
    else:
        ind = (theta >= Dtheta[0]) | (theta <= Dtheta[1])

    disttheta = np.full((theta.size,), np.nan)

    # phi within Dphi
    if noccur > 0:
        indphi = np.zeros((nphi,), dtype=bool)
        for ii in range(0, noccur):
            for jj in range(0, nDphi[ii]):
                indphi |= (phi >= Dphi[ii, jj]) & (phi <= Dphi[ii, jj + 1])
        if not np.any(indphi):
            return disttheta, indphi
    else:
        indphi = np.ones((nphi,), dtype=bool)

    # No theta within Dtheta
    if not np.any(ind):
        return disttheta, indphi

    # Check for non-parallel AB / u pairs
    u = np.array([np.cos(theta), np.sin(theta)])
    AB = np.diff(poly_closed, axis=1)
    detABu = AB[0, :, None] * u[1, None, :] - AB[1, :, None] * u[0, None, :]
    inddet = ind[None, :] & (np.abs(detABu) > 1.0e-9)
    if not np.any(inddet):
        return disttheta, indphi

    nseg = poly_closed.shape[1] - 1
    k = np.full((nseg, ntheta), np.nan)

    OA = poly_closed[:, :-1] - ax[:, None]
    detOAu = (OA[0, :, None] * u[1, None, :] - OA[1, :, None] * u[0, None, :])[
        inddet
    ]
    ss = -detOAu / detABu[inddet]
    inds = (ss >= 0.0) & (ss < 1.0)
    inddet[inddet] = inds

    if not np.any(inds):
        return disttheta, indphi

    scaOAu = (OA[0, :, None] * u[0, None, :] + OA[1, :, None] * u[1, None, :])[
        inddet
    ]
    scaABu = (AB[0, :, None] * u[0, None, :] + AB[1, :, None] * u[1, None, :])[
        inddet
    ]
    k[inddet] = scaOAu + ss[inds] * scaABu
    indk = k[inddet] > 0.0
    inddet[inddet] = indk

    if not np.any(indk):
        return disttheta, indphi

    k[~inddet] = np.nan
    indok = np.any(inddet, axis=0)
    disttheta[indok] = np.nanmin(k[:, indok], axis=0)

    return disttheta, indphi

# ==============================================================================
# =  LOS functions
# ==============================================================================
def LOS_PRMin(Ds, us, kOut=None, Eps=1.0e-12, squeeze=True, Test=True):
    """  Compute the point on the LOS where the major radius is minimum """
    if Test:
        assert Ds.ndim in [1, 2, 3] and 3 in Ds.shape and Ds.shape == us.shape
    if kOut is not None:
        kOut = np.atleast_1d(kOut)
        assert kOut.size == Ds.size / 3

    if Ds.ndim == 1:
        Ds, us = Ds[:, None, None], us[:, None, None]
    elif Ds.ndim == 2:
        Ds, us = Ds[:, :, None], us[:, :, None]
    if kOut is not None:
        if kOut.ndim == 1:
            kOut = kOut[:, None]
    _, nlos, nref = Ds.shape

    kRMin = np.full((nlos, nref), np.nan)
    uparN = np.sqrt(us[0, :, :] ** 2 + us[1, :, :] ** 2)

    # Case with u vertical
    ind = uparN > Eps
    kRMin[~ind] = 0.0

    # Else
    kRMin[ind] = (
        -(us[0, ind] * Ds[0, ind] + us[1, ind] * Ds[1, ind]) / uparN[ind] ** 2
    )

    # Check
    kRMin[kRMin <= 0.0] = 0.0
    if kOut is not None:
        kRMin[kRMin > kOut] = kOut[kRMin > kOut]

    # squeeze
    if squeeze:
        if nref == 1 and nlos == 11:
            kRMin = kRMin[0, 0]
        elif nref == 1:
            kRMin = kRMin[:, 0]
        elif nlos == 1:
            kRMin = kRMin[0, :]
    return kRMin

def LOS_CrossProj(
    VType,
    Ds,
    us,
    kOuts,
    proj="All",
    multi=False,
    num_threads=16,
    return_pts=False,
    Test=True,
):
    """ Compute the parameters to plot the poloidal projection of the LOS  """
    assert type(VType) is str and VType.lower() in ["tor", "lin"]
    dproj = {
        "cross": ("R", "Z"),
        "hor": ("x,y"),
        "all": ("R", "Z", "x", "y"),
        "3d": ("x", "y", "z"),
    }
    assert type(proj) in [str, tuple]
    if type(proj) is tuple:
        assert all([type(pp) is str for pp in proj])
        lcoords = proj
    else:
        proj = proj.lower()
        assert proj in dproj.keys()
        lcoords = dproj[proj]
    if return_pts:
        assert proj in ["cross", "hor", "3d"]

    lc = [Ds.ndim == 3, Ds.shape == us.shape]
    if not all(lc):
        msg = "Ds and us must have the same shape and dim in [2,3]:\n"
        msg += "    - provided Ds.shape: %s\n" % str(Ds.shape)
        msg += "    - provided us.shape: %s" % str(us.shape)
        raise Exception(msg)
    lc = [kOuts.size == Ds.size / 3, kOuts.shape == Ds.shape[1:]]
    if not all(lc):
        msg = "kOuts must have the same shape and ndim = Ds.ndim-1:\n"
        msg += "    - Ds.shape    : %s\n" % str(Ds.shape)
        msg += "    - kOutss.shape: %s" % str(kOuts.shape)
        raise Exception(msg)

    # Prepare inputs
    _, nlos, nseg = Ds.shape

    # Detailed sampling for 'tor' and ('cross' or 'all')
    R, Z = None, None
    if "R" in lcoords or "Z" in lcoords:
        angcross = np.arccos(
            np.sqrt(us[0, ...] ** 2 + us[1, ...] ** 2)
            / np.sqrt(np.sum(us ** 2, axis=0))
        )
        resnk = np.ceil(25.0 * (1 - (angcross / (np.pi / 4) - 1) ** 2) + 5)
        resnk = 1.0 / resnk.ravel()

        # Use optimized get sample
        DL = np.vstack((np.zeros((nlos * nseg,), dtype=float), kOuts.ravel()))

        k, reseff, lind = _GG.LOS_get_sample(
            nlos * nseg,
            resnk,
            DL,
            dmethod="rel",
            method="simps",
            num_threads=num_threads,
            Test=Test,
        )

        assert lind.size == nseg * nlos - 1
        ind = lind[nseg - 1 :: nseg]  # noqa
        nbrep = np.r_[lind[0], np.diff(lind), k.size - lind[-1]]
        pts = np.repeat(Ds.reshape((3, nlos * nseg)), nbrep, axis=1) + k[
            None, :
        ] * np.repeat(us.reshape((3, nlos * nseg)), nbrep, axis=1)

        if return_pts:
            pts = np.array([np.hypot(pts[0, :], pts[1, :]), pts[2, :]])
            if multi:
                pts = np.split(pts, ind, axis=1)
            else:
                pts = np.insert(pts, ind, np.nan, axis=1)
        else:
            if multi:
                if "R" in lcoords:
                    R = np.split(np.hypot(pts[0, :], pts[1, :]), ind)
                if "Z" in lcoords:
                    Z = np.split(pts[2, :], ind)
            else:
                if "R" in lcoords:
                    R = np.insert(np.hypot(pts[0, :], pts[1, :]), ind, np.nan)
                if "Z" in lcoords:
                    Z = np.insert(pts[2, :], ind, np.nan)

    # Normal sampling => pts
    # unnecessary only if 'tor' and 'cross'
    x, y, z = None, None, None
    if "x" in lcoords or "y" in lcoords or "z" in lcoords:
        pts = np.concatenate(
            (Ds, Ds[:, :, -1:] + kOuts[None, :, -1:] * us[:, :, -1:]), axis=-1
        )

        if multi:
            ind = np.arange(1, nlos) * (nseg + 1)
            pts = pts.reshape((3, nlos * (nseg + 1)))
        else:
            nancoords = np.full((3, nlos, 1), np.nan)
            pts = np.concatenate((pts, nancoords), axis=-1)
            pts = pts.reshape((3, nlos * (nseg + 2)))

        if return_pts:
            assert proj in ["hor", "3d"]
            if multi:
                if proj == "hor":
                    pts = np.split(pts[:2, :], ind, axis=1)
                else:
                    pts = np.split(pts, ind, axis=1)
            elif proj == "hor":
                pts = pts[:2, :]

        else:
            if multi:
                if "x" in lcoords:
                    x = np.split(pts[0, :], ind)
                if "y" in lcoords:
                    y = np.split(pts[1, :], ind)
                if "z" in lcoords:
                    z = np.split(pts[2, :], ind)
            else:
                if "x" in lcoords:
                    x = pts[0, :]
                if "y" in lcoords:
                    y = pts[1, :]
                if "z" in lcoords:
                    z = pts[2, :]

    if return_pts:
        return pts
    else:
        return R, Z, x, y, z

# ==============================================================================
# =  Meshing & signal
# ==============================================================================
def LOS_get_sample(D, u, dL, DL=None, dLMode="abs", method="sum", Test=True):
    """ Return the sampled line, with the specified method

    'linspace': return the N+1 edges, including the first and last point
    'sum' : return the N middle of the segments
    'simps': return the N+1 egdes, where N has to be even
             (scipy.simpson requires an even number of intervals)
    'romb' : return the N+1 edges, where N+1 = 2**k+1
             (fed to scipy.romb for integration)
    """
    if Test:
        assert all(
            [type(dd) is np.ndarray and dd.shape == (3,) for dd in [D, u]]
        )
        assert not hasattr(dL, "__iter__")
        assert DL is None or all(
            [
                hasattr(DL, "__iter__"),
                len(DL) == 2,
                all([not hasattr(dd, "__iter__") for dd in DL]),
            ]
        )
        assert dLMode in ["abs", "rel"]
        assert type(method) is str and method in [
            "linspace",
            "sum",
            "simps",
            "romb",
        ]
    # Compute the min number of intervals to satisfy the specified resolution
    N = (
        int(np.ceil((DL[1] - DL[0]) / dL))
        if dLMode == "abs"
        else int(np.ceil(1.0 / dL))
    )
    # Modify N according to the desired method
    if method == "simps":
        N = N if N % 2 == 0 else N + 1
    elif method == "romb":
        N = 2 ** int(np.ceil(np.log(N) / np.log(2.0)))

    # Derive k and dLr
    if method == "sum":
        dLr = (DL[1] - DL[0]) / N
        k = DL[0] + (0.5 + np.arange(0, N)) * dLr
    else:
        k, dLr = np.linspace(
            DL[0], DL[1], N + 1, endpoint=True, retstep=True, dtype=float
        )

    Pts = D[:, np.newaxis] + k[np.newaxis, :] * u[:, np.newaxis]
    return Pts, k, dLr

def LOS_calc_signal(
    ff, D, u, dL, DL=None, dLMode="abs", method="romb", Test=True
):
    assert hasattr(ff, "__call__"), (
        "Arg ff must be a callable (function) taking at least 1 positional ",
        "Pts (a (3,N) np.ndarray of cartesian (X,Y,Z) coordinates) !",
    )
    assert not method == "linspace"
    Pts, k, dLr = LOS_get_sample(
        D, u, dL, DL=DL, dLMode=dLMode, method=method, Test=Test
    )
    out = insp(ff)
    N = np.sum(
        [
            (
                pp.kind == pp.POSITIONAL_OR_KEYWORD
                and pp.default is pp.empty
            )
            for pp in out.parameters.values()
        ]
    )

    if N == 1:
        Vals = ff(Pts)
    elif N == 2:
        Vals = ff(Pts, np.tile(-u, (Pts.shape[1], 1)).T)
    else:
        raise ValueError(
            "The function (ff) assessing the emissivity locally "
            + "must take a single positional argument: Pts a (3,N)"
            + " np.ndarray of (X,Y,Z) cartesian coordinates !"
        )

    Vals[np.isnan(Vals)] = 0.0
    if method == "sum":
        Int = np.sum(Vals) * dLr
    elif method == "simps":
        Int = scpintg.simps(Vals, x=None, dx=dLr)
    elif method == "romb":
        Int = scpintg.romb(Vals, dx=dLr, show=False)
    return Int

max_stars_repo_name: WinstonLHS/tofu

max_issues_repo_name: WinstonLHS/tofu


text:
"""
Copyright 2016 Jaehyun Park, 2017 Robin Verschueren

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

import cvxpy.lin_ops.lin_utils as lu

class InverseData(object):
    """Data useful for retrieving a solution from a problem."""

    def __init__(self, problem):
        varis = problem.variables()
        self.id_map, self.var_offsets, self.x_length, self.var_shapes = (
                                                self.get_var_offsets(varis))
        self.id2var = {var.id: var for var in varis}
        self.real2imag = {var.id: lu.get_id() for var in varis
                          if var.is_complex()}
        constr_dict = {cons.id: lu.get_id() for cons in problem.constraints
                       if cons.is_complex()}
        self.real2imag.update(constr_dict)
        self.id2cons = {cons.id: cons for cons in problem.constraints}
        self.cons_id_map = dict()

    def get_var_offsets(self, variables):
        var_shapes = {}
        var_offsets = {}
        id_map = {}
        vert_offset = 0
        for x in variables:
            var_shapes[x.id] = x.shape
            var_offsets[x.id] = vert_offset
            id_map[x.id] = (vert_offset, x.size)
            vert_offset += x.size
        return (id_map, var_offsets, vert_offset, var_shapes)

ext: py

reason: import cvxpy


text:
#!/usr/bin/python

import numpy as np
import math

import astropy.units as u
import astropy.constants as const

def TauAcc(E, B, eff=1):
    # Eq 22 of Dubus 2013
    # E in TeV
    # B in G
    return 0.1 * eff * E / B

def TauIC(E, U, T):
    # E in TeV
    # U in erg / cm3
    # T in K
    tt = const.k_B * T * u.K / (const.m_e * const.c**2)
    gamma = E * u.TeV / (const.m_e * const.c**2)
    g_dot = 5.5e17 * (tt**3) * gamma
    g_dot *= math.log(1 + 0.55 * gamma * tt) / (1 + 25 * tt * gamma)
    g_dot *= 1 + 1.4 * gamma * tt / (1 + 12 * (gamma**2)*(tt**2))
    Ufac = (U * u.erg / u.cm**3) / (4 * const.sigma_sb * (T * u.K)**4 / const.c)
    return gamma.to(1) / g_dot.to(1) / Ufac

def TauSyn(E, B):
    # E in TeV
    # B in G
    # From Eq 23 of Dubus 2013
    return 400 * (1. / E) * (1 / B)**2

def Ebreak(B, T, U):
    # B in G
    # T in K
    # U in erg /cm**3
    ####
    # 1st step
    en = np.logspace(-7, 3, 100)
    tdif = [abs(TauSyn(E=e, B=B) - TauIC(E=e, T=T, U=U)) for e in en]
    idx = np.argmin(tdif)
    ebr = en[idx]
    # 2nd step
    lge = math.log10(ebr)
    en = np.logspace(lge - 0.05, lge + 0.05, 100)
    tdif = [abs(TauSyn(E=e, B=B) - TauIC(E=e, T=T, U=U)) for e in en]
    idx = np.argmin(tdif)
    ebr = en[idx]
    # 3rd step
    lge = math.log10(ebr)
    en = np.logspace(lge - 0.0005, lge + 0.0005, 100)
    tdif = [abs(TauSyn(E=e, B=B) - TauIC(E=e, T=T, U=U)) for e in en]
    idx = np.argmin(tdif)
    ebr = en[idx]
    return en[idx] if tdif[idx] < 10 else 0

def Emax(B, eff=1):
    # Eq. 24 from Dubus 2013
    return 60 * math.pow(eff * B, -0.5)

def Esyn(E, B):
    # E in TeV
    # B in G
    return 5 * 1e-8 * B * (E**2)

def bracket_sig(x):
    r = 1 - 2 * (x + 1) / x**2
    r *= math.log(2 * x + 1)
    r += 0.5 + 4 / x - 1 / (2 * ((2*x + 1)**2))
    return r

def sigma_kn(E, T):
    eph = 2.7 * const.k_B * T * u.K
    gamma = E * u.TeV / (const.m_e * const.c**2)
    x = eph * gamma / (const.m_e * const.c**2)
    return (3./8) * const.sigma_T * bracket_sig(x.to(1)) / x.to(1)
    # return (3./8) * const.sigma_T * (math.log(2*x.to(1))+0.5) / x.to(1)

def tau_ic_old(E, U):
    t = 3e7 * u.yr.to(u.s)
    t /= U * u.erg.to(u.eV)
    t /= 100 * E
    return t

def tau_kn_new(E, U, T):
    gamma = E * u.TeV / (const.m_e * const.c**2)
    t = tau_ic(E=E, U=U) * const.sigma_T / sigma_kn(E=E, T=T)
    #
    # eph = 2.7 * const.k_B * T * u.K
    # t = 1e6 * u.yr.to(u.s) * const.sigma_T / sigma_kn(E=E, T=T)
    # t /= U / eph.to(u.erg).value
    #
    # eph = 2.7 * const.k_B * T * u.K
    # t = 1e-20 * 3e7 * u.yr.to(u.s) / U / eph.to(u.eV).value
    # t *= const.sigma_T / sigma_kn(E=E, T=T)
    return t

def tau_kn(E, T, U):
    # From BG70
    # E in TeV
    # T in K
    # U in erg / c**3
    # return t in s
    de = (const.sigma_T / 16) * (const.m_e * const.c * const.k_B * T * u.K)**2 / const.hbar**3
    a = 4 * E * u.TeV * const.k_B * T * u.K / (const.m_e**2 * const.c**4)
    de *= (math.log(a.to(1)) - 5/6. - 1.147)
    # The calculation by BG70 includes the energy density
    # In order to multiply by the energy density, we need
    # to correct the unit
    # density_unit is the unit of density from BG70 calculations
    # that we need to convert to erg / cm3
    # (This was cross-checked with Eq 28 of Dubus 2013)
    density_unit = const.c * const.k_B**2 * u.K**2 * const.m_e**2 / (const.hbar**3)
    fac = density_unit.to(u.erg / u.cm**3).value
    de *= U / fac
    t = E * u.TeV / de
    return t.to(u.s).value

max_stars_repo_name: RaulRPrado/tev-binaries-model

max_issues_repo_name: RaulRPrado/tev-binaries-model

max_forks_repo_name: RaulRPrado/tev-binaries-model

reason: import numpy,import astropy


text:
from sklearn.metrics import f1_score
from sklearn.metrics import jaccard_score
from sklearn.metrics import accuracy_score

from sklearn.neighbors import KernelDensity
import numpy as np

def _kerneldensity(*, train, test, x_predict=None, metrics, bandwidth=1.0, algorithm='auto', kernel='gaussian', metric='euclidean', atol=0, rtol=0, breadth_first=True, leaf_size=40, metric_params=None):
    """
    For more info visit :
    https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html#sklearn.neighbors.KernelDensity
    """

    model = KernelDensity(bandwidth=bandwidth, algorithm=algorithm, kernel=kernel, metric=metric, atol=atol, rtol=rtol,
breadth_first=breadth_first, leaf_size=leaf_size, metric_params=metric_params)
    model.fit(train[0], train[1])
    model_name = 'Kernel Density'
    y_hat = model.predict(test[0])

    if metrics == 'accuracy':
        accuracy = accuracy_score(test[1], y_hat)

    if metrics == 'f1':
        accuracy = f1_score(test[1], y_hat)

    if metrics == 'jaccard':
        accuracy = jaccard_score(test[1], y_hat)

    if x_predict is None:
        return (model_name, accuracy, None)

    y_predict = model.predict(x_predict)
    return (model_name, accuracy, y_predict)

max_stars_repo_name: kritikaparmar-programmer/DataScientist

max_issues_repo_name: kritikaparmar-programmer/DataScientist

max_forks_repo_name: kritikaparmar-programmer/DataScientist

reason: import numpy


text:
import string
import re
from numpy import array, argmax, random, take
import pandas as pd
from keras.models import Sequential
from keras.layers import Dense, LSTM, Embedding, RepeatVector
from keras.preprocessing.text import Tokenizer
from keras.callbacks import ModelCheckpoint
from keras.preprocessing.sequence import pad_sequences
from keras.models import load_model
from keras import optimizers
import matplotlib.pyplot as plt
%matplotlib inline
pd.set_option('display.max_colwidth', 200)

# function to read raw text file
def read_text(filename):
        # open the file
        file = open(filename, mode='rt', encoding='utf-8')

        # read all text
        text = file.read()
        file.close()
        return text
# split a text into sentences
def to_lines(text):
      sents = text.strip().split('\n')
      sents = [i.split('\t') for i in sents]
      return sents

data = read_text("deu.txt")
deu_eng = to_lines(data)
deu_eng = array(deu_eng)

deu_eng = deu_eng[:50000,:]

deu_eng

# Remove punctuation
deu_eng[:,0] = [s.translate(str.maketrans('', '', string.punctuation)) for s in deu_eng[:,0]]
deu_eng[:,1] = [s.translate(str.maketrans('', '', string.punctuation)) for s in deu_eng[:,1]]

deu_eng

# convert text to lowercase
for i in range(len(deu_eng)):
    deu_eng[i,0] = deu_eng[i,0].lower()
    deu_eng[i,1] = deu_eng[i,1].lower()

deu_eng

# empty lists
eng_l = []
deu_l = []

# populate the lists with sentence lengths
for i in deu_eng[:,0]:
      eng_l.append(len(i.split()))

for i in deu_eng[:,1]:
      deu_l.append(len(i.split()))

length_df = pd.DataFrame({'eng':eng_l, 'deu':deu_l})

length_df.hist(bins = 30)
plt.show()

# function to build a tokenizer
def tokenization(lines):
      tokenizer = Tokenizer()
      tokenizer.fit_on_texts(lines)
      return tokenizer

# prepare english tokenizer
eng_tokenizer = tokenization(deu_eng[:, 0])
eng_vocab_size = len(eng_tokenizer.word_index) + 1

eng_length = 8
print('English Vocabulary Size: %d' % eng_vocab_size)

# prepare Deutch tokenizer
deu_tokenizer = tokenization(deu_eng[:, 1])
deu_vocab_size = len(deu_tokenizer.word_index) + 1

deu_length = 8
print('Deutch Vocabulary Size: %d' % deu_vocab_size)

def encode_sequences(tokenizer, length, lines):
         # integer encode sequences
         seq = tokenizer.texts_to_sequences(lines)
         # pad sequences with 0 values
         seq = pad_sequences(seq, maxlen=length, padding='post')
         return seq
from sklearn.model_selection import train_test_split

# split data into train and test set
train, test = train_test_split(deu_eng, test_size=0.2, random_state = 12)

# prepare training data
trainX = encode_sequences(deu_tokenizer, deu_length, train[:, 1])
trainY = encode_sequences(eng_tokenizer, eng_length, train[:, 0])

# prepare validation data
testX = encode_sequences(deu_tokenizer, deu_length, test[:, 1])
testY = encode_sequences(eng_tokenizer, eng_length, test[:, 0])

def define_model(in_vocab,out_vocab, in_timesteps,out_timesteps,units):
      model = Sequential()
      model.add(Embedding(in_vocab, units, input_length=in_timesteps, mask_zero=True))
      model.add(LSTM(units))
      model.add(RepeatVector(out_timesteps))
      model.add(LSTM(units, return_sequences=True))
      model.add(Dense(out_vocab, activation='softmax'))
      return model

model = define_model(deu_vocab_size, eng_vocab_size, deu_length, eng_length, 512)
rms = optimizers.RMSprop(lr=0.001)
model.compile(optimizer=rms, loss='sparse_categorical_crossentropy')

filename = 'model.h1.Akbar'
checkpoint = ModelCheckpoint(filename, monitor='val_loss', verbose=1, save_best_only=True, mode='min')

# train model
history = model.fit(trainX, trainY.reshape(trainY.shape[0], trainY.shape[1], 1),
                    epochs=30, batch_size=512, validation_split = 0.2,callbacks=[checkpoint],
                    verbose=1)

plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.legend(['train','validation'])
plt.show()

model = load_model('model.h1.Akbar')
preds = model.predict_classes(testX.reshape((testX.shape[0],testX.shape[1])))

def get_word(n, tokenizer):
      for word, index in tokenizer.word_index.items():
          if index == n:
              return word
      return None

preds_text = []
for i in preds:
       temp = []
       for j in range(len(i)):
            t = get_word(i[j], eng_tokenizer)
            if j > 0:
                if (t == get_word(i[j-1], eng_tokenizer)) or (t == None):
                     temp.append('')
                else:
                     temp.append(t)
            else:
                   if(t == None):
                          temp.append('')
                   else:
                          temp.append(t)

       preds_text.append(' '.join(temp))

pred_df = pd.DataFrame({'actual' : test[:,0], 'predicted' : preds_text})

pred_df.sample(20)

reason: from numpy


text:
""" Contains the functions to be used to produce a fit from a taught Gaussian process or ROM.

**Contents**:
    **Input_Grid**: A function that returns a grid of known inputs -- Can be computationally expensive for high dimensions.

    **Generate_Inputs**: A function used to generate many inputs.

    **Test_Data**: Loads test_data that can be used as inputs for predictions and to test predictions.

    **Rotate_Inputs**: Rotates the standardized inputs by Theta ready to produce a fit.

    **Predict**: Prediction using a GaussianBundle.

    **Create_CSV**: Saves the inputs, rotated inputs and predictions into a CSV file.

    **Create_CSV_with_Observed**: Saves the inputs, rotated inputs, predictions and observed output into a CSV file.
"""

import numpy as np
import pandas as pd
from pathlib import Path
from ROMCOMMA.archivedtst.romcomma.data import Store, Fold, Frame
import ROMCOMMA.archivedtst.romcomma.model as model
from ROMCOMMA.archivedtst.romcomma.typing_ import NP, Union, Tuple
from shutil import rmtree, copytree
from pyDOE import lhs
from scipy.stats.distributions import norm
from sklearn.metrics import r2_score
import time

def Input_Grid(bound: float = 1.0, n_samples: int = 3, dim: int = 2) -> NP.Matrix:
    """ A function to generate a mesh grid of inputs that can be used to produce a fit.

    Args:
        bound: The upper and lower bound in which the inputs can be.
        n_samples: The number of values between the upper and lower bound.
        dim: The dimensions of an input grid.
    Returns: A NP.matrix of size (dim x n_samples**dim)
    """
    space = np.linspace(-bound, bound, n_samples)
    a_tuple = (space,)*dim
    mesh = np.ravel(np.meshgrid(*a_tuple))
    mesh.shape = (dim, n_samples**dim)
    Stand_Inputs = np.transpose(mesh)
    return Stand_Inputs

def Generate_Inputs(M: int, N: int, mean: float = 0, std: float = 0.5) -> NP.Matrix:
    """ Generates a large number of inputs.

    Args:
        M: The dimensions of the inputs.
        N: The number of samples / amount of inputs generated.
        mean: The samples will be standardized around this mean. The default is 0.
        std: The samples will be standardized by this standard deviation. 95% of the samples will have be within 2 stdv's from the mean.
    Returns: A NP.matrix of size (M x N)
    """
    Stand_Inputs = norm(loc=mean, scale=std).ppf(lhs(M, samples=N, criterion='cm'))
    return Stand_Inputs

def Test_Data(gp_path: Union[str, Path], test_data_file: str, standard_data_file: str) -> Tuple[NP.Matrix, NP.Vector]:
    """ Loads and saves the test data as standardized inputs (X) and standardized observed outputs (Y).

    Args:
        gp_path: Path to the folder (inside a split) where the test data and standard data files are stored.
        test_data_file: The name of the test data file, e.g. 'test_data.csv'.
        standard_data_file: The name of the standard data file which contains the mean and standard deviation used to standardize the test data.
    Returns: A tuple containing the NP.Matrix of inputs (X) and the NP.Vector of outputs (Y).
    """
    GP_PATH = Path(gp_path)
    test_data = pd.read_csv(GP_PATH / test_data_file, index_col=0)
    mean_std = pd.read_csv(GP_PATH / standard_data_file, index_col=0)
    Dataset = test_data.values
    AVG_STD = mean_std.values
    Stand_Dataset = (Dataset - AVG_STD[0]) / AVG_STD[1]
    Stand_Inputs = Stand_Dataset[:, :-1]
    Stand_Observed_Outputs = Stand_Dataset[:, -1]
    return Stand_Inputs, Stand_Observed_Outputs

def Rotate_Inputs(gp_path: Union[str, Path], Stand_Inputs: NP.Matrix, Us_taken: int, k: int = 0) -> NP.Matrix:
    """ Rotates the standardized inputs by theta to produce the rotated inputs that can be used when predicting with a ROM.

    Args:
        gp_path: Path to a model.GaussianBundle. The extension of this filename is the number of input dimensions M.
            An extension of 0 or a missing extension means full order, taking M from the training data.
        Stand_Inputs: The standardized input values - an (N,M) numpy array, consisting of N test inputs, each of dimension M.
        Us_taken: The amount of Us taken from the ROM to be used in predicting.
        k: The fold that contains the data that has been used to train the GP's.
    Returns: The rotated inputs, U, that can be used for predicting using a ROM - a numpy array of dimensions (N x Xs_taken).
    """
    GP_PATH = Path(gp_path)
    Mu = Us_taken
    Split_Store = Store(GP_PATH)
    fold = Fold(Split_Store, k)
    fold_dir = fold.dir
    sobol = fold_dir / "ROM.optimized" / "sobol"
    theta_T = np.transpose(Frame(sobol / "Theta.csv", csv_parameters={'header': [0]}).df.values)
    k = int(fold_dir.suffix[1:])
    M = Fold(fold_dir.parent, k, Mu).M +1
    if 0 < Mu < M:
        Rotated_Inputs = Stand_Inputs @ theta_T[:, 0:Mu]
    else:
        Rotated_Inputs = Stand_Inputs @ theta_T
    return Rotated_Inputs

def Predict(gp_path: Union[str, Path], gb_source: str, gb_destination: str, Rotated_Inputs: NP.Matrix, Us_taken: int ,k: int = 0) -> Tuple[NP.Vector, NP.Vector]:
    """ Prediction using a GaussianBundle.

    Args:
        gp_path: The path that takes you into the trained Gaussian process' Split folder where the Fold can be found.
        gb_source: The name of the trained GB e.g. "ROM.optimized" or "ARD".
        gb_destination: The name of the new folder where the GB will be copied and predictions made e.g. "ROM_Predictions".
        Rotated_Inputs: The rotated inputs that are predicted on.
        Us_taken: The amount of Us taken from the ROM to be used in predicting.
        k: The fold that contains the data that has been used to train the GP's.
    Returns: A pair (predictive_mean, predictive_std) of (N,1) numpy arrays.

    """
    GP_PATH = Path(gp_path)
    Split_Store = Store(GP_PATH)
    fold = Fold(Split_Store, k, Us_taken)
    gb_dir = fold.dir / gb_destination
    rmtree(gb_dir, ignore_errors=True)
    copytree(src=fold.dir / gb_source, dst=gb_dir)
    gb = model.mygpy.GaussianBundle(fold=fold, name=gb_destination, parameters=None, overwrite=True, reset_log_likelihood=False)
    Predicted_Output = gb.predict(Rotated_Inputs)
    return Predicted_Output

def Create_CSVs(gp_path: Union[str, Path], gb_destination: str, Predicted_Outputs: Tuple[NP.Vector, NP.Vector], Stand_Inputs: NP.Matrix,
                Rotated_Inputs: NP.Matrix, k: int = 0, csv_name: str = "__predictions__"):
    """ Saves the inputs, rotated inputs and predictions into a CSV file.

    Args:
        gp_path: The path that takes you into the trained Gaussian process' Split folder where the Fold can be found.
        gb_destination: The name of the new folder where the GB will be copied and predictions made e.g. "ROM_Predictions".
        Predicted_Outputs: A pair (predictive_mean, predictive_std) of (N,1) numpy arrays.
        Stand_Inputs: The standardized input values - an (N,M) numpy array, consisting of N test inputs, each of dimension M.
        Rotated_Inputs: The rotated inputs that are predicted on.
        k: The fold that contains the data that has been used to train the GP's.
        csv_name: The name of the csv file that has been created. Default is "__predictions__".
    Returns: A frame containing the inputs, the rotated inputs and predictions for each split directory.

    """
    GP_PATH = Path(gp_path)
    Split_Store = Store(GP_PATH)
    fold = Fold(Split_Store, k)
    gb_dir = fold.dir / gb_destination
    mean_prediction_T = np.transpose(Predicted_Outputs[0, :, [0]])
    std_prediction_T = np.sqrt(np.transpose(Predicted_Outputs[1, :, [0]]))
    df_X = pd.DataFrame(Stand_Inputs)
    inputs_label = ['x' + str(i - 1) for i in np.arange(1, len(df_X.columns) + 1)]
    inputs_dict = dict(zip(df_X.columns, inputs_label))
    df_X = df_X.rename(columns=inputs_dict)
    DF_U = pd.DataFrame(Rotated_Inputs)
    rotated_label = ['u' + str(i - 1) for i in np.arange(1, len(DF_U.columns) + 1)]
    rotated_dict = dict(zip(DF_U.columns, rotated_label))
    DF_U = DF_U.rename(columns=rotated_dict)
    df_mean = pd.DataFrame(mean_prediction_T)
    mean_label = ['Predicted_Mean']
    mean_dict = dict(zip(df_mean.columns, mean_label))
    df_mean = df_mean.rename(columns=mean_dict)
    df_std = pd.DataFrame(std_prediction_T)
    std_label = ['Predicted_STD']
    std_dict = dict(zip(df_std.columns, std_label))
    df_std = df_std.rename(columns=std_dict)
    df = pd.concat([df_X, DF_U, df_mean, df_std], axis=1)
    frame = Frame(gb_dir / csv_name, df)
    return frame

def Create_CSV_with_Observed(gp_path: Union[str, Path], gb_destination: str, Predicted_Outputs: Tuple[NP.Vector, NP.Vector], Stand_Inputs: NP.Matrix,
                             Rotated_Inputs: NP.Matrix, Stand_Observed_Outputs: NP.Vector, k: int = 0, csv_name: str = "__predictions__"):
    """ Saves the inputs, rotated inputs, predictions and observed output into a CSV file.

    Args:
        gp_path: The path that takes you into the trained Gaussian process' Split folder where the Fold can be found.
        gb_destination: The name of the new folder where the GB will be copied and predictions made e.g. "ROM_Predictions".
        Predicted_Outputs: A pair (predictive_mean, predictive_std) of (N,1) numpy arrays.
        Stand_Inputs: The standardized input values - an (N,M) numpy array, consisting of N test inputs, each of dimension M.
        Rotated_Inputs: The rotated inputs that are predicted on.
        Stand_Observed_Outputs: A NP.Vector of standardized observed outputs from the inputs.
        k: The fold that contains the data that has been used to train the GP's.
        csv_name: The name of the csv file that has been created. Default is "__predictions__".
    Returns: A frame containing the inputs, the rotated inputs, the predictions and the observed outputs.

    """
    GP_PATH = Path(gp_path)
    Split_Store = Store(GP_PATH)
    fold = Fold(Split_Store, k)
    gb_dir = fold.dir / gb_destination
    mean_prediction_T = np.transpose(Predicted_Outputs[0, :, [0]])
    std_prediction_T = np.sqrt(np.transpose(Predicted_Outputs[1, :, [0]]))
    df_X = pd.DataFrame(Stand_Inputs)
    inputs_label = ['x' + str(i - 1) for i in np.arange(1, len(df_X.columns) + 1)]
    inputs_dict = dict(zip(df_X.columns, inputs_label))
    df_X = df_X.rename(columns=inputs_dict)
    DF_U = pd.DataFrame(Rotated_Inputs)
    rotated_label = ['u' + str(i - 1) for i in np.arange(1, len(DF_U.columns) + 1)]
    rotated_dict = dict(zip(DF_U.columns, rotated_label))
    DF_U = DF_U.rename(columns=rotated_dict)
    df_mean = pd.DataFrame(mean_prediction_T)
    mean_label = ['Predicted_Mean']
    mean_dict = dict(zip(df_mean.columns, mean_label))
    df_mean = df_mean.rename(columns=mean_dict)
    df_std = pd.DataFrame(std_prediction_T)
    std_label = ['Predicted_STD']
    std_dict = dict(zip(df_std.columns, std_label))
    df_std = df_std.rename(columns=std_dict)
    df_Y = pd.DataFrame(Stand_Observed_Outputs)
    Y_label = ['Observed_Output']
    Y_dict = dict(zip(df_Y.columns, Y_label))
    df_Y = df_Y.rename(columns=Y_dict)
    df = pd.concat([df_X, DF_U, df_mean, df_std, df_Y], axis=1)
    frame = Frame(gb_dir / csv_name, df)
    return frame

if __name__ == '__main__':
    start = time.time()
    Base_Path = Path('X:\\comma_group1\\Rom\\dat\\AaronsTraining\\MOF\\MOF_Data_Reduced_6D\\MOF_Rotated_1.1')
    store = Store(Base_Path)
    splits = store.splits
    for split_dir in splits:
        split_start = time.time()
        split = Store(split_dir)
        GP_Path = Path(Base_Path / split_dir)
        Mu = 6
        X, Y = Test_Data(gp_path=GP_Path, test_data_file="__test_data__.csv", standard_data_file="__test_standard__.csv")
        U = Rotate_Inputs(gp_path=GP_Path, Stand_Inputs=X, Us_taken=Mu)
        Prediction = Predict(gp_path=GP_Path, gb_source="ROM.optimized", gb_destination="ROM.Predictions", Rotated_Inputs=U, Us_taken=Mu)
        predict_frame = Create_CSV_with_Observed(gp_path=GP_Path, gb_destination="ROM.Predictions", Predicted_Outputs=Prediction,
                                                 Stand_Inputs=X, Rotated_Inputs=U, Stand_Observed_Outputs=Y, csv_name="__predictions.Mu.6__")
        print('Goodness of fit = ', r2_score(Y, np.transpose(Prediction[0, :, [0]])))
        split_time_mins = (time.time() - split_start) / 60
        print(split_dir.name, "has finished in {:.2f} minutes.".format(split_time_mins))
    print()
    time_mins = (time.time() - start) / 60
    print("It took {:.2f} minutes to execute this code.".format(time_mins))

max_stars_repo_name: judejeh/rom-comma

max_issues_repo_name: judejeh/rom-comma

max_forks_repo_name: judejeh/rom-comma

reason: import numpy,from scipy


text:
"""
Module for computing the Local Binary Pattern in EOPatch

Credits:
Copyright (c) 2018-2019 Hugo Fournier (Magellium)
Copyright (c) 2017-2019 Matej Aleksandrov, Devis Peresutti (Sinergise)

This source code is licensed under the MIT license found in the LICENSE
file in the root directory of this source tree.
"""

import skimage.feature
import numpy as np

from eolearn.core import EOTask, FeatureType

class LocalBinaryPatternTask(EOTask):
    """ Task to compute the Local Binary Pattern images

    LBP looks at points surrounding a central point and tests whether the surrounding points are greater than or
    less than the central point

    The task uses skimage.feature.local_binary_pattern to extract the texture features.
    """
    def __init__(self, feature, nb_points=24, radius=3):
        """
        :param feature: A feature that will be used and a new feature name where data will be saved. If new name is not
            specified it will be saved with name '<feature_name>_LBP'.

            Example: (FeatureType.DATA, 'bands') or (FeatureType.DATA, 'bands', 'lbp')
        :param nb_points: Number of point to use
        :type nb_points: int
        :param radius: Radius of the circle of neighbors
        :type radius: int
        """
        self.feature = self._parse_features(feature, default_feature_type=FeatureType.DATA, new_names=True,
                                            rename_function='{}_LBP'.format)

        self.nb_points = nb_points
        self.radius = radius
        if nb_points < 1 or radius < 1:
            raise ValueError('Local binary pattern task parameters must be positives')

    def _compute_lbp(self, data):
        result = np.empty(data.shape, dtype=float)
        for time in range(data.shape[0]):
            for band in range(data.shape[-1]):
                image = data[time, :, :, band]
                result[time, :, :, band] = skimage.feature.local_binary_pattern(image, self.nb_points, self.radius,
                                                                                method='uniform')
        return result

    def execute(self, eopatch):
        """ Execute computation of local binary patterns on input eopatch

            :param eopatch: Input eopatch
            :type eopatch: eolearn.core.EOPatch
            :return: EOPatch instance with new key holding the LBP image.
            :rtype: eolearn.core.EOPatch
        """
        for feature_type, feature_name, new_feature_name in self.feature:
            eopatch[feature_type][new_feature_name] = self._compute_lbp(eopatch[feature_type][feature_name])

        return eopatch

max_stars_repo_name: mohammadrezabk/eo-learn

max_issues_repo_name: mohammadrezabk/eo-learn

max_forks_repo_name: mohammadrezabk/eo-learn


text:
"""
(weighted) random walks for walk-based NE methods:
DeepWalk, Node2Vec, TriDNR, and ABRW;
alias sampling; walks by multiprocessors; etc.

by Chengbin Hou & Zeyu Dong 2018
"""

import random
import time

import numpy as np
from networkx import nx
from scipy import sparse

# ===========================================ABRW-weighted-walker============================================
class WeightedWalker:
    ''' Weighted Walker for Attributed Biased Randomw Walks (ABRW) method
    '''

    def __init__(self, node_id_map, transition_mat, workers):
        assert sparse.isspmatrix_csr(transition_mat)  # currently support csr format
        self.T = transition_mat
        self.look_back_list = node_id_map
        self.workers = workers

    # alias sampling for ABRW-------------------------
    def simulate_walks(self, num_walks, walk_length):
        t1 = time.time()
        print('333')
        self.preprocess_transition_probs()  # construct alias table; adapted from node2vec
        print('444')
        t2 = time.time()
        print(f'Time for construct alias table: {(t2-t1):.2f}')

        walks = []
        nodes = [i for i in range(self.T.shape[0])]
        for walk_iter in range(num_walks):
            t1 = time.time()
            random.shuffle(nodes)
            for node in nodes:
                walks.append(self.weighted_walk(walk_length=walk_length, start_node=node))
            t2 = time.time()
            print(f'Walk iteration: {walk_iter+1}/{num_walks}; time cost: {(t2-t1):.2f}')

        for i in range(len(walks)):  # use ind to retrive original node ID
            for j in range(len(walks[0])):
                walks[i][j] = self.look_back_list[int(walks[i][j])]
        return walks

    def weighted_walk(self, walk_length, start_node):  # more efficient way instead of copy from node2vec
        walk = [start_node]
        while len(walk) < walk_length:
            cur = walk[-1]  # the last node
            cur_nbrs = self.T.getrow(cur).nonzero()[1]
            if len(cur_nbrs) > 0:  # if non-isolated node
                walk.append( cur_nbrs[ alias_draw(self.alias_nodes[cur][0], self.alias_nodes[cur][1]) ] )  # alias sampling in O(1) time to get the index of
            else:  # if isolated node                                                                  # 1) randomly choose a nbr; 2) judge if use nbr or its alias
                break
        return walk

    def preprocess_transition_probs(self):
        T = self.T
        alias_nodes = {}
        # alias_nodes = [alias_setup(T.getrow(i).data) for i in range(T.shape[0])]

        for i in range(T.shape[0]):
            sparse_row = T.getrow(i)
            # r_index, c_index = sparse_row.nonzero()
            # c_index = sparse_row.nonzero()[1]
            c_value = sparse_row.data
            # assert np.sum(c_value) == 1.0  # error if not prob dist
            # print(np.sum(c_value))
            alias_node = alias_setup(c_value)  # where array0 gives alias node indexes; array1 gives its prob
            # alias_index = [c_index[j] for j in alias_prob_index] # use alias index to find matrix index; and replace
            # alias_nodes[i] = (alias_index, alias_prob)
            alias_nodes[i] = alias_node

        self.alias_nodes = alias_nodes

def deepwalk_walk_wrapper(class_instance, walk_length, start_node):
    class_instance.deepwalk_walk(walk_length, start_node)

# ===========================================deepWalk-walker============================================

class BasicWalker:
    def __init__(self, g, workers):
        self.g = g
        self.node_size = g.get_num_nodes()
        self.look_up_dict = g.look_up_dict

    def deepwalk_walk(self, walk_length, start_node):
        '''
        Simulate a random walk starting from start node.
        '''
        G = self.g.G

        walk = [start_node]

        while len(walk) < walk_length:
            cur = walk[-1]
            cur_nbrs = list(G.neighbors(cur))
            if len(cur_nbrs) > 0:
                walk.append(random.choice(cur_nbrs))
            else:
                break
        return walk

    def simulate_walks(self, num_walks, walk_length):
        '''
        Repeatedly simulate random walks from each node.
        '''
        G = self.g.G
        walks = []
        nodes = list(G.nodes())
        for walk_iter in range(num_walks):
            t1 = time.time()
            random.shuffle(nodes)
            for node in nodes:
                walks.append(self.deepwalk_walk(walk_length=walk_length, start_node=node))
            t2 = time.time()
            print(f'Walk iteration: {walk_iter+1}/{num_walks}; time cost: {(t2-t1):.2f}')
        return walks

# ===========================================node2vec-walker============================================
class Walker:
    def __init__(self, g, p, q, workers):
        self.g = g
        self.p = p
        self.q = q

        if self.g.get_isweighted():
            # print('is weighted graph: ', self.g.get_isweighted())
            pass
        else:  # otherwise, add equal weights 1.0 to all existing edges
            # print('is weighted graph: ', self.g.get_isweighted())
            self.g.add_edge_weight(equal_weight=1.0)  # add 'weight' to networkx graph

        self.node_size = g.get_num_nodes()
        self.look_up_dict = g.look_up_dict

    def node2vec_walk(self, walk_length, start_node):
        '''
        Simulate a random walk starting from start node.
        '''
        G = self.g.G
        alias_nodes = self.alias_nodes
        alias_edges = self.alias_edges

        walk = [start_node]

        while len(walk) < walk_length:
            cur = walk[-1]
            cur_nbrs = list(G.neighbors(cur))
            if len(cur_nbrs) > 0:
                if len(walk) == 1:
                    walk.append(cur_nbrs[alias_draw(alias_nodes[cur][0], alias_nodes[cur][1])])
                else:
                    prev = walk[-2]
                    next = cur_nbrs[alias_draw(alias_edges[(prev, cur)][0], alias_edges[(prev, cur)][1])]
                    walk.append(next)
            else:
                break
        return walk

    def simulate_walks(self, num_walks, walk_length):
        '''
        Repeatedly simulate random walks from each node.
        '''
        G = self.g.G
        walks = []
        nodes = list(G.nodes())
        for walk_iter in range(num_walks):
            t1 = time.time()
            random.shuffle(nodes)
            for node in nodes:
                walks.append(self.node2vec_walk(walk_length=walk_length, start_node=node))
            t2 = time.time()
            print(f'Walk iteration: {walk_iter+1}/{num_walks}; time cost: {(t2-t1):.2f}')
        return walks

    def get_alias_edge(self, src, dst):
        '''
        Get the alias edge setup lists for a given edge.
        '''
        G = self.g.G
        p = self.p
        q = self.q

        unnormalized_probs = []
        for dst_nbr in G.neighbors(dst):
            if dst_nbr == src:
                unnormalized_probs.append(G[dst][dst_nbr]['weight']/p)
            elif G.has_edge(dst_nbr, src):
                unnormalized_probs.append(G[dst][dst_nbr]['weight'])
            else:
                unnormalized_probs.append(G[dst][dst_nbr]['weight']/q)
        norm_const = sum(unnormalized_probs)
        normalized_probs = [float(u_prob)/norm_const for u_prob in unnormalized_probs]
        return alias_setup(normalized_probs)

    def preprocess_transition_probs(self):
        '''
        Preprocessing of transition probabilities for guiding the random walks.
        '''
        G = self.g.G
        alias_nodes = {}
        for node in G.nodes():
            unnormalized_probs = [G[node][nbr]['weight'] for nbr in G.neighbors(node)]  # pick prob of neighbors with non-zero weight
            norm_const = sum(unnormalized_probs)
            normalized_probs = [float(u_prob)/norm_const for u_prob in unnormalized_probs]
            alias_nodes[node] = alias_setup(normalized_probs)

        alias_edges = {}

        if self.g.get_isdirected():
            for edge in G.edges():
                alias_edges[edge] = self.get_alias_edge(edge[0], edge[1])
        else:  # if undirected, duplicate the reverse direction; otherwise may get key error
            for edge in G.edges():
                alias_edges[edge] = self.get_alias_edge(edge[0], edge[1])
                alias_edges[(edge[1], edge[0])] = self.get_alias_edge(edge[1], edge[0])

        self.alias_nodes = alias_nodes
        self.alias_edges = alias_edges

# ========================================= utils: alias sampling method ====================================================
def alias_setup(probs):
    '''
    Compute utility lists for non-uniform sampling from discrete distributions.
    Refer to https://hips.seas.harvard.edu/blog/2013/03/03/the-alias-method-efficient-sampling-with-many-discrete-outcomes/
    for details
    '''
    K = len(probs)
    q = np.zeros(K, dtype=np.float32)
    J = np.zeros(K, dtype=np.int32)

    smaller = []
    larger = []
    for kk, prob in enumerate(probs):
        q[kk] = K*prob
        if q[kk] < 1.0:
            smaller.append(kk)
        else:
            larger.append(kk)

    while len(smaller) > 0 and len(larger) > 0:  # it is all about use large prob to compensate small prob untill reach the average
        small = smaller.pop()
        large = larger.pop()

        J[small] = large
        q[large] = q[large] + q[small] - 1.0
        if q[large] < 1.0:
            smaller.append(large)
        else:
            larger.append(large)

    return J, q  # the values in J are indexes; it is possible to have repeated indexes if that that index have large prob to compensate others

def alias_draw(J, q):
    '''
    Draw sample from a non-uniform discrete distribution using alias sampling.
    '''
    K = len(J)

    kk = int(np.floor(np.random.rand()*K))  # randomly choose a nbr (an index)
    if np.random.rand() < q[kk]:  # use alias table to choose
        return kk  # either that nbr node (an index)
    else:
        return J[kk]  # or the nbr's alias node (an index)

max_stars_repo_name: houchengbin/OpenANE

max_issues_repo_name: houchengbin/OpenANE

max_forks_repo_name: houchengbin/OpenANE

reason: import numpy,from scipy,from networkx


text:
import numpy as np
from scipy.fftpack import rfft, irfft, rfftfreq
from ....routines import rescale

def fourier_filter(data: np.ndarray, fs: float,
                   lp_freq: float = None, hp_freq: float = None, bs_freqs: list = [],
                   trans_width: float = 1, band_width: float = 1) -> np.ndarray:

    """
    Fourer filter along last axis of ``data`` with lowpass, highpass and bandstop options.

    Parameters
    ----------
    ``data`` : np.ndarray
    ``fs``: float
        sampling frequency
    ``lp_freq``: float, optional
        lowpass frequency (default is None)
    ``hp_freq``: float, optional
        highpass frequency (default is None)
    ``bs_freqs``: list, optional
        bandstop frequencies (default is [])
    ``trans_width``: float, optional
        width of the transition region between bands (default is 1)
        in frequency units
    ``band_width``: float, optional
        width of the band to remove (default is 1)
        in frequency units

    Returns
    -------
    np.ndarray
        filtered ``data``
    """

    T = data.shape[-1]
    d = 1. / fs
    freq = rfftfreq(T, d)
    f_data = rfft(data, axis=-1)

    freq_resp = create_freq_resp(freq, lp_freq,
                                 hp_freq, bs_freqs,
                                 trans_width, band_width)

    f_data = np.apply_along_axis(lambda x: x * freq_resp, -1, f_data)
    data_filtered = irfft(f_data)

    return data_filtered

def create_freq_resp(freq: np.ndarray, lp_freq: float,
                     hp_freq: float, bs_freqs: list,
                     trans_width: float, band_width: float) -> np.ndarray:
    """Calculates frequency responce for given ``freq``

    Parameters
    ----------
    ``freq``: np.ndarray, shape=(N)
        frequency array
    ``lp_freq``: float
        lowpass frequency
    ``hp_freq``: float
        highpass frequency
    ``bs_freqs``: list
        bandstop frequencies
    ``trans_width``: float
        width of the transition region between bands
    ``band_width``: float
        width of the band to remove

    Returns
    --------
    np.ndarray, shape=(N)
        frequency responce array in range form 0 to 1
    """

    freq_resp = np.ones_like(freq)

    if lp_freq is not None:
        freq_resp *= FR_lowpass(freq, lp_freq, trans_width)

    if hp_freq is not None:
        freq_resp *= FR_highpass(freq, hp_freq, trans_width)

    for bs_freq in bs_freqs:
        freq_resp *= FR_bandstop(freq, bs_freq, trans_width, band_width)

    return freq_resp

def FR_lowpass(freq: np.ndarray, lp_freq: float,
               trans_width: float) -> np.ndarray:
    """Frequency responce for lowpass filter

    Parameters
    ----------
    ``freq``: np.ndarray
        frequency array
    ``lp_freq``: float
        lowpass frequency
    ``trans_width``: float
        width of the transition region between bands

    Returns
    -------
    np.ndarray
        with values in [0, 1]
    """

    sigma = trans_width / 6.
    return 1 / (1 + np.exp((freq - lp_freq) / sigma))

def FR_highpass(freq: np.ndarray, hp_freq: float,
                trans_width: float) -> np.ndarray:
    """Frequency responce for highpass filter

    Parameters
    ----------
    ``freq``: np.ndarray
        frequency array
    ``hp_freq``: float
        highpass frequency
    ``trans_width``: float
        width of the transition region between bands

    Returns
    -------
    np.ndarray
        with values in [0, 1]
    """

    sigma = trans_width / 6.
    return 1 / (1 + np.exp((hp_freq - freq) / sigma))

def FR_bandstop(freq: np.ndarray, bs_freq: float,
                trans_width: float, band_width: float) -> np.ndarray:
    """Frequency responce for bandstop filter

    Parameters
    ----------
    ``freq``: np.ndarray
        frequency array
    ``bs_freq``: float
        bandstop frequency
    ``trans_width``: float
        width of the transition region between bands

    Returns
    -------
    np.ndarray
        with values in [0, 1]
    """

    left = FR_lowpass(freq, bs_freq - band_width / 2., trans_width)
    right = FR_highpass(freq, bs_freq + band_width / 2., trans_width)
    return rescale(left + right)

max_stars_repo_name: humanphysiologylab/pyCardiac

max_issues_repo_name: humanphysiologylab/pyCardiac

max_forks_repo_name: humanphysiologylab/pyCardiac


text:
import unittest

import numpy as np
import numpy.ma as ma

import cv2
from PIL import Image

import scipy.io as scio

#######################################
#######################################

import sys
sys.path.append('../../..')

#######################################
#######################################

from affpose.YCB_Aff import cfg as config
from affpose.YCB_Aff.dataset import dataloader
from affpose.YCB_Aff.dataset import ycb_aff_dataset_utils

#######################################
#######################################

class TestARLAffPoseDataloader(unittest.TestCase):

    def __init__(self, *args, **kwargs):
        super(TestARLAffPoseDataloader, self).__init__(*args, **kwargs)
        # load real images.
        self.dataloader = dataloader.YCBAff(split='train',
                                            select_random_images=True,
                                            num_images=100)

    def load_images(self):

        for image_idx, image_addr in enumerate(self.dataloader.img_files):
            data = self.dataloader.get_item(image_idx)

            #####################
            #####################

            rgb = data["rgb"]
            obj_label = data["obj_label"]
            aff_label = data["aff_label"]
            depth_8bit = data["depth_8bit"]

            #####################
            #####################

            colour_obj_label = ycb_aff_dataset_utils.colorize_obj_mask(obj_label)
            colour_obj_label = cv2.addWeighted(rgb, 0.35, colour_obj_label, 0.65, 0)

            colour_aff_label = ycb_aff_dataset_utils.colorize_aff_mask(aff_label)
            colour_aff_label = cv2.addWeighted(rgb, 0.35, colour_aff_label, 0.65, 0)

            #####################
            # PLOTTING
            #####################

            cv2.imshow('colour_obj_label', cv2.cvtColor(colour_obj_label, cv2.COLOR_BGR2RGB))
            cv2.imshow('colour_obj_part_label', cv2.cvtColor(colour_aff_label, cv2.COLOR_BGR2RGB))
            cv2.imshow('depth', depth_8bit)

            cv2.waitKey(0)

    def load_gt_pose(self):

        for image_idx, image_addr in enumerate(self.dataloader.img_files):
            data = self.dataloader.draw_gt_obj_pose(image_idx, project_mesh_on_image=False, verbose=True)

            #####################
            #####################

            depth_8bit = data["depth_8bit"]
            cv2_obj_pose_img = data["cv2_obj_pose_img"]
            cv2_obj_part_pose_img = data["cv2_obj_part_pose_img"]

            #####################
            # PLOTTING
            #####################

            cv2.imshow('depth', depth_8bit)
            cv2.imshow('cv2_obj_pose_img', cv2.cvtColor(cv2_obj_pose_img, cv2.COLOR_BGR2RGB))
            cv2.imshow('cv2_obj_part_pose_img', cv2.cvtColor(cv2_obj_part_pose_img, cv2.COLOR_BGR2RGB))

            cv2.waitKey(0)

    def get_occlusion_metrics(self):

        total_obj_occlusion = np.zeros(shape=(len(self.dataloader.img_files), config.NUM_OBJECTS))
        total_obj_part_occlusion = np.zeros(shape=(len(self.dataloader.img_files), config.NUM_OBJECTS))

        for image_idx, image_addr in enumerate(self.dataloader.img_files):
            data = self.dataloader.draw_gt_obj_pose(image_idx, get_occlusion_metrics=True)

            #####################
            #####################

            obj_part_label = data["obj_part_label"]
            obj_part_occlusion_mask = data["obj_part_occlusion_mask"]

            obj_occlusion_dict = {obj_id: 0 for obj_id in range(1, config.NUM_OBJECTS + 1)}
            obj_part_occlusion_dict = {obj_part_id_1: 0 for obj_part_id_1 in range(1, config.NUM_OBJECTS_PARTS + 1)}

            ####################
            # occlusion from overlap
            ####################

            obj_part_ids = list(obj_part_occlusion_mask.keys())
            for obj_part_id_1 in obj_part_ids:
                if obj_part_id_1 in ycb_aff_dataset_utils.DRAW_OBJ_PART_POSE:
                    obj_id_1 = ycb_aff_dataset_utils.map_obj_part_ids_to_obj_id(obj_part_id_1)
                    obj_part_mask_1 = obj_part_occlusion_mask[obj_part_id_1]

                    for obj_part_id_2 in obj_part_ids:
                        if obj_part_id_1 != obj_part_id_2:
                            obj_id_2 = ycb_aff_dataset_utils.map_obj_part_ids_to_obj_id(obj_part_id_2)
                            obj_part_mask_2 = obj_part_occlusion_mask[obj_part_id_2]

                            #####################
                            # intersection
                            #####################

                            intersection = cv2.bitwise_and(obj_part_mask_1, obj_part_mask_2).astype(np.uint8)
                            intersection_points = np.count_nonzero(intersection)

                            if intersection_points:

                                # remove duplicate.
                                idx = int(np.where(np.array(obj_part_ids) == obj_part_id_2)[0])
                                obj_part_ids.pop(idx)

                                # get obj_part_id.
                                idxs = np.nonzero(intersection)
                                rows, cols = idxs[0], idxs[1]
                                mask_obj_part_ids = obj_part_label[rows, cols]
                                counts = np.bincount(mask_obj_part_ids)
                                mask_obj_part_id = np.argmax(counts)

                                # get occlusion
                                if mask_obj_part_id == obj_part_id_2:
                                    occlusion = intersection_points / np.count_nonzero(obj_part_mask_1)
                                    obj_part_occlusion_dict[obj_part_id_1] += occlusion
                                    # add to object masks.
                                    # expected_obj_mask[obj_id_1] += np.count_nonzero(obj_part_mask_1)
                                    # actual_obj_mask[obj_id_1] += intersection_points
                                    obj_occlusion_dict[obj_id_1] += occlusion

                                elif mask_obj_part_id == obj_part_id_1:
                                    occlusion = intersection_points / np.count_nonzero(obj_part_mask_2)
                                    obj_part_occlusion_dict[obj_part_id_2] += occlusion
                                    # add to object masks.
                                    # expected_obj_mask[obj_id_2] += np.count_nonzero(obj_part_mask_2)
                                    obj_occlusion_dict[obj_id_2] += occlusion

            # for obj_id_1, expected_points in expected_obj_mask.items():
            #     intersection_points = actual_obj_mask[obj_id_1]
            #     if intersection_points > 0:
            #         occlusion = intersection_points / expected_points
            #         obj_occlusion_dict[obj_id_1] += occlusion

            #####################
            #####################

            for obj_id, obj_occlusion in obj_occlusion_dict.items():
                if obj_occlusion > 0:
                    obj_name = "{}".format(ycb_aff_dataset_utils.map_obj_id_to_name(obj_id))
                    total_obj_occlusion[image_idx, obj_id - 1] = obj_occlusion
                    print(f'{obj_name}: {obj_occlusion * 100:.2f} [%]')
                    obj_part_ids = ycb_aff_dataset_utils.map_obj_ids_to_obj_part_ids(obj_id)
                    for obj_part_id in obj_part_ids:
                        if obj_part_id in ycb_aff_dataset_utils.DRAW_OBJ_PART_POSE:
                            aff_id = ycb_aff_dataset_utils.map_obj_part_ids_to_aff_ids(obj_part_id)
                            aff_name = "{}".format(ycb_aff_dataset_utils.map_aff_id_to_name(aff_id))
                            obj_part_occlusion = obj_part_occlusion_dict[obj_part_id]
                            if obj_part_occlusion > 0:
                                total_obj_part_occlusion[image_idx, obj_id - 1] = obj_part_occlusion
                                print(f'\t{aff_name}: {obj_part_occlusion * 100:.2f} [%]')

            cv2.imshow('rgb', cv2.cvtColor(data["rgb"], cv2.COLOR_BGR2RGB))
            cv2.waitKey(0)

        print()
        for obj_id in range(config.NUM_OBJECTS):
            obj_name = "{:<25}".format(ycb_aff_dataset_utils.map_obj_id_to_name(obj_id + 1))
            mean_obj_occlusion = np.mean(total_obj_occlusion[:, obj_id])
            max_obj_occlusion = np.max(total_obj_occlusion[:, obj_id])
            mean_obj_part_occlusion = np.mean(total_obj_part_occlusion[:, obj_id])
            max_obj_part_occlusion = np.max(total_obj_part_occlusion[:, obj_id])
            obj_part_ids = ycb_aff_dataset_utils.map_obj_ids_to_obj_part_ids(obj_id + 1)
            for obj_part_id in obj_part_ids:
                if obj_part_id in ycb_aff_dataset_utils.DRAW_OBJ_PART_POSE:
                    aff_id = ycb_aff_dataset_utils.map_obj_part_ids_to_aff_ids(obj_part_id)
                    aff_name = "{:<15}".format(ycb_aff_dataset_utils.map_aff_id_to_name(aff_id))
            print(f'{obj_name}',
                  f'mean:{mean_obj_occlusion * 100:.2f} [%], max:{max_obj_occlusion * 100:.2f} [%]\t\t'
                  f'{aff_name}'
                  f'mean:{mean_obj_part_occlusion * 100:.2f} [%], max:{max_obj_part_occlusion * 100:.2f} [%]',
                  )

if __name__ == '__main__':
    # run all test.
    # unittest.main()

    # run desired test.
    suite = unittest.TestSuite()
    suite.addTest(TestARLAffPoseDataloader("load_gt_pose"))
    runner = unittest.TextTestRunner()
    runner.run(suite)

max_stars_repo_name: UW-Advanced-Robotics-Lab/densefusion

max_issues_repo_name: akeaveny/DenseFusion

max_forks_repo_name: akeaveny/DenseFusion

reason: import numpy,import scipy