Run normal FBA and flux minimization on model.
541788
def run_fba_minimized(self, reaction):
        """"""

        epsilon = self._args.epsilon
        solver = self._get_solver()

        p = fluxanalysis.FluxBalanceProblem(self._mm, solver)

        start_time = time.time()

        # Maximize reaction flux
        try:
            p.maximize(reaction)
        except fluxanalysis.FluxBalanceError as e:
            self.report_flux_balance_error(e)

        fluxes = {r: p.get_flux(r) for r in self._mm.reactions}

        # Run flux minimization
        flux_var = p.get_flux_var(reaction)
        p.prob.add_linear_constraints(flux_var == p.get_flux(reaction))
        p.minimize_l1()

        logger.info('Solving took {:.2f} seconds'.format(
            time.time() - start_time))

        count = 0
        for reaction_id in self._mm.reactions:
            flux = p.get_flux(reaction_id)
            if abs(flux - fluxes[reaction_id]) > epsilon:
                count += 1
            yield reaction_id, flux
        logger.info('Minimized reactions: {}'.format(count))
1259550
def optimize(self, x0, f=None, df=None, f_df=None):
        
        import scipy.optimize
        if f_df is None and df is not None: f_df = lambda x: float(f(x)), df(x)
        if f_df is not None:
            def _f_df(x):
                return f(x), f_df(x)[1][0]
        if f_df is None and df is None:
            res = scipy.optimize.fmin_l_bfgs_b(f, x0=x0, bounds=self.bounds,approx_grad=True, maxiter=self.maxiter)
        else:
            res = scipy.optimize.fmin_l_bfgs_b(_f_df, x0=x0, bounds=self.bounds, maxiter=self.maxiter)

        ### --- We check here if the the optimizer moved. If it didn't we report x0 and f(x0) as scipy can return NaNs
        if res[2]['task'] == b'ABNORMAL_TERMINATION_IN_LNSRCH':
            result_x  = np.atleast_2d(x0)
            result_fx =  np.atleast_2d(f(x0))
        else:
            result_x = np.atleast_2d(res[0])
            result_fx = np.atleast_2d(res[1])

        return result_x, result_fx
164633
def minimize():
    
    ########################### HYPER PARAMETERS ###############################

    # Stop when r2 is this much worse than full set of positions
    r2_stopping_percent = 0.96
    # for this many iterations
    stopping_iterations = 5

    # Limit SVM to a smaller number of positions to speed up code.
    max_positions_fit = 300
    # Filter any position that "contributes" less than this percent of max.
    filter_contribution_percent = 0.3
    # Never filter more than this many positions in one iterations
    filter_limit = 25

    ########################### HYPER PARAMETERS ###############################

    # Load positons
    model_paths = oneoff_utils.get_model_paths(FLAGS.model_dir)
    num_models = len(model_paths)
    assert num_models > 0, FLAGS.model_dir

    # Load model ratings
    # wget https://cloudygo.com/v12-19x19/json/ratings.json
    ratings = json.load(open(FLAGS.rating_json))
    raw_ratings = {int(r[0]): float(r[1]) for r in ratings}

    model_ratings = []
    for model in model_paths:
        model_idx = get_model_idx(model)
        if model_idx < FLAGS.min_idx:
            continue

        model_ratings.append(raw_ratings[model_idx])
    model_ratings = np.array(model_ratings)

    assert 0 < len(model_ratings) <= num_models, len(model_ratings)
    num_models = len(model_ratings)

    sgf_names, all_positions = get_final_positions()
    # Trim off common path prefix.
    common_path = os.path.commonpath(sgf_names)
    sgf_names = [name[len(common_path) + 1:] for name in sgf_names]

    print("Considering {} positions, {} models".format(
        len(all_positions), num_models))
    print()

    # Load model data
    top_n = FLAGS.top_n
    positions = defaultdict(list)
    with open(FLAGS.results) as results:
        headers = results.readline().strip()
        assert headers.count(",") + 1 == len(sgf_names)

        # Row is <model_name> + positions x [value, top_n x [move, move_policy]]
        for row in tqdm(results.readlines(), desc="result line"):
            data = row.split(",")
            model_idx = get_model_idx(data.pop(0))
            if model_idx < FLAGS.min_idx:
                continue

            data_per = 1 + top_n * 2
            assert len(data) % data_per == 0, len(data)

            for position, position_data in enumerate(grouper(data_per, data)):
                value = float(position_data.pop(0))
                moves = list(map(int, position_data[0::2]))
                move_policy = list(map(float, position_data[1::2]))

                positions[position].append([value, moves, move_policy])

    def one_hot(n, i):
        one_hot = [0] * n
        if 0 <= i < n:
            one_hot[i] += 1
        return one_hot

    # NOTE: top_n isn't the same semantic value here and can be increased.
    one_hot_moves = top_n
    num_features = 1 + 5 + (one_hot_moves + 1)

    # Features by position
    features = []
    pos_top_moves = []
    for position, data in tqdm(positions.items(), desc="featurize"):
        assert len(data) == num_models, len(data)

        top_moves = Counter([d[1][0] for d in data])
        top_n_moves = [m for m, c in top_moves.most_common(one_hot_moves)]
        if len(top_n_moves) < one_hot_moves:
            top_n_moves.extend([-1] * (one_hot_moves - len(top_n_moves)))
        assert len(top_n_moves) == one_hot_moves, "pad with dummy moves"
        pos_top_moves.append(top_n_moves)

        # Eventaully we want
        # [model 1 position 1 features, m1 p2 features, m1 p3 features, ... ]
        # [model 2 position 1 features, m2 p2 features, m2 p3 features, ... ]
        # [model 3 position 1 features, m3 p2 features, m3 p3 features, ... ]
        # ...
        # [model m position 1 features, mm p2 features, mm p3 features, ... ]

        # We'll do position selection by joining [model x position_feature]

        feature_columns = []
        for model, (v, m, mv) in enumerate(data):
            # Featurization (for each positions):
            #   * Value (-1 to 1), Bucketed value
            #   * Cluster all model by top_n moves (X,Y,Z or other)?
            #     * value of that move for model
            #   * policy value of top move
            model_features = []

            model_features.append(2 * v - 1)
            # NOTE(sethtroisi): Consider bucketize value by value percentiles.
            value_bucket = np.searchsorted((0.2, 0.4, 0.6, 0.8), v)
            model_features.extend(one_hot(5, value_bucket))

            # Policy weight for most common X moves (among all models).
            policy_weights = [0] * (one_hot_moves + 1)
            for move, policy_value in zip(m, mv):
                if move in top_n_moves:
                    policy_weights[top_n_moves.index(move)] = policy_value
                else:
                    policy_weights[-1] += policy_value
            model_features.extend(policy_weights)

            assert len(model_features) == num_features

            feature_columns.append(model_features)
        features.append(feature_columns)

    features = np.array(features)
    print("Feature shape", features.shape)
    print()

    # Split the models to test / train
    train_size = int(num_models * 0.9)
    train_models = sorted(np.random.permutation(num_models)[:train_size])
    test_models = sorted(set(range(num_models)) - set(train_models))
    assert set(train_models + test_models) == set(range(num_models))
    features_train = features[:, train_models, :]
    features_test  = features[:, test_models, :]

    labels_train = model_ratings[train_models]
    labels_test = model_ratings[test_models]

    # Choose some set of positions and see how well they explain ratings
    positions_to_use = set(positions.keys())
    linearSVM = svm.LinearSVR()
    best_test_r2 = 0
    below_threshold = 0

    for iteration in itertools.count(1):
        iter_positions = np.random.permutation(list(positions_to_use))
        iter_positions = sorted(iter_positions[:max_positions_fit])

        # Take this set of positions and build X
        X = np.concatenate(features_train[iter_positions], axis=1)
        Xtest = np.concatenate(features_test[iter_positions], axis=1)
        assert X.shape == (train_size, num_features * len(iter_positions))

        linearSVM.fit(X, labels_train)

        score_train = linearSVM.score(X, labels_train)
        score_test = linearSVM.score(Xtest, labels_test)
        print("iter {}, {}/{} included, R^2: {:.4f} train, {:.3f} test".format(
            iteration, len(iter_positions), len(positions_to_use),
            score_train, score_test))

        # Determine the most and least useful position:
        # TODO(amj,brilee): Validate this math.
        assert len(linearSVM.coef_) == num_features * len(iter_positions)

        # The intercepts tell us how much this contributes to overall rating
        # but coef tell us how much different answers differentiate rating.
        coef_groups = list(grouper(num_features, linearSVM.coef_))
        position_coefs = [abs(sum(c)) for c in coef_groups]

        pos_value_idx = np.argsort(position_coefs)
        max_pos = pos_value_idx[-1]
        most_value = position_coefs[max_pos]

        print("\tMost value {} => {:.1f} {}".format(
            max_pos, most_value, sgf_names[iter_positions[max_pos]]))

        # Drop any positions that aren't very useful
        for dropped, pos_idx in enumerate(pos_value_idx[:filter_limit], 1):
            contribution = position_coefs[pos_idx]
            positions_to_use.remove(iter_positions[pos_idx])
            print("\t\tdropping({}): {:.1f} {}".format(
                dropped, contribution, sgf_names[iter_positions[pos_idx]]))

            if contribution > filter_contribution_percent * most_value:
                break
        print()

        best_test_r2 = max(best_test_r2, score_test)
        if score_test > r2_stopping_percent * best_test_r2:
            below_threshold = 0
        else:
            below_threshold += 1
            if below_threshold == stopping_iterations:
                print("{}% decrease in R^2, stopping".format(
                    100 - int(100 * r2_stopping_percent)))
                break

    # Write down the differentiating positions and their answers.
    svm_data = []
    for position_idx in list(reversed(pos_value_idx)):
        coefs = coef_groups[position_idx]

        # Global position index.
        position = iter_positions[position_idx]
        sgf_name = sgf_names[position]
        top_moves = pos_top_moves[position]

        svm_data.append([sgf_name, [top_moves, coefs.tolist()]])

    with open(FLAGS.SVM_json, "w") as svm_json:
        json.dump(svm_data, svm_json)
    print("Dumped data about {} positions to {}".format(
        len(svm_data), FLAGS.SVM_json))
1065015
def run_tfba(self, reaction):
        """"""

        solver = self._get_solver(integer=True)
        p = fluxanalysis.FluxBalanceProblem(self._mm, solver)

        start_time = time.time()

        p.add_thermodynamic()

        try:
            p.maximize(reaction)
        except fluxanalysis.FluxBalanceError as e:
            self.report_flux_balance_error(e)

        logger.info('Solving took {:.2f} seconds'.format(
            time.time() - start_time))

        for reaction_id in self._mm.reactions:
            yield reaction_id, p.get_flux(reaction_id)
1210858
def model_exchange(model):
    """"""
    # Determine the default flux limits. If the value is already at the
    # default it does not need to be included in the output.
    lower_default, upper_default = None, None
    if model.default_flux_limit is not None:
        lower_default = -model.default_flux_limit
        upper_default = model.default_flux_limit

    compounds = []
    for compound, reaction_id, lower, upper in sorted(
            itervalues(model.exchange)):
        d = OrderedDict([('id', compound.name)])
        if reaction_id is not None:
            d['reaction'] = reaction_id

        lower = _get_output_limit(lower, lower_default)
        upper = _get_output_limit(upper, upper_default)
        d.update(_generate_limit_items(lower, upper))

        compounds.append(d)

    return OrderedDict([('compounds', compounds)])
706755
def admm_linearized_simple(x, f, g, L, tau, sigma, niter, **kwargs):
    
    callback = kwargs.pop('callback', None)
    z = L.range.zero()
    u = L.range.zero()
    for _ in range(niter):
        x[:] = f.proximal(tau)(x - tau / sigma * L.adjoint(L(x) + u - z))
        z = g.proximal(sigma)(L(x) + u)
        u = L(x) + u - z
        if callback is not None:
            callback(x)
553179
def flux_v2(v_vars: List[fl.Var], i: int):
    """"""
    # Index with the base of (v_x, v_y, v_z) for body i
    k = 3*i
    # The speed squared of body i
    return fl.square(v_vars[k+0]) + fl.square(v_vars[k+1]) + fl.square(v_vars[k+2])
1321509
def fit(self, inputs, outputs, epochs=1):
        
        self.model.fit(inputs, outputs, epochs=epochs)
1287251
def _run_opf(self, P, q, AA, ll, uu, xmin, xmax, x0, opt):
        
        AAcvx = tocvx(AA)
        nx = x0.shape[0] # number of variables
        # add var limits to linear constraints
        eyex = spmatrix(1.0, range(nx), range(nx))
        A = sparse([eyex, AAcvx])
        l = matrix([-xmin, ll])
        u = matrix([ xmax, uu])

        Ae, be, Ai, bi = split_linear_constraints(A, l, u)


        if len(P.V) > 0:
            cvx_sol = solvers.qp(P, q, Ai, bi, Ae, be, self.solver, {"x": x0})
        else:
            cvx_sol = solvers.lp(q, Ai, bi, Ae, be, self.solver, {"x": x0})

        return cvx_sol
360705
def _minimize_constraints_fun_summation(x):
    '''
    
    '''
    summation = sum([x[i] for i in CONSTRAINT_PARAMS_IDX])
    return CONSTRAINT_UPPERBOUND >= summation >= CONSTRAINT_LOWERBOUND
333530
def _adapt_params(self, stepsize, stepsize_bar, h_bar, mu, index_i, alpha, n_alpha=1):
        
        gamma = 0.05  # free parameter that controls the amount of shrinkage towards mu
        t0 = 10.0  # free parameter that stabilizes the initial iterations of the algorithm
        kappa = 0.75
        # See equation (6) section 3.2.1 for details

        estimate = 1.0 / (index_i + t0)
        h_bar = (1 - estimate) * h_bar + estimate * (self.delta - alpha / n_alpha)

        stepsize = np.exp(mu - sqrt(index_i) / gamma * h_bar)
        i_kappa = index_i ** (-kappa)
        stepsize_bar = np.exp(i_kappa * np.log(stepsize) + (1 - i_kappa) * np.log(stepsize_bar))

        return stepsize, stepsize_bar, h_bar
1072266
def get_fine_tune_model(symbol, arg_params, num_classes, layer_name, dtype='float32'):
    
    all_layers = symbol.get_internals()
    net = all_layers[layer_name+'_output']
    net = mx.symbol.FullyConnected(data=net, num_hidden=num_classes, name='fc')
    if dtype == 'float16':
        net = mx.sym.Cast(data=net, dtype=np.float32)
    net = mx.symbol.SoftmaxOutput(data=net, name='softmax')
    new_args = dict({k:arg_params[k] for k in arg_params if 'fc' not in k})
    return (net, new_args)
1136080
def CDPP(flux, mask=[], cadence='lc'):
    '''
    

    '''

    # 13 cadences is 6.5 hours
    rmswin = 13
    # Smooth the data on a 2 day timescale
    svgwin = 49

    # If short cadence, need to downbin
    if cadence == 'sc':
        newsize = len(flux) // 30
        flux = Downbin(flux, newsize, operation='mean')

    flux_savgol = SavGol(np.delete(flux, mask), win=svgwin)
    if len(flux_savgol):
        return Scatter(flux_savgol / np.nanmedian(flux_savgol),
                       remove_outliers=True, win=rmswin)
    else:
        return np.nan
774009
def fc(self, output_nodes, keep_prob=1, activation_fn=tf.nn.relu, b_value=0.0, s_value=1.0, bn=True,
           trainable=True):
        
        self.count['fc'] += 1
        scope = 'fc_' + str(self.count['fc'])
        with tf.variable_scope(scope):

            # Flatten if necessary
            if len(self.input.get_shape()) == 4:
                input_nodes = tf.Dimension(
                    self.input.get_shape()[1] * self.input.get_shape()[2] * self.input.get_shape()[3])
                output_shape = tf.stack([-1, input_nodes])
                self.input = tf.reshape(self.input, output_shape)

            # Matrix Multiplication Function
            input_nodes = self.input.get_shape()[1]
            output_shape = [input_nodes, output_nodes]
            w = self.weight_variable(name='weights', shape=output_shape, trainable=trainable)
            self.input = tf.matmul(self.input, w)

            if bn is True:  # batch normalization
                self.input = self.batch_norm(self.input, 'fc')
            if b_value is not None:  # bias value
                b = self.const_variable(name='bias', shape=[output_nodes], value=b_value, trainable=trainable)
                self.input = tf.add(self.input, b)
            if s_value is not None:  # scale value
                s = self.const_variable(name='scale', shape=[output_nodes], value=s_value, trainable=trainable)
                self.input = tf.multiply(self.input, s)
            if activation_fn is not None:  # activation function
                self.input = activation_fn(self.input)
            if keep_prob != 1:  # dropout function
                self.input = tf.nn.dropout(self.input, keep_prob=keep_prob)
        print(scope + ' output: ' + str(self.input.get_shape()))
327560
private void minimizeWithGeometricConstraints() {
		extractEpipoles.setTensor(solutionN);
		extractEpipoles.extractEpipoles(e2,e3);

		// encode the parameters being optimized
		param[0] = e2.x; param[1] = e2.y; param[2] = e2.z;
		param[3] = e3.x; param[4] = e3.y; param[5] = e3.z;

		// adjust the error function for the current inputs
		errorFunction.init();

		// set up the optimization algorithm
		optimizer.setFunction(errorFunction,null);
		optimizer.initialize(param,gtol,ftol);

		// optimize until convergence or the maximum number of iterations
		UtilOptimize.process(optimizer, maxIterations);

		// get the results and compute the trifocal tensor
		double found[] = optimizer.getParameters();
		paramToEpipoles(found,e2,e3);

		enforce.process(e2,e3,A);
		enforce.extractSolution(solutionN);
	}
1354590
def _get_flux_bounds(self, r_id, model, flux_limits, equation):
        
        if r_id not in flux_limits or flux_limits[r_id][0] is None:
            if equation.direction == Direction.Forward:
                lower = 0
            else:
                lower = -model.default_flux_limit
        else:
            lower = flux_limits[r_id][0]

        if r_id not in flux_limits or flux_limits[r_id][1] is None:
            if equation.direction == Direction.Reverse:
                upper = 0
            else:
                upper = model.default_flux_limit
        else:
            upper = flux_limits[r_id][1]

        if lower % 1 == 0:
            lower = int(lower)
        if upper % 1 == 0:
            upper = int(upper)
        return text_type(lower), text_type(upper)
86235
def minimise(objective_function, p0, **kwargs):
    

    p0 = np.array(p0)

    op_kwargs = kwargs.copy()

    # Which optimisation algorithm?
    available_methods = ("BFGS", "Nelder-Mead", "Powell", "CG", "TNC")
    method = op_kwargs.pop("method", available_methods[0])
    logger.debug("Optimising with {0} algorithm".format(method))

    op_kwargs.update({ "callback": None, "full_output": True, "retall": False })

    def _drop_forbidden_keywords(op_kwds, allowed):
        for keyword in set(op_kwds).difference(allowed):
            logger.debug("Dropping optimisation keyword {0}: {1}".format(
                keyword, op_kwds.pop(keyword)))

    t_init = time()
    if method == "Nelder-Mead":
        # Need to force some arguments.
        # Optional arguments:
        # xtol, ftol, maxiter, maxfun
        _drop_forbidden_keywords(op_kwargs, ("xtol", "ftol", "maxiter", 
            "maxfun", "full_output", "retall", "callback", "disp"))

        # Set defaults.
        op_kwargs.setdefault("disp", False)

        x_opt, f_opt, num_iter, num_funcalls, warnflag \
            = op.fmin(objective_function, p0, **op_kwargs)

        logger.debug("Number of iterations: {0}, number of function calls: {1}"\
            .format(num_iter, num_funcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Maximum number of function evaluations.",
                "Maximum number of iterations."][warnflag - 1]))

    elif method == "Powell":
        # Optional arguments: xtol, ftol, maxiter, maxfun
        _drop_forbidden_keywords(op_kwargs, ("xtol", "ftol", "maxiter", 
            "maxfun", "retall", "callback", "full_output", "disp"))
        
        # Set defaults.
        op_kwargs.setdefault("disp", False)

        x_opt, f_opt, direc, num_iter, num_funcalls, warnflag \
            = op.fmin_powell(objective_function, p0, **op_kwargs)

        logger.debug("Number of iterations: {0}, number of function calls: {1}"\
            .format(num_iter, num_funcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Maximum number of function evaluations.",
                "Maximum number of iterations."][warnflag - 1]))

    elif method == "CG":
        # Optional arguments: gtol, norm, epsilon, maxiter.
        _drop_forbidden_keywords(op_kwargs, ("gtol", "norm", "epsilon", 
            "maxiter", "retall", "callback", "full_output"))

        x_opt, f_opt, num_funcalls, num_gradcalls, warnflag \
            = op.fmin_cg(objective_function, p0, **op_kwargs)

        logger.debug("Number of function calls: {0}, gradient calls: {1}"\
            .format(num_funcalls, num_gradcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.debug("Optimisation failed: {0}".format([
                "Maximum number of iterations exceeded",
                "Gradient and/or function calls were not changing."
                ][warnflag - 1]))

    elif method == "BFGS":
        # Since we have at least some boundaries, this will actually call
        # the L-BFGS-B algorithm.
        # Optional arguments: m, factr, pgtol, epsilon, maxfun, maxiter

        _drop_forbidden_keywords(op_kwargs, ("m", "factr", "pgtol", 
            "approx_grad", "epsilon", "disp", "iprint", "maxfun", "maxiter"))

        # Default/required:
        op_kwargs.setdefault("factr", 10.0)
        op_kwargs["approx_grad"] = True
        
        # Because the parameters will vary by orders of magnitude, here we
        # scale everything to the initial value so that the epsilon keyword
        # makes some sense.

        scale = p0.copy()
        
        def scaled_objective_function(theta):
            return objective_function(theta.copy() * scale)

        x_opt, f_opt, info_dict = op.fmin_l_bfgs_b(scaled_objective_function,
            np.ones(p0.size, dtype=float), **op_kwargs)

        # Rescale.
        x_opt *= scale

        g_opt = info_dict["grad"]
        num_iter = info_dict["nit"]
        num_funcalls = info_dict["funcalls"]
        warnflag = info_dict["warnflag"]

        logger.debug("Number of iterations: {0}, function calls: {1}".format(
            num_iter, num_funcalls))
        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Too many function evaluations or iterations!",
                "{0}".format(info_dict.get("task", None))][warnflag - 1]))
    
    elif method == "TNC":
        # Optional arguments:   bounds, epsilon, scale, offset, messages, 
        #                       maxCGit, maxfun, eta, stepmx, accuracy, 
        #                       fmin, ftol, xtol, pgtol, rescale

        _drop_forbidden_keywords(op_kwargs, ("approx_grad", "bounds", "epsilon",
            "scale", "offset", "messages", "maxCGit", "maxfun", "eta", "stepmx",
            "accuracy", "fmin", "ftol", "xtol", "pgtol", "rescale", "disp"))

        # Required:
        op_kwargs["approx_grad"] = True

        x_opt, num_funcalls, rc = op.fmin_tnc(objective_function, p0, 
            **op_kwargs)
        
        rcstring = {
            -1: "Infeasible (lower bound > upper bound)",
            0: "Local minimum reached (|pg| ~= 0)",
            1: "Converged (|f_n-f_(n-1)| ~= 0)",
            2: "Converged (|x_n-x_(n-1)| ~= 0)",
            3: "Max. number of function evaluations reached",
            4: "Linear search failed",
            5: "All lower bounds are equal to the upper bounds",
            6: "Unable to progress",
            7: "User requested end of minimization"
        }[rc]

        logger.debug("Number of function calls: {0}, result: {1}".format(
            num_funcalls, rcstring))
        if rc in (1, 2):
            logger.debug("Optimisation converged after {0:.1f} seconds: {1}"\
                .format(time() - t_init, rcstring))
        else:
            logger.warn("Optimisation failed: {0}".format(rcstring))

    else:
        raise ValueError("optimisation algorithm {0} is not available "\
            "(available methods are {1})".format(
                method, ", ".join(available_methods)))


    return x_opt