function to plot two estimates of Bcr against each other
352802
def plot_bcr(fignum, Bcr1, Bcr2):
    
    plt.figure(num=fignum)
    plt.plot(Bcr1, Bcr2, 'ro')
    plt.xlabel('Bcr1')
    plt.ylabel('Bcr2')
    plt.title('Compare coercivity of remanence')
600915
function brcomp(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      if (Math.min(a, b) < 8) {
         return a * Math.log(x) + b * Math.log(y) - lbeta(a, b);
      }
      return 0.5 * (Math.log(a) + Math.log(b) - C.log2pi - Math.log(a + b)) -
             a * phi(1 - lambda / a) - b * phi(1 + lambda / b) -
             bcorr(a, b);
   }
456732
def function(vec1, vec2)
		case @func
		when :linear
			return normalized_radius(vec1, vec2)
		when :cubic_alt
			return normalized_radius(vec1, vec2)**(1.5)
		when :thin_plate_splines
			return 0.0 if radius(vec1, vec2) == 0.0
			return normalized_radius(vec1, vec2)**2.0 * Math.log(normalized_radius(vec1, vec2))
		when :thin_plate_splines_alt
			rnorm = ((@r0.prod.abs)**(2.0/@r0.size)*(((vec1-vec2).square / @r0.square).sum + 1.0))
			return rnorm * Math.log(rnorm)
		when :multiquadratic
# 			return Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
			(@r0.prod.abs)**(1.0/@r0.size)*Math.sqrt(((vec1-vec2).square / @r0.square).sum + 1.0)
		when :inverse_multiquadratic
			1.0 / ((@r0.prod.abs)**(1.0/@r0.size)*Math.sqrt(((vec1-vec2).square / @r0.square).sum + 1.0))
		when :cubic
			((@r0.prod.abs)**(2.0/@r0.size)*(((vec1-vec2).square / @r0.square).sum + 1.0))**(1.5)
# 			invs = ((vec1-vec2).square + @r0.square).sqrt**(-1)
# 			invs.sum
# 				p @ro
# 			return 1.0 / Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
# 		when :inverse_multiquadratic
# # 				p @ro
# 			return 1.0 / Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
		else
			raise ArgumentError.new("Bad radial basis function: #{@func}")
		end
	end
899487
function(a, b)
{
    return (a.bb_l <= b.bb_r && b.bb_l <= a.bb_r && a.bb_b <= b.bb_t && b.bb_b <= a.bb_t);
}
627123
def _interp_function(self, y_ip1, y_i, t_ip1, t_i, imt_per):
        
        return y_i + (y_ip1 - y_i) / (t_ip1 - t_i) * (imt_per - t_i)
104461
def bukin(self, x):
        """"""
        s = 0
        for k in xrange((1+len(x)) // 2):
            z = x[2 * k]
            y = x[min((2*k + 1, len(x)-1))]
            s += 100 * np.abs(y - 0.01 * z**2)**0.5 + 0.01 * np.abs(z + 10)
        return s
1009060
def _function(self, x, a, b, c, d):
        
        return c/((a**2-x**2)**2 + (b*x)**2)+d
1143979
def _to_bfloat16_unbiased(x, noise):
  
  x_sign = tf.sign(x)
  # Make sure x is positive.  If it is zero, the two candidates are identical.
  x = x * x_sign + 1e-30
  cand1 = tf.to_bfloat16(x)
  cand1_f = tf.to_float(cand1)
  # This relies on the fact that for a positive bfloat16 b,
  # b * 1.005 gives you the next higher bfloat16 and b*0.995 gives you the
  # next lower one. Both 1.005 and 0.995 are ballpark estimation.
  cand2 = tf.to_bfloat16(
      tf.where(tf.greater(x, cand1_f), cand1_f * 1.005, cand1_f * 0.995))
  ret = _randomized_roundoff_to_bfloat16(x, noise, cand1, cand2)
  return ret * tf.to_bfloat16(x_sign)
184722
def fun_residuals(params, xnor, ynor, w, bbox, k, ext):
    """"""

    spl = LSQUnivariateSpline(
        x=xnor,
        y=ynor,
        t=[item.value for item in params.values()],
        w=w,
        bbox=bbox,
        k=k,
        ext=ext,
        check_finite=False
    )
    return spl.get_residual()
1002326
def plot_xbt(fignum, XB, T, e, b):
    
    plt.figure(num=fignum)
    plt.clf()
    if not isServer:
        plt.figtext(.02, .01, version_num)
    plt.xlabel('Field (T)')
    plt.ylabel('Susceptibility (m^3/kg)')
    k = 0
    B, X = [], []
    for xb in XB:
        X.append(xb[0])
        B.append(xb[1])
    plt.plot(B, X)
    plt.legend(['%i' % (int(T)) + ' K'])
    plt.title(e + ': f = ' + '%i' % (int(f)) + ' Hz')
997185
double psfunc2(double x, double t1, int dim, int type, int css) {
    // type:  the type of curve 
    // css:   the class of index
    double beta;
    beta = 0.0;

    dim++;

    if (type == 21) {
      double xy = 2 * (x - 0.5);
      beta = xy - Math.pow(t1, 0.5 * (nvar + 3 * dim - 8) / (nvar - 2));
    }

    if (type == 22) {
      double theta = 6 * Math.PI * t1 + dim * Math.PI / nvar;
      double xy = 2 * (x - 0.5);
      beta = xy - Math.sin(theta);
    }

    if (type == 23) {
      double theta = 6 * Math.PI * t1 + dim * Math.PI / nvar;
      double ra = 0.8 * t1;
      double xy = 2 * (x - 0.5);
      if (css == 1) {
        beta = xy - ra * Math.cos(theta);
      } else {
        beta = xy - ra * Math.sin(theta);
      }
    }

    if (type == 24) {
      double theta = 6 * Math.PI * t1 + dim * Math.PI / nvar;
      double xy = 2 * (x - 0.5);
      double ra = 0.8 * t1;
      if (css == 1) {
        beta = xy - ra * Math.cos(theta / 3);
      } else {
        beta = xy - ra * Math.sin(theta);
      }
    }

    if (type == 25) {
      double rho = 0.8;
      double phi = Math.PI * t1;
      double theta = 6 * Math.PI * t1 + dim * Math.PI / nvar;
      double xy = 2 * (x - 0.5);
      if (css == 1) {
        beta = xy - rho * Math.sin(phi) * Math.sin(theta);
      } else if (css == 2) {
        beta = xy - rho * Math.sin(phi) * Math.cos(theta);
      } else {
        beta = xy - rho * Math.cos(phi);
      }
    }

    if (type == 26) {
      double theta = 6 * Math.PI * t1 + dim * Math.PI / nvar;
      double ra = 0.3 * t1 * (t1 * Math.cos(4 * theta) + 2);
      double xy = 2 * (x - 0.5);
      if (css == 1) {
        beta = xy - ra * Math.cos(theta);
      } else {
        beta = xy - ra * Math.sin(theta);
      }
    }

    return beta;
  }
252038
def py_bisect(f, a, b, args=(), xtol=_xtol, rtol=_rtol, maxiter=_iter,
              ytol=None, full_output=False, disp=True):
    '''
    '''
    fa = f(a, *args)
    fb = f(b, *args)
    if fa*fb > 0.0:
        raise ValueError("f(a) and f(b) must have different signs") 
    elif fa == 0.0:
        return a
    elif fb == 0.0:
        return b

    dm = b - a
    iterations = 0.0
    
    for i in range(maxiter):
        dm *= 0.5
        xm = a + dm
        fm = f(xm, *args)
        if fm*fa >= 0.0:
            a = xm
        abs_dm = fabs(dm)
        
        if fm == 0.0:
            return xm
        elif ytol is not None:
            if (abs_dm < xtol + rtol*abs_dm) and abs(fm) < ytol:
                return xm
        elif (abs_dm < xtol + rtol*abs_dm):
            return xm
    raise UnconvergedError("Failed to converge after %d iterations" %maxiter)
1202362
function cubicBezierAtTime(t,p1x,p1y,p2x,p2y,duration) {
      var ax = 0,bx = 0,cx = 0,ay = 0,by = 0,cy = 0;
      function sampleCurveX(t) {return ((ax * t + bx) * t + cx) * t;}
      function sampleCurveY(t) {return ((ay * t + by) * t + cy) * t;}
      function sampleCurveDerivativeX(t) {return (3.0 * ax * t + 2.0 * bx) * t + cx;}
      function solveEpsilon(duration) {return 1.0 / (200.0 * duration);}
      function solve(x,epsilon) {return sampleCurveY(solveCurveX(x,epsilon));}
      function fabs(n) {if (n >= 0) {return n;}else {return 0 - n;}}
      function solveCurveX(x,epsilon) {
        var t0,t1,t2,x2,d2,i;
        for (t2 = x, i = 0; i < 8; i++) {x2 = sampleCurveX(t2) - x; if (fabs(x2) < epsilon) {return t2;} d2 = sampleCurveDerivativeX(t2); if (fabs(d2) < 1e-6) {break;} t2 = t2 - x2 / d2;}
        t0 = 0.0; t1 = 1.0; t2 = x; if (t2 < t0) {return t0;} if (t2 > t1) {return t1;}
        while (t0 < t1) {x2 = sampleCurveX(t2); if (fabs(x2 - x) < epsilon) {return t2;} if (x > x2) {t0 = t2;}else {t1 = t2;} t2 = (t1 - t0) * 0.5 + t0;}
        return t2; // Failure.
      }
      cx = 3.0 * p1x; bx = 3.0 * (p2x - p1x) - cx; ax = 1.0 - cx - bx; cy = 3.0 * p1y; by = 3.0 * (p2y - p1y) - cy; ay = 1.0 - cy - by;
      return solve(t, solveEpsilon(duration));
    }
826744
function(a, b){
    return new BB(
            min(a.l, b.l),
            min(a.b, b.b),
            max(a.r, b.r),
            max(a.t, b.t)
        );
}
884310
def _twofilter_smoothing_ON(self, t, ti, info, phi, lwinfo, return_ess,
                               modif_forward, modif_info):
        
        if modif_info is not None:
            lwinfo += modif_info
        Winfo = rs.exp_and_normalise(lwinfo)
        I = rs.multinomial(Winfo)
        if modif_forward is not None:
            lw = self.wgt[t].lw + modif_forward
            W = rs.exp_and_normalise(lw)
        else:
            W = self.wgt[t].W
        J = rs.multinomial(W)
        log_omega = self.model.logpt(t + 1, self.X[t][J], info.hist.X[ti][I])
        if modif_forward is not None:
            log_omega -= modif_forward[J]
        if modif_info is not None:
            log_omega -= modif_info[I]
        Om = rs.exp_and_normalise(log_omega)
        est = np.average(phi(self.X[t][J], info.hist.X[ti][I]), axis=0,
                         weights=Om)
        if return_ess:
            return (est, 1. / np.sum(Om**2))
        else:
            return est
40423
function invIncBeta(p, a, b) {
  if(x < 0 || x > 1) {
    throw new RangeError('First argument must be between 0 and 1.');
	}

  // Special cases, there can make trouble otherwise
  else if (a === 1 && b === 1) return p;
  else if (p === 1) return 1;
  else if (p === 0) return 0;
  else if (a === 0) return 0;
  else if (b === 0) return 1;

  else {
    var EPS = 1e-8,
        a1 = a - 1,
        b1 = b - 1,
        j = 0,
        lna, lnb, pp, t, u, err, x, al, h, w, afac;

	if(a >= 1 && b >= 1) {
    pp = (p < 0.5) ? p : 1 - p;
    t = Math.sqrt(-2 * Math.log(pp));

		x = (2.30753 + t * 0.27061) / (1 + t* (0.99229 + t * 0.04481)) - t;
		if(p < 0.5) x = -x;
		al = (x * x - 3) / 6;
		h = 2 / (1 / (2 * a - 1)  + 1 / (2 * b - 1));
		w = (x * Math.sqrt(al + h) / h) - (1 / (2 * b - 1) - 1 / (2 * a - 1)) * (al + 5 / 6 - 2 / (3 * h));
		x = a / (a + b * Math.exp(2 * w));
	} else {
		lna = Math.log(a / (a + b));
		lnb = Math.log(b / (a + b));
		t = Math.exp(a * lna) / a;
		u = Math.exp(b * lnb) / b;
		w = t + u;
		if (p < t / w) x = Math.pow(a * w * p, 1 / a);
		else x = 1 - Math.pow(b * w * (1 - p), 1 / b);
	}

	afac = -gammaCollection.logGamma(a) - gammaCollection.logGamma(b) + gammaCollection.logGamma(a + b);

  for(; j < 10; j++) {
		if(x === 0 || x === 1) return x;
		err = regualizedBeta(x, a, b) - p;

    t = Math.exp(a1 * Math.log(x) + b1 * log1p(-x) + afac);
		u = err / t;
		x -= (t = u / (1 - 0.5 * Math.min(1, u * (a1 / x - b1 / (1 - x)))));

    if (x <= 0) x = 0.5 * (x + t);
		if (x >= 1) x = 0.5 * (x + t + 1);

		if (Math.abs(t) < EPS * x && j > 0) break;
	}

	return x;
  }
}
1285875
def plot_multi(func:Callable[[int,int,plt.Axes],None], r:int=1, c:int=1, figsize:Tuple=(12,6)):
    ""
    axes = plt.subplots(r, c, figsize=figsize)[1]
    for i in range(r):
        for j in range(c): func(i,j,axes[i,j])