@@ -18,7 +18,7 @@ flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel} =

...

@@ -18,7 +18,7 @@ flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel} =

Run flux balance analysis (FBA) on the `model` optionally specifying `objective_rxn(s)` and their `weights` (empty `weights` mean equal weighting per reaction).

Run flux balance analysis (FBA) on the `model` optionally specifying `objective_rxn(s)` and their `weights` (empty `weights` mean equal weighting per reaction).

Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.

Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.

Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.

Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.

The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.

The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.

This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.

This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.

Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.

Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.

@@ -116,7 +116,7 @@ Run flux variability analysis (FVA) on the `model` with `objective_rxn(s)` and o

...

@@ -116,7 +116,7 @@ Run flux variability analysis (FVA) on the `model` with `objective_rxn(s)` and o

It runs fba on the model once to determine the optimum of the objective.

It runs fba on the model once to determine the optimum of the objective.

Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.

Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (ub, lb) flux constraints.

The model is then constrained to produce objective flux bounded by `optimum_bound` from below (set to slightly less than 1.0 for stability) and each flux in the model sequentially minimized and maximized.

The model is then constrained to produce objective flux bounded by `optimum_bound` from below (set to slightly less than 1.0 for stability) and each flux in the model sequentially minimized and maximized.

Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.

Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.

The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.

The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.

This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.

This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.

Returns two dictionaries (`fva_max` and `fva_min`) that each reaction `id`s to dictionaries of the resultant flux distributions (if solved successfully) when that `id` is optimized.

Returns two dictionaries (`fva_max` and `fva_min`) that each reaction `id`s to dictionaries of the resultant flux distributions (if solved successfully) when that `id` is optimized.

Perform basic hit and run sampling for `N` iterations using `model` with `optimizer` from `JuMP`.

Perform basic hit and run sampling for `N` iterations using `model` with `optimizer` from `JuMP`.

Additional constraints supplied by `constraints` as a dictionary of reaction `id`s mapped to a tuple of `(lb, ub)` of fluxes.

Additional constraints supplied by `constraints` as a dictionary of reaction `id`s mapped to a tuple of `(lb, ub)` of fluxes.

Every `keepevery` iteration is logged as a sample, where the sample size matrix has `samplesize` columns.

Every `keepevery` iteration is logged as a sample, where the sample size matrix has `samplesize` columns.

Solver specific settings can be set using `solver_attributes`.

Solver specific settings can be set using `solver_attributes`.

Warm up points are generated in a flux variability sense, unless `random_objective` is true, in which case a randomly weighted objective is used 2*number of reactions to define the warmup points.

Warm up points are generated in a flux variability sense, unless `random_objective` is true, in which case a randomly weighted objective is used 2*number of reactions to define the warmup points.

Note that N needs to be >> samplesize.

Note that N needs to be >> samplesize.

Sample size is the size of the samples kept in memory.

Sample size is the size of the samples kept in memory.

The larger samplesize is the better the approximation becomes, but the more memory the sampler requires.

The larger samplesize is the better the approximation becomes, but the more memory the sampler requires.

See also: [`achr`](@ref)

See also: [`achr`](@ref)

...

@@ -126,13 +126,13 @@ function hit_and_run(

...

@@ -126,13 +126,13 @@ function hit_and_run(

current_point=zeros(size(wpoints,1))

current_point=zeros(size(wpoints,1))

current_point.=wpoints[:,rand(1:nwpts)]# pick random initial point

current_point.=wpoints[:,rand(1:nwpts)]# pick random initial point

δdirtol=1e-6# too small directions get ignored ≈ 0 (solver precision issue)

δdirtol=1e-6# too small directions get ignored ≈ 0 (solver precision issue)

sample_num=0

sample_num=0

samplelength=0

samplelength=0

updatesamplesizelength=true

updatesamplesizelength=true

forn=1:N

forn=1:N

# direction = random point - current point

# direction = random point - current point

ifupdatesamplesizelength

ifupdatesamplesizelength

direction_point=(@viewwpoints[:,rand(1:nwpts)])-(@viewcurrent_point[:])# use warmup points to find direction in warmup phase

direction_point=(@viewwpoints[:,rand(1:nwpts)])-(@viewcurrent_point[:])# use warmup points to find direction in warmup phase

else

else

...

@@ -157,7 +157,7 @@ function hit_and_run(

...

@@ -157,7 +157,7 @@ function hit_and_run(

upper=1e10

upper=1e10

end

end

lower>λmin&&(λmin=lower)# max min step size that satisfies all bounds

lower>λmin&&(λmin=lower)# max min step size that satisfies all bounds

upper<λmax&&(λmax=upper)# min max step size that satisfies all bounds

upper<λmax&&(λmax=upper)# min max step size that satisfies all bounds

end

end

ifλmax<=λmin||λmin==-1e10||λmax==1e10# this sometimes can happen

ifλmax<=λmin||λmin==-1e10||λmax==1e10# this sometimes can happen

...

@@ -260,12 +260,12 @@ function achr(

...

@@ -260,12 +260,12 @@ function achr(

shat=mean(wpoints,dims=2)[:]# mean point

shat=mean(wpoints,dims=2)[:]# mean point

δdirtol=1e-6# too small directions get ignored ≈ 0 (solver precision issue)

δdirtol=1e-6# too small directions get ignored ≈ 0 (solver precision issue)

sample_num=0

sample_num=0

samplelength=0

samplelength=0

updatesamplesizelength=true

updatesamplesizelength=true

forn=1:N

forn=1:N

ifupdatesamplesizelength# switch to samples

ifupdatesamplesizelength# switch to samples

direction_point=(@viewwpoints[:,rand(1:nwpts)])-(@viewcurrent_point[:])# use warmup points to find direction in warmup phase

direction_point=(@viewwpoints[:,rand(1:nwpts)])-(@viewcurrent_point[:])# use warmup points to find direction in warmup phase

else

else

direction_point=(@viewsamples[:,rand(1:(samplelength))])-(@viewshat[:])# after warmup phase, only find directions in sampled space

direction_point=(@viewsamples[:,rand(1:(samplelength))])-(@viewshat[:])# after warmup phase, only find directions in sampled space

...

@@ -287,7 +287,7 @@ function achr(

...

@@ -287,7 +287,7 @@ function achr(

upper=1e10

upper=1e10

end

end

lower>λmin&&(λmin=lower)# max min step size that satisfies all bounds

lower>λmin&&(λmin=lower)# max min step size that satisfies all bounds

upper<λmax&&(λmax=upper)# min max step size that satisfies all bounds

upper<λmax&&(λmax=upper)# min max step size that satisfies all bounds

end

end

ifλmax<=λmin||λmin==-1e10||λmax==1e10# this sometimes can happen

ifλmax<=λmin||λmin==-1e10||λmax==1e10# this sometimes can happen