/*
auxiliary function for approxKdist methods.
510011
private double ssqerr(int k0, int kmax, double[] logk, double[] log_kDist, double m, double t) {
    int k = kmax - k0;
    double result = 0;
    for(int i = 0; i < k; i++) {
      // double h = log_kDist[i] - (m * (logk[i] - logk[0]) + t); ???
      double h = log_kDist[i] - m * logk[i] - t;
      result += h * h;
    }
    return result;
  }
882826
def _call(self, x):
        """"""
        if isinstance(self.domain, ProductSpace):
            norm = PointwiseNorm(self.domain, 2)(x)
        else:
            norm = x.ufuncs.absolute()

        if self.gamma > 0:
            tmp = norm.ufuncs.square()
            tmp *= 1 / (2 * self.gamma)

            index = norm.ufuncs.greater_equal(self.gamma)
            tmp[index] = norm[index] - self.gamma / 2
        else:
            tmp = norm

        return tmp.inner(tmp.space.one())
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
1102730
def _D(self, R):
        """"""
        return (0.3 * self._HNn**2 / self._sin_alpha / R
                + self._HNn + R * self._sin_alpha) / (0.3 * self._HNn + R * self._sin_alpha)
1028992
def numbaGaussian2d(psf, sy, sx):
    '''
    
    '''
    ps0, ps1 = psf.shape
    c0,c1 = ps0//2, ps1//2
    ssx = 2*sx**2
    ssy = 2*sy**2
    for i in range(ps0):
        for j in range(ps1):
            psf[i,j]=exp( -( (i-c0)**2/ssy
                            +(j-c1)**2/ssx) )
    psf/=psf.sum()
975887
def interpolate(x, y, z, interp_type='linear', hres=50000,
                minimum_neighbors=3, gamma=0.25, kappa_star=5.052,
                search_radius=None, rbf_func='linear', rbf_smooth=0,
                boundary_coords=None):
    """"""
    return interpolate_to_grid(x, y, z, interp_type=interp_type, hres=hres,
                               minimum_neighbors=minimum_neighbors, gamma=gamma,
                               kappa_star=kappa_star, search_radius=search_radius,
                               rbf_func=rbf_func, rbf_smooth=rbf_smooth,
                               boundary_coords=boundary_coords)
710602
def mape(y, p):
    

    filt = np.abs(y) > EPS
    return np.mean(np.abs(1 - p[filt] / y[filt]))
687079
def _dK_dR(self, R):
        """"""
        return -self._ns * self._N / R**2 / self._sin_alpha
1063892
protected ApproximationLine progressiveKnnDistanceApproximation(int k_max) {
    if(!isLeaf()) {
      throw new UnsupportedOperationException("Progressive KNN-distance approximation " + "is only vailable in leaf nodes!");
    }

    // determine k_0, y_1, y_kmax
    int k_0 = 0;
    double y_1 = Double.POSITIVE_INFINITY;
    double y_kmax = Double.POSITIVE_INFINITY;

    for(int i = 0; i < getNumEntries(); i++) {
      MkCoPLeafEntry entry = (MkCoPLeafEntry) getEntry(i);
      ApproximationLine approx = entry.getProgressiveKnnDistanceApproximation();
      k_0 = Math.max(approx.getK_0(), k_0);
    }

    for(int i = 0; i < getNumEntries(); i++) {
      MkCoPLeafEntry entry = (MkCoPLeafEntry) getEntry(i);
      ApproximationLine approx = entry.getProgressiveKnnDistanceApproximation();
      y_1 = Math.min(approx.getValueAt(k_0), y_1);
      y_kmax = Math.min(approx.getValueAt(k_max), y_kmax);
    }

    // determine m and t
    double m = (y_kmax - y_1) / (FastMath.log(k_max) - FastMath.log(k_0));
    double t = y_1 - m * FastMath.log(k_0);

    return new ApproximationLine(k_0, m, t);
  }
365365
def approximate_column(self, i):
        

        """
        return np.dot(self._C_k, self._R_k[:, i])
1351087
function findClosestPointIndexWithinRadius(points, x, y, radiusThreshold) {
    let closestIdx;
    let closestDist = Number.MAX_VALUE;
    points.forEach((d, idx) => {
      let distSqrd = Math.pow(d.x - x, 2) + Math.pow(d.y - y, 2); // changeto: (d.x - x)**2 + (d.y - y)**2;
      if (distSqrd < closestDist && distSqrd < radiusThreshold) {
        closestIdx = idx;
        closestDist = distSqrd;
      }
    });
    return closestIdx;
  }
924908
def _ndtri(y):
    
    # approximation for 0 <= abs(z - 0.5) <= 3/8
    P0 = [
        -5.99633501014107895267E1,
        9.80010754185999661536E1,
        -5.66762857469070293439E1,
        1.39312609387279679503E1,
        -1.23916583867381258016E0,
    ]

    Q0 = [
        1.95448858338141759834E0,
        4.67627912898881538453E0,
        8.63602421390890590575E1,
        -2.25462687854119370527E2,
        2.00260212380060660359E2,
        -8.20372256168333339912E1,
        1.59056225126211695515E1,
        -1.18331621121330003142E0,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 2 and 8
    # i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    P1 = [
        4.05544892305962419923E0,
        3.15251094599893866154E1,
        5.71628192246421288162E1,
        4.40805073893200834700E1,
        1.46849561928858024014E1,
        2.18663306850790267539E0,
        -1.40256079171354495875E-1,
        -3.50424626827848203418E-2,
        -8.57456785154685413611E-4,
    ]

    Q1 = [
        1.57799883256466749731E1,
        4.53907635128879210584E1,
        4.13172038254672030440E1,
        1.50425385692907503408E1,
        2.50464946208309415979E0,
        -1.42182922854787788574E-1,
        -3.80806407691578277194E-2,
        -9.33259480895457427372E-4,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 8 and 64
    # i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    P2 = [
        3.23774891776946035970E0,
        6.91522889068984211695E0,
        3.93881025292474443415E0,
        1.33303460815807542389E0,
        2.01485389549179081538E-1,
        1.23716634817820021358E-2,
        3.01581553508235416007E-4,
        2.65806974686737550832E-6,
        6.23974539184983293730E-9,
    ]

    Q2 = [
        6.02427039364742014255E0,
        3.67983563856160859403E0,
        1.37702099489081330271E0,
        2.16236993594496635890E-1,
        1.34204006088543189037E-2,
        3.28014464682127739104E-4,
        2.89247864745380683936E-6,
        6.79019408009981274425E-9,
    ]

    sign_flag = 1

    if y > (1 - EXP_NEG2):
        y = 1 - y
        sign_flag = 0

    # Shortcut case where we don't need high precision
    # between -0.135 and 0.135
    if y > EXP_NEG2:
        y -= 0.5
        y2 = y ** 2
        x = y + y * (y2 * _polevl(y2, P0, 4) / _p1evl(y2, Q0, 8))
        x = x * ROOT_2PI
        return x

    x = math.sqrt(-2.0 * math.log(y))
    x0 = x - math.log(x) / x

    z = 1.0 / x
    if x < 8.0:                 # y > exp(-32) = 1.2664165549e-14
        x1 = z * _polevl(z, P1, 8) / _p1evl(z, Q1, 8)
    else:
        x1 = z * _polevl(z, P2, 8) / _p1evl(z, Q2, 8)

    x = x0 - x1
    if sign_flag != 0:
        x = -x

    return x
95046
def _assess_approximation(self, boundary_points, derivs, funcs, nodes, problem):
        
        interior_residuals = self._compute_interior_residuals(derivs, funcs,
                                                              nodes, problem)
        boundary_residuals = self._compute_boundary_residuals(boundary_points,
                                                              funcs, problem)
        return np.hstack(interior_residuals + boundary_residuals)
465660
def eps_func(self):
        '''
        
        '''
        interp_real = interpolate.interp2d(self.x, self.y, self.eps.real)
        interp_imag = interpolate.interp2d(self.x, self.y, self.eps.imag)
        interp = lambda x, y: interp_real(x, y) + 1.j*interp_imag(x, y)
        return interp
64916
def smooth(x, y, weights):
    '''
    
    '''
    # Spline hard to smooth properly, therefore solfed with
    # bounded polynomal interpolation
    # ext=3: no extrapolation, but boundary value
#     return UnivariateSpline(x, y, w=weights,
#                             s=len(y)*weights.max()*100, ext=3)

#     return np.poly1d(np.polyfit(x,y,w=weights,deg=2))
    p = np.polyfit(x, y, w=weights, deg=2)
    if np.any(np.isnan(p)):
        # couldn't even do polynomial fit
        # as last option: assume constant noise
        my = np.average(y, weights=weights)
        return lambda x: my
    return lambda xint: np.poly1d(p)(np.clip(xint, x[0], x[-1]))
968974
def _K(self, R):
        """"""
        return self._ns * self._N / R / self._sin_alpha
828209
function findClosestPointOnSpecialAttractor(x,y) {
  var nSpecialAttractor = specialAttractors.length;
  var currentMinDistance = 0;
  var currentSpecialAttractor;
  var ox = 0;
  var oy = 0;
  for(var a=0; a<nSpecialAttractor; a++) {
    var specialAttractor = specialAttractors[a];
    var closestSegmentPoint = distanceToSegment(specialAttractor.x1, specialAttractor.y1, specialAttractor.x2, specialAttractor.y2, x, y);
    if(a == 0 || (closestSegmentPoint.distance) < currentMinDistance) {
      currentMinDistance = closestSegmentPoint.distance;
      currentSpecialAttractor = specialAttractor;
      ox = closestSegmentPoint.originX;
      oy = closestSegmentPoint.originY;
    }
  }
  return {
    distance: currentMinDistance,
    specialAttractor: currentSpecialAttractor,
    originX: ox,
    originY: oy
  }
}