Inverse Generalised Pareto distribution function.
511391
def gpinv(p, k, sigma):
    """"""
    x = np.empty(p.shape)
    x.fill(np.nan)
    if sigma <= 0:
        return x
    ok = (p > 0) & (p < 1)
    if np.all(ok):
        if np.abs(k) < np.finfo(float).eps:
            np.negative(p, out=x)
            np.log1p(x, out=x)
            np.negative(x, out=x)
        else:
            np.negative(p, out=x)
            np.log1p(x, out=x)
            x *= -k
            np.expm1(x, out=x)
            x /= k
        x *= sigma
    else:
        if np.abs(k) < np.finfo(float).eps:
            # x[ok] = - np.log1p(-p[ok])
            temp = p[ok]
            np.negative(temp, out=temp)
            np.log1p(temp, out=temp)
            np.negative(temp, out=temp)
            x[ok] = temp
        else:
            # x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
            temp = p[ok]
            np.negative(temp, out=temp)
            np.log1p(temp, out=temp)
            temp *= -k
            np.expm1(temp, out=temp)
            temp /= k
            x[ok] = temp
        x *= sigma
        x[p == 0] = 0
        if k >= 0:
            x[p == 1] = np.inf
        else:
            x[p == 1] = -sigma / k
    return x
705367
def weighted_axioms(self, x, y, xg):
        
        # TODO: allow additional weighting
        # weights are log odds w=log(p/(1-p))
        # (Sub,Sup,Eq,Other)
        scope_pairs = [
            ('label',   'label',   0.0, 0.0, 3.0,-0.8), 
            ('label',   'exact',   0.0, 0.0, 2.5,-0.5), 
            ('label',   'broad',  -1.0, 1.0, 0.0, 0.0), 
            ('label',   'narrow',  1.0,-1.0, 0.0, 0.0), 
            ('label',   'related', 0.0, 0.0, 0.0, 0.0), 
            ('exact',   'exact',   0.0, 0.0, 2.5,-0.5), 
            ('exact',   'broad',  -1.0, 1.0, 0.0, 0.0), 
            ('exact',   'narrow',  1.0,-1.0, 0.0, 0.0), 
            ('exact',   'related', 0.0, 0.0, 0.0, 0.0), 
            ('related', 'broad',  -0.5, 0.5, 0.0, 0.0), 
            ('related', 'narrow',  0.5,-0.5, 0.0, 0.0), 
            ('related', 'related', 0.0, 0.0, 0.0, 0.0), 
            ('broad',   'broad',   0.0, 0.0, 0.0, 1.0), 
            ('broad',   'narrow', -0.5, 0.5, 0.0, 0.2), 
            ('narrow',  'narrow',  0.0, 0.0, 0.0, 0.0)
        ]
        # populate symmetric lookup matrix
        scope_map = defaultdict(dict)
        for (l,r,w1,w2,w3,w4) in scope_pairs:
            l = l.upper()
            r = r.upper()
            scope_map[l][r] = np.array((w1,w2,w3,w4))
            scope_map[r][l] = np.array((w2,w1,w3,w4))

        # TODO: get prior based on ontology pair
        # cumulative sum of weights
        WS = None
        pfx1 = self._id_to_ontology(x)
        pfx2 = self._id_to_ontology(y)
        for mw in self.config.get('match_weights', []):
            mpfx1 = mw.get('prefix1','')
            mpfx2 = mw.get('prefix2','')
            X = np.array(mw['weights'])
            if mpfx1 == pfx1 and mpfx2 == pfx2:
                WS = X
            elif mpfx2 == pfx1 and mpfx1 == pfx2:
                WS = self._flipweights(X)
            elif mpfx1 == pfx1 and mpfx2 == '' and WS is None:
                WS = X
            elif mpfx2 == pfx1 and mpfx1 == '' and WS is None:
                WS = self._flipweights(X)

        if WS is None:
            WS = np.array((0.0, 0.0, 0.0, 0.0))
        # defaults
        WS += np.array(self.config.get('default_weights', [0.0, 0.0, 1.5, -0.1]))
        logging.info('WS defaults={}'.format(WS))

        for xw in self.config.get('xref_weights', []):
            left = xw.get('left','')
            right = xw.get('right','')
            X = np.array(xw['weights'])
            if x == left and y == right:
                WS += X
                logging.info('MATCH: {} for {}-{}'.format(X, x, y))
            elif y == left and x == right:
                WS += self._flipweights(X)
                logging.info('IMATCH: {}'.format(X))

        smap = self.smap
        # TODO: symmetrical
        WT = np.array((0.0, 0.0, 0.0, 0.0))
        WBESTMAX = np.array((0.0, 0.0, 0.0, 0.0))
        n = 0
        for sx in smap[x]:
            WBEST, _ = self._best_match_syn(sx, smap[y], scope_map)
            if WBEST is not None:
                WT += WBEST
                n += 1
                if max(abs(WBEST)) > max(abs(WBESTMAX)):
                    WBESTMAX = WBEST
        for sy in smap[y]:
            WBEST, _ = self._best_match_syn(sy, smap[x], scope_map)
            if WBEST is not None:
                WT += WBEST
                n += 1
        # average best match
        if n > 0:
            logging.info('Adding BESTMAX={}'.format(WBESTMAX))
            WS += WBESTMAX
                    
        # TODO: xref, many to many
        WS += self._graph_weights(x, y, xg)
        # TODO: include additional defined weights, eg ORDO
        logging.info('Adding WS, gw={}'.format(WS))

        # jaccard similarity
        (ss1,ss2) = xg[x][y][self.SIMSCORES]
        WS[3] += ((1-ss1) + (1-ss2)) / 2

        # reciprocal best hits are higher confidence of equiv
        rs = xg[x][y]['reciprocal_score']
        if rs == 4:
            WS[2] += 0.5
        if rs == 0:
            WS[2] -= 0.2
        
        #P = np.expit(WS)
        P = 1/(1+np.exp(-WS))
        logging.info('Final WS={}, init P={}'.format(WS, P))
        # probs should sum to 1.0
        P = P / np.sum(P)
        return P
295969
private static double score1static(Chunk[] chks, int treeIdx, double offset, double[] fs, int row, Distribution dist, int nClasses) {
    double f = chks[treeIdx].atd(row) + offset;
    double p = dist.linkInv(f);
    if (dist.distribution == DistributionFamily.modified_huber || dist.distribution == DistributionFamily.bernoulli || dist.distribution == DistributionFamily.quasibinomial) {
      fs[2] = p;
      fs[1] = 1.0 - p;
      return 1;                 // f2 = 1.0 - f1; so f1+f2 = 1.0
    } else if (dist.distribution == DistributionFamily.multinomial) {
      if (nClasses == 2) {
        // This optimization assumes the 2nd tree of a 2-class system is the
        // inverse of the first.  Fill in the missing tree
        fs[1] = p;
        fs[2] = 1 / p;
        return fs[1] + fs[2];
      }
      // Multinomial loss function; sum(exp(data)).  Load tree data
      assert (offset == 0);
      fs[1] = f;
      for (int k = 1; k < nClasses; k++)
        fs[k + 1] = chks[treeIdx + k].atd(row);
      // Rescale to avoid Infinities; return sum(exp(data))
      return hex.genmodel.GenModel.log_rescale(fs);
    } else {
      return fs[0] = p;
    }
  }
1330207
def partial_fit(self, X):
        
        opt, cost = self.sess.run((self.optimizer, self.cost), 
                                  feed_dict={self.x: X})
        return cost
852694
def create(entropy_coefficient, value_coefficient, cliprange, max_grad_norm, discount_factor,
           normalize_advantage=True, gae_lambda=1.0):
    
    return PpoPolicyGradient(
        entropy_coefficient, value_coefficient, cliprange, max_grad_norm,
        discount_factor=discount_factor,
        normalize_advantage=normalize_advantage,
        gae_lambda=gae_lambda
    )
736270
def get_p_vals(self, X):
		'''
		

		'''
		z_scores = self.get_scores(X[:, 0], X[:, 1])
		return norm.cdf(z_scores)
31101
def cumulative_distribution(self, X, U=0):
        
        self.check_fit()

        low_bounds = self.model.dataset.mean() - (5 * self.model.dataset.std())
        return self.model.integrate_box_1d(low_bounds, X) - U
1354308
def _preoptimize_model(self):
        
        gaussian_model = DynReg(formula=self.formula, data=self.data_original)
        gaussian_model.fit()

        for i in range(self.z_no-self.family_z_no):
            self.latent_variables.z_list[i].start = gaussian_model.latent_variables.get_z_values()[i+1]

        if self.model_name2 == 't':

            def temp_function(params):
                return -np.sum(ss.t.logpdf(x=self.data, df=np.exp(params[0]), 
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2])))

            p = optimize.minimize(temp_function,np.array([2.0,0.0,-1.0]),method='L-BFGS-B')
            self.latent_variables.z_list[-2].start = p.x[2]
            self.latent_variables.z_list[-1].start = p.x[0]

        elif self.model_name2 == 'Skewt':

            def temp_function(params):
                return -np.sum(fam.Skewt.logpdf_internal(x=self.data,df=np.exp(params[0]),
                    loc=np.ones(self.data.shape[0])*params[1], scale=np.exp(params[2]),gamma=np.exp(params[3])))

            p = optimize.minimize(temp_function,np.array([2.0,0.0,-1.0,0.0]),method='L-BFGS-B')
            self.latent_variables.z_list[-3].start = p.x[3]
            self.latent_variables.z_list[-2].start = p.x[2]
            self.latent_variables.z_list[-1].start = p.x[0]

        return gaussian_model.latent_variables
1043978
def reeval(self, X, fit, func, ask, args=()):
        
        self.fit = list(fit)
        self.fitre = list(fit)
        self.idx = self.indices(fit)
        if not len(self.idx):
            return self.idx
        evals = int(self.evaluations) if self.f_aggregate else 1
        fagg = np.median if self.f_aggregate is None else self.f_aggregate
        for i in self.idx:
            X_i = X[i]
            if self.epsilon:
                if self.parallel:
                    self.fitre[i] = fagg(func(ask(evals, X_i, self.epsilon), *args))
                else:
                    self.fitre[i] = fagg([func(ask(1, X_i, self.epsilon)[0], *args)
                                            for _k in xrange(evals)])
            else:
                self.fitre[i] = fagg([func(X_i, *args) for _k in xrange(evals)])
        self.evaluations_just_done = evals * len(self.idx)
        return self.fit, self.fitre, self.idx
647825
def optimize_NxN(x, y, fit_offset=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = get_weight(x, y, perc).astype(bool)
        if issparse(weights): weights = weights.A
    else:
        weights = None

    x, y = x.astype(np.float64), y.astype(np.float64)

    n_vars = x.shape[1]
    offset, gamma = np.zeros(n_vars), np.zeros(n_vars)

    for i in range(n_vars):
        xi = x[:, i] if weights is None else x[:, i][weights[:, i]]
        yi = y[:, i] if weights is None else y[:, i][weights[:, i]]

        if fit_offset:
            offset[i], gamma[i] = minimize(lambda m: np.sum((-yi + xi * m[1] + m[0])**2), method="L-BFGS-B",
                                           x0=(0, 0.1), bounds=[(0, None), (None, None)]).x
        else:
            gamma[i] = minimize(lambda m: np.sum((-yi + xi * m) ** 2), x0=0.1, method="L-BFGS-B").x
    offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
    return offset, gamma
716147
def combined_loss_given_predictions(log_probab_actions_new,
                                    log_probab_actions_old,
                                    value_prediction,
                                    padded_actions,
                                    padded_rewards,
                                    reward_mask,
                                    gamma=0.99,
                                    lambda_=0.95,
                                    epsilon=0.2,
                                    c1=1.0,
                                    c2=0.01):
  """"""
  loss_value = value_loss_given_predictions(
      value_prediction, padded_rewards, reward_mask, gamma=gamma)
  loss_ppo = ppo_loss_given_predictions(log_probab_actions_new,
                                        log_probab_actions_old,
                                        value_prediction,
                                        padded_actions,
                                        padded_rewards,
                                        reward_mask,
                                        gamma=gamma,
                                        lambda_=lambda_,
                                        epsilon=epsilon)
  # TODO(afrozm): Add the entropy bonus, but since we don't do that in T2T
  # we'll skip if for now.
  entropy_bonus = 0.0
  return (loss_ppo + (c1 * loss_value) - (c2 * entropy_bonus), loss_ppo,
          loss_value, entropy_bonus)
742312
def proper(self,x):
        """"""
        x[x < 0] = self.max_fit
        x[np.isnan(x)] = self.max_fit
        x[np.isinf(x)] = self.max_fit
        return x
486078
def random_offspring(self):
        ""
        function_set = self.function_set
        function_selection = self._function_selection_ins
        function_selection.density = self.population.density
        function_selection.unfeasible_functions.clear()
        for i in range(self._number_tries_feasible_ind):
            if self._function_selection:
                func_index = function_selection.tournament()
            else:
                func_index = function_selection.random_function()
            func = function_set[func_index]
            args = self.get_args(func)
            if args is None:
                continue
            args = [self.population.population[x].position for x in args]
            f = self._random_offspring(func, args)
            if f is None:
                function_selection.unfeasible_functions.add(func_index)
                continue
            function_selection[func_index] = f.fitness
            return f
        raise RuntimeError("Could not find a suitable random offpsring")
394987
def _M(self, k, t, tips=None, gaps=None):
        """"""
        if self._distributionmodel:
            return self.model.M(k, t, tips, gaps)
        else:
            return self.model.M(t, tips, gaps)
1248388
function norminv(p) {

  function rateval(a, b, x) {
    var u = a[a.length - 1];
    for (var i = a.length - 1; i > 0; --i) {
      u = x * u + a[i - 1];
    }
    var v = b[b.length - 1];
    for (var j = b.length - 1; j > 0; --j) {
      v = x * v + b[j - 1];
    }
    r = u / v;
    return r;
  };

  var a = [
    3.3871328727963666080,     1.3314166789178437745e+2,
    1.9715909503065514427e+3,  1.3731693765509461125e+4,
    4.5921953931549871457e+4,  6.7265770927008700853e+4,
    3.3430575583588128105e+4,  2.5090809287301226727e+3];
  var b = [
    1.0,                       4.2313330701600911252e+1,
    6.8718700749205790830e+2,  5.3941960214247511077e+3,
    2.1213794301586595867e+4,  3.9307895800092710610e+4,
    2.8729085735721942674e+4,  5.2264952788528545610e+3];
  var c = [
    1.42343711074968357734,     4.63033784615654529590,
    5.76949722146069140550,     3.64784832476320460504,
    1.27045825245236838258,     2.41780725177450611770e-1,
    2.27238449892691845833e-2,  7.74545014278341407640e-4];
  var d = [
    1.0,                        2.05319162663775882187,
    1.67638483018380384940,     6.89767334985100004550e-1,
    1.48103976427480074590e-1,  1.51986665636164571966e-2,
    5.47593808499534494600e-4,  1.05075007164441684324e-9];
  var e = [
    6.65790464350110377720,     5.46378491116411436990,
    1.78482653991729133580,     2.96560571828504891230e-1,
    2.65321895265761230930e-2,  1.24266094738807843860e-3,
    2.71155556874348757815e-5,  2.01033439929228813265e-7];
  var f = [
    1.0,                        5.99832206555887937690e-1,
    1.36929880922735805310e-1,  1.48753612908506148525e-2,
    7.86869131145613259100e-4,  1.84631831751005468180e-5,
    1.42151175831644588870e-7,  2.04426310338993978564e-15];

  if (p <= 0.0) {
    return Number.NEGATIVE_INFINITY;
  }
  if (p >= 1.0) {
    return Number.POSITIVE_INFINITY;
  }

  var q = p - 0.5;
  if (Math.abs(q) <= 0.425) {
    return q * rateval(a, b, 0.180625 - q*q);
  }

  var r = Math.sqrt(-Math.log(q < 0.0 ? p : 1.0 - p));
  var x;
  if (r <= 5.0) {
    x = rateval(c, d, r - 1.6);
  } else {
    x = rateval(e, f, r - 0.5);
  }
  if (p < 0.5) {
    return -x;
  }
  return x;
}
1195902
def lab2xyz(l__, a__, b__):
    


    def f_inv(arr):
        """Inverse of the f function.
        """
        return np.where(arr > 6.0/29.0,
                        arr ** 3,
                        3 * (6.0 / 29.0) * (6.0 / 29.0) * (arr - 4.0 / 29.0))

    new_l = (l__ + 16.0) / 116.0
    x__ = 95.047 * f_inv(new_l + a__ / 500.0)
    y__ = 100.0 * f_inv(new_l)
    z__ = 108.883 * f_inv(new_l - b__ / 200.0)

    return x__, y__, z__
53289
def min_func(fun, bounds_min, bounds_max, d=None, rmax=10,
             n0=64, nswp=10, verb=True, smooth_fun=None):
    """"""
    if d is None:
        d = len(bounds_min)
        a = np.asanyarray(bounds_min).copy()
        b = np.asanyarray(bounds_max).copy()
    else:
        a = np.ones(d) * bounds_min
        b = np.ones(d) * bounds_max

    if smooth_fun is None:
        smooth_fun = lambda p, lam: (math.pi / 2 - np.arctan(p - lam))
    #smooth_fun = lambda p, lam: np.exp(-10*(p - lam))

    # We do not need to store the cores, only the interfaces!
    Rx = [[]] * (d + 1)  # Python list for the interfaces
    Rx[0] = np.ones((1, 1))
    Rx[d] = np.ones((1, 1))
    Jy = [np.empty(0, dtype=np.int)] * (d + 1)
    ry = rmax * np.ones(d + 1, dtype=np.int)
    ry[0] = 1
    ry[d] = 1
    n = n0 * np.ones(d, dtype=np.int)
    fun_evals = 0

    grid = [np.reshape(np.linspace(a[i], b[i], n[i]), (n[i], 1))
            for i in xrange(d)]
    for i in xrange(d - 1):
        #cr1 = y[i]
        ry[i + 1] = min(ry[i + 1], n[i] * ry[i])
        cr1 = np.random.randn(ry[i], n[i], ry[i + 1])
        cr1 = reshape(cr1, (ry[i] * n[i], ry[i + 1]))
        q, r = np.linalg.qr(cr1)
        ind = maxvol(q)
        w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
        w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
        Jy[i + 1] = np.hstack((w1, w2))
        Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
        Jy[i + 1] = Jy[i + 1][ind, :]

        # Jy{i+1} = [kron(ones(n(i),1), Jy{i}), kron((1:n(i))', ones(ry(i),1))];
        # Jy{i+1} = Jy{i+1}(ind,:);
    swp = 0
    dirn = -1
    i = d - 1
    lm = float('Inf')
    while swp < nswp:
        # Right-to-left sweep
        # The idea: compute the current core; compute the function of it;
        # Shift locally or globally? Local shift would be the first try
        # Compute the current core

        if np.size(Jy[i]) == 0:
            w1 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w1 = mkron(np.ones((n[i] * ry[i + 1], 1), dtype=np.int), Jy[i])
        w2 = mkron(mkron(np.ones((ry[i + 1], 1), dtype=np.int),
                         grid[i]), np.ones((ry[i], 1), dtype=np.int))
        if np.size(Jy[i + 1]) == 0:
            w3 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w3 = mkron(Jy[i + 1], np.ones((ry[i] * n[i], 1), dtype=np.int))

        J = np.hstack((w1, w2, w3))
        # Just add some random indices to J, which is rnr x d, need to make rn (r + r0) x add,
        # i.e., just generate random r, random n and random multiindex

        cry = fun(J)
        fun_evals += cry.size
        cry = reshape(cry, (ry[i], n[i], ry[i + 1]))
        min_cur = np.min(cry.flatten("F"))
        ind_cur = np.argmin(cry.flatten("F"))
        if lm > min_cur:
            lm = min_cur
            x_full = J[ind_cur, :]
            val = fun(x_full)
            if verb:
                print('New record:', val, 'Point:', x_full, 'fevals:', fun_evals)
        cry = smooth_fun(cry, lm)
        if (dirn < 0 and i > 0):
            cry = reshape(cry, (ry[i], n[i] * ry[i + 1]))
            cry = cry.T
            #q, r = np.linalg.qr(cry)
            u, s, v = mysvd(cry, full_matrices=False)
            ry[i] = min(ry[i], rmax)
            q = u[:, :ry[i]]
            ind = rect_maxvol(q)[0]  # maxvol(q)
            ry[i] = ind.size
            w1 = mkron(np.ones((ry[i + 1], 1), dtype=np.int), grid[i])
            if np.size(Jy[i + 1]) == 0:
                w2 = np.zeros((n[i] * ry[i + 1], 0), dtype=np.int)
            else:
                w2 = mkron(Jy[i + 1], np.ones((n[i], 1), dtype=np.int))
            Jy[i] = np.hstack((w1, w2))
            Jy[i] = reshape(Jy[i], (n[i] * ry[i + 1], -1))
            Jy[i] = Jy[i][ind, :]

        if (dirn > 0 and i < d - 1):
            cry = reshape(cry, (ry[i] * n[i], ry[i + 1]))
            q, r = np.linalg.qr(cry)
            #ind = maxvol(q)
            ind = rect_maxvol(q)[0]
            ry[i + 1] = ind.size
            w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
            w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
            Jy[i + 1] = np.hstack((w1, w2))
            Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
            Jy[i + 1] = Jy[i + 1][ind, :]

        i += dirn
        if i == d or i == -1:
            dirn = -dirn
            i += dirn
            swp = swp + 1
    return val, x_full