-
Breaking:
define(),insert(),repar()andbranch(apply = )write a substitution as"lhs <- rhs", e.g.insert(trafo, "x <- exp(x)", x = pars). A"lhs ~ rhs"string or a formula is an error that names the new form;~is left to model formulas. -
evalConditionResidual()withdfgives the Student-t likelihood ofnormT(). -
Random effects as formulas:
insert()on a branched transformation takes effect formulas such ascl ~ wt + (1 + dose | p | ID). The symbol is replaced by its linear predictor on the scale the transformation gives it, so afterinsert(trafo, "x <- exp(x)", x = "cl")the random effects are normal on the log scale. Fixed effects are coded as bymodel.matrix()over the columns of the branch table, random terms follow lme4 with the IDs of brms for correlation across formulas,gr(ID, dist = )sets the distribution.x ~ 1 + (1 | ID), x = parsandka + v ~ ...apply one right-hand side to many symbols.randomEffects()returns the structure from the transformation or from theP()built from it. -
New
normal(),laplace(lambda),spikeSlab(prob, tau)anddirichletProcess(concentration, tau), the distributions of a random term. A spike-and-slab deviation is exactly 0 or, with probabilityprob, drawn from N(0, tau^2). UnderdirichletProcess()the groups share the population value or a common deviation in clusters, with the Ewens prior over the partitions and N(0, tau^2) as base measure. Every hyperparameter leftNULLis estimated by marginal maximum likelihood with the other parameters. -
New
em()andmsem(), mixed-effects fits astrust()andmstrust():em(obj, random = randomEffects(p), init, fixed, method)takes the objective of the data and builds the prior itself. A group may own several conditions, conditions without a group enter the fit, random terms over several groupings, nested or crossed, are integrated per connected component, and normal and sparse random effects combine in one fit. Methods:"focei", the Laplace approximation of the marginal with exact gradient, censored data by M3 and random effects in the error model; Laplace, spike-and-slab and Dirichlet-process random effects are integrated on the Schur complement of the normal ones, the partitions of a Dirichlet process summed exactly (at most 14 groups), and sparse random effects of one group that the data couple by expectation consistent (control$laplace$coupling, or mean field as a bound). Their hyperparameters move with the other parameters in one trust region: the Schur complement of the fixed parameters takes the sparse random effects at the curvature of their posterior, the curvature of a hyperparameter puts its Newton step at its M-step, and a BFGS correction learns the missing information. Such fits have no Hessian of the marginal: with Laplace and spike-and-slab random effectsvcov()inverts the empirical information of the exact scores of the groups,fit$information, which needs more groups than parameters; with Dirichlet-process ones, whose groups are not independent,vcov()is an error andconfint(method = "profile")gives intervals;"quadrature", from the FOCEI fit, on nested sparse grids whose level per group followscontrol$quadrature$tol, with an importance-sampling certificate per group;"laplaceEM", the same nodes in an ECM;"saem", stochastic approximation with the mu-referenced fixed effects and Omega in closed form; for Dirichlet-process random effects of any number of groups by Gibbs sweeps over the clusters, -2 log L by sequential importance sampling withvalueSE.
msem()perturbs the starts on the scale on which each parameter is unbounded and returns the best fit with the table of its starts, whichplotValues()draws as a waterfall;coressplits as inmstrust(). The fit frominitruns first, and every other start ends afterbudget(default 5) times its outer iterations, withconverged = FALSE. -
em()fits reportomega(standard deviations and correlations with standard errors),random(the modes per group with conditional standard errors, for the sparse distributions the cluster and the probability of 0),strengths,slab,logLik()(soAIC()andBIC()),residuals()of type IRES, IWRES, CWRES, PRES or PWRES,predict(type = ),confint(method = "profile"),anova(),simulate(),summary(),coef(),vcov(),plot(),plotResiduals(),plotTrace(),plotCombined()(data, IPRED and PRED per condition) and the new genericplotRandomEffects(), which compares draws from the conditional distributions of the random effects with their estimated normal, so the shrinkage of the modes does not distort it (type = "modes"for the modes). -
New
normT(): the Student-t likelihood with fixed degrees of freedom, with exact derivatives and censored data by M3, throughtrust()and every method ofem(). -
em()andmsem()takesparse = "auto"(default),TRUEorFALSE, astrust()does. -
em()by FOCEI over normal random effects runs its outer trust region on the sparse Hessian wheresparsechooses it: the joint Gauss-Newton Hessian minus the Schur complement of every group plus the volume term, each local to the parameters of the group, for one grouping and for several.summary()and the standard errors of Omega solve the sparse Hessian for the columns they need. -
The normal prior of
em()returns a sparse Hessian onhessian = "sparse". -
FOCEI, with normal and with sparse random effects, gets the second derivatives of its volume term by forward over reverse (
control$derivMode = "forward-reverse", the default, for models built withderivMode = c("forward", "reverse", "forward-reverse")): one pass forward along the random effects and a few backward sweeps give the exact gradient, the mode shift and the derivatives of the curvature of every group, at a cost that does not grow with the square of the parameters. The outer trust region takes the Gauss-Newton Schur complement with the volume curvature of the variance parameters from these derivatives; the reported Hessian is exact forward over forward where the model hasderiv2, else the exact Schur complement of the joint Hessian with that volume curvature, exact for models linear in the random effects, andsummary()says so."forward-forward"keeps the second-order sensitivities; a model with only one of the two takes it. In both modes, dense and sparse, the outer model of FOCEI gets an SR1 correction from the changes of the full outer gradient between the last three accepted points: it adds the curvature Gauss-Newton leaves out, negative curvature included, so fewer fits stop at saddles of the Gauss-Newton model. -
FOCEI in both derivative modes: the inner trust region of the random effects steps on the Gauss-Newton Hessian with an SR1 correction from its accepted steps where that stays positive definite; the outer gradient has the first-order correction for the eta gradient the inner solve leaves at the modes; a rejected step shrinks the radius to a quarter of its own length where that is shorter; the fit ends on the evaluation of its last accepted point. Forward over reverse answers all second-order seeds of a step in one backward sweep.
-
em()with Laplace, spike-and-slab or Dirichlet-process random effects takes the terms ofobjbesides the data, such as priors and constraints; a term on a hyperparameter or a random effect is an error. -
Sparse random effects by forward over reverse: the marginal alone needs no second-order sensitivities, so models built without
deriv2also return to expectation consistent after a fall back to mean field. -
em()says in a message when its models give second derivatives forward over forward only and namesderivMode = c("forward", "reverse", "forward-reverse");msem()says it once. -
Removed from the devel-EM layer:
EM()andmsEM()(nowem()andmsem()),omega(),penaltyL1(),emInit(),updateOmegaChol(),makeSubjectNodes(),emObjfn(),sparsify(),scanL1(),gateL1(),selectEB()and the vignette "Structure selection with an estimated penalty". -
New
alignmentModel()andalign(): alignment of replicate readings with free expressions for the prediction and its standard deviation. Every symbol is a linear predictor given by a formula inmodel.matrix()coding with random terms, a column with one level per value, numeric or not, andI(dose)a covariate.trafo = c("sj <- 10^sj")sets the scale of a symbol,error = NULLan error model without symbols,prioraconstraintL2()on any parameter, the spread of a random term included; censored readings by M3, a Student-t option. As in blotIt, the experiments connected by shared values form blocks, and per block the mean of the predictors of the first biological symbol is held at 0. Absolute readings, marked in a logical columnabsolute, such as concentrations by targeted mass spectrometry, are predicted by the expressionabsolutein the biological symbols and data columns (the first biological symbol by default, units by the expression), depend on no scaling parameter and may beNAin the columns of the scaling formulas; their block is not normalized but comes out in their units. Every experiment within its connected component of the groups of the random terms is a condition of its own, so the Hessian is sparse over them; designs and coefficients are sparse matrices and the standard errors of the values come from a sparse solve.align()starts the coefficients of the first scaling symbol where its predictor is, per experiment, the median of the predictors at which the prediction meets the relative readings, runs the fits of a multistart in parallel oncores, one core each, and stops at tolerances ofsqrt(.Machine$double.eps), inner fits included. Methodsprint(),as.datalist(),plotValues()(the starts by rank),plotResiduals(),anova()andsimulate(), plotsplotAligned()andplotScaled(). The demo "alignment" aligns qPCR experiments with undetected wells, and Western blots with a random background per gel together with targeted mass spectrometry, which puts the values in nM. -
alignmentModel()takes random terms of several groupings, nested or crossed, such as a lane effect per gel and lane shared by the antibodies of a gel beside a background per gel and antibody. They are normal and integrated per connected component; a condition is an experiment within its component, and the normalization blocks follow the fixed parts of the scaling formulas.fit$omegaandfit$randomhave a columngrouping, andplotRandomEffects()draws a panel per grouping and random effect. -
The objective of an alignment answers all its conditions by one prediction and one error-model function: the compiled model is evaluated once at the readings of all conditions, and the derivatives of every condition in its parameters and the contractions of the backward sweep run in C++, in parallel over the conditions. The error model of an alignment is taken from the prediction, which evaluates it alongside. Its backward sweeps answer several second-order seeds at once, so FOCEI fits it by forward over reverse.
-
An alignment holds the evaluation that the backward sweeps of prediction and error model share only within one sweep, so a fitted alignment keeps none of it; an error model that depends on no symbol is not evaluated backwards.
-
The reverse mode: the Hessian of
normL2()withsweep = "reverse"has the curvature of an estimated error model, an error model that is a sum of per-condition ones goes backwards, and so does an observation function without parameters of its own at second order.Y(),P()andXs()answer several second-order seeds in one backward sweep. -
A prediction or observation function whose kernel has a
conditionargument gets the condition of every call, its batch entries getconditions; an observation function answers the backward sweep through its batch entry where it has one.normL2()withtimesas a list evaluates a function of several conditions at the times of each, and its backward sweep sums the cotangents of the conditions in one pass. An error model whose output the prediction holds is taken from there. -
A cppDE model keeps the batch handles of its last eight request sets.
-
A sum of objectives holds its terms, not the frame it is built in.
-
distributedComputing(build = "compile")builds a model of eight or more C++ sources that lie outside the working directory, such as ingetOption("dMod2.outdir").
-
evalConditionResidual()has the parameters that enter only the error model in its gradient and Hessian. -
An objective called with
hessian = "sparse"returns its Gauss-Newton Hessian as a sparseMatrixassembled from the blocks of the conditions;trust(sparse = TRUE)factorises it by sparse Cholesky andvcov(which = )solves it for the columns named.datapointL2()returns a sparse Hessian there too. -
trust()andprofile()takesparse = "auto"(default),TRUEorFALSE."auto"takes the sparse Gauss-Newton Hessian for a minimisation without phases, bounds or exact Hessian when there are at least 500 free parameters and the Hessian at the start fills at most 5% of its matrix (options(dMod2.sparseAuto = )), and starts from that evaluation.trust(sparse = TRUE)on an objective that returns a dense Hessian is an error that names the terms of its sum that do.mstrust()passessparsetotrust(). -
profile()on a sparse Hessian decides by the pivots of its LU whether the path equation is singular, by the rule ofvcov(): a pivot at most 1e-14 times the largest. A singular step then takes the pseudo-inverse and reoptimises as on a dense Hessian, up to 5000 parameters; above, the error names the non-identifiable parameters. The pseudo-inverse step keeps all its components.algoControl$W = "identity"moves the profiled parameter alone and reoptimises the others at every step, without a matrix in the path equation and without a warning on singularity. -
Requires cppDE 0.13.4: a second-order solve whose parameters have no second derivatives crashed R with an older cppDE.
-
distributedComputing()keeps the manifest of a job also in its folder on every remote node.recover = TRUEwithongiven takes it from there, so a job submitted from another machine is recovered, and a job of the same name that still runs on the nodes is an error also when it was submitted from another machine. -
distributedComputing()runs on Windows, with the local commands in the sh of Rtools.detectFreeCores()measures the processor load of the local Windows machine, without a warning. -
compile()compiles in parallel on Windows too, with or without a combinedoutput, andP()andY()pass theircoresto it there.
-
distributedComputing()reports the start of a job as a message, on a SLURM cluster with the job ID and the array tasks. A build writes its output tobuild.login the job folder, and a failed build is an error that shows its end instead of starting the tasks. -
plotProfile()leaves more space between its panels, and a narrow panel leaves out the x labels that would overlap.
-
Breaking:
distributedComputing()takes where its tasks run ason:localNode(),sshNodes()orslurmCluster(), each with thecoresof a task. A node is one worker, and a SLURM job hasarrayworkers, its array tasks; the tasks are dealt to the workers in turn.slurmCluster()takeswalltime,memPerCore, the job script andpartition,account,maxRunningand further options ofsbatch.build = c("auto", "compile", "link", "copy")replacescompileandlink, and"auto"compiles on remote nodes. The argumentsmachine,slurm,jobscript,memPerCore,walltime,compileandlinkare gone, andslurmTemplate()names its first argumentlogin. -
distributedComputing()ships the model files of the objects ininputinstead of every file ofdMod2.outdir. A job of the same name that still runs where the new one would run is an error. The manifest is written before the start, and every job is registered untilpurge(local = TRUE);dMod2Jobs()lists them, andrecoverfinds a job from any working directory. The job haswait()andlog(i), andget()joins parlists withc()and data frames withrbind(), leaving out tasks without a result with a warning. -
rbind()of parameter frames keeps their attributes and joins the parts of a profile, such as its two sides, into one.confint()sorts the right side of a profile as it does the left. -
Breaking:
trust()takesboundaryas an argument instead ofstepControl$boundary.trustL1()takes its settings astrust()does:ftol,mtol,gtol,xtolandrminintolControl,thetaMax(wastheta.max) instepControl, andparupper,parlowerandboundaryafteriterlim; a flat name is refused with the list it moved to, anditerlimdefaults to 1000. -
The
check()of adistributedComputing()job counts a task whose process remains as a zombie as finished, as in a container whose init does not reap. -
distributedComputing(slurm = TRUE)reads the job ID ofsbatchalso when the build on the login node prints before it. -
plotCombined(aesthetics = )replaces a default aesthetic also under its other spelling, such ascolourforcolor. -
wide2long()of a matrix builds its data frame without a row-name warning.
-
Breaking:
steadyStates()has one solver, the former version"1.4": every unknown is solved from a combination of balances, so every result is a ratio of positive sums. The versions"1.0"to"1.3"and the argumentsversion,sparsifyLevel,outputFormat,positiveandbranchesare gone;testSteadyis"fast"or"skip". -
Breaking: the options are named
dMod2.*:dMod2.cores,dMod2.outdir,dMod2.batch.checkanddMod2.compile.cmdlimit; thedMod.*names are no longer read.dMod.quietanddMod.sym.<name>are gone:suppressPackageStartupMessages()silences the startup message, andSYMIDENT_<NAME>sets a switch of symident.?dMod2lists the options and environment variables. -
The solver options of
Xs()andXf()on thecppDEandSundialsbackends takemaxOrder, an upper bound on the order of the multistep methods, andstald, stability limit detection for"bdf". -
Breaking:
distributedComputing()is rebuilt and takes overrunbg();runbg()andprofileParsPerNode()are removed. The code runs once per task:tasksis a number, a list or a data frame, and the code reads its element or row astaskand its position astaskId, and the cores it may use ascores. Withoutslurmthe tasks run as detached background processes on"localhost", without ssh, or on machines reached by ssh, one after the other per machine and each with the machine's cores: thosemachinenames, as inc(server1 = 20, server2 = 10), or its free cores, and never more than are free. Withslurm = TRUEthey form an array job with the coresmachinenames, submitted through the job script template of the login node, with{{ }}fields.inputis taken from the calling environment, every task gets the seedseed + taskId,recoverneeds only thejobname, and the job gainsterminate();check()returnsTRUEonce no task runs any more and warns about failed tasks.cores,nRep,varValues,partition,nodes,sshPasswd,customFolders,returnAll,resetSeedsandpurgeLocalare gone.demo("BA_transport")sends its multistarts this way, andexamples/distributedComputing.mddescribes it. -
slurmTemplate()sets up the job script template of a SLURM cluster once per login node, kept across sessions, from a job script of the cluster or from the one dMod2 ships for the bwForCluster Helix. -
mstrust()draws every random start, those of retries included, before it spreads the fits overcores, soset.seed()reproduces a run whatevercores. -
unique()on a parframe returns the rows with distinct values sorted by value, also for a parframe not sorted by value; rows are selected withsubset()first. dMod2 needs ggplot2 3.5.0. -
rref()no longer takesverbose, which it did not use. -
detectFreeCores()names its values by machine and reaches"localhost"without ssh. -
symmetryDetection(unknownInputs = )declares unknown functions of time; an exact direction lists the change of each input asinputChange. Needs symident 0.12.0. -
Loading dMod2 no longer requires symident in the Python environment; symmetryDetection() and symmetryReduction() require it on first use, so a symident that cannot be resolved leaves the code generation of models intact.
-
Breaking:
getElements(),blockdiagSymb(),getLinVars(),long2wide(),lsdMod(),compare(),forcingsSymb()anddot()are removed, as is thecharactermethod ofreduceReplicates();reduceReplicates(read.csv(file))replaces it, and openxlsx and purrr are no longer suggested.progressBar()is no longer exported. -
Breaking: a dMod function refuses a named argument it does not take, more unnamed arguments than it has inputs and a missing input, so
x(times, pars, condition = "C1")is an error instead of a prediction for all conditions.Xs()andXf()refuse arguments they do not take. -
Breaking: the subset expressions of
plotPars(),plotValues()andplotArray()come right after the data,tol,direction,covtableandnSimafter them;covtabledefaults toNULL, and the subset applies with and without it.plotArray()drawsnSimtrajectories evenly spaced from the optimum, at most one per profile point, grouped by parameter set. -
Xd()applies the chain rule of a parameter transformation and differentiates by the free parameters only; a missing grid parameter is an error. -
predict()on a parframe labels each prediction with its own row, also for a parframe not sorted by value, and adds a condition grid with one covariate;plotResiduals()uses the row of each fit. -
addReaction()reads any numeric coefficient, such as3/2*A,(2)*Aor2*3*A, and refuses a term that is nocoefficient*species, a reaction without species and a state namedRateorDescription.getReactions()writes the coefficients as they are.subset()of an eqnlist takes%in%on a vector as base R does. -
rref()sets entries belowtolto zero and returns the pivots as an integer vector;nullZ()passestolon and returns the identity for a matrix without pivots. Conserved quantities have no rounding residues. -
lloq = NAin data means no limit of quantification. -
as.datalist()takes data without a condition or covariate column, as one condition"1", tibbles and data.tables. Renaming a datalist renames itsconditioncolumn,subset()of a datalist keeps the condition grid, andcovariates()of a data frame returns a data frame for any number of covariates. -
reduceReplicates()names conditions by the levels of factor columns. -
profile()checkswhichPar,stepControlandalgoControl, names a failed profile with its error, reports an objective error that ended a side and singular Hessians once per profile, and returns the optimum alone when no step is accepted. Each path step solves the linear system directly and takes the pseudoinverse only for an ill-conditioned one. -
confint()takes a profile without adatacolumn throughval.column = "value". -
mstrust()takes an anonymoussamplefunor one with a long name, passes the arguments of an optimiser other thantrust(), refusesfits < 1, and writes one log entry per fit also from PSOCK workers;load.parlist()reads the fits in the order of their numbers.mstrust()andprofile()share one parallel backend. -
A scaled objective
a %.*% objkeeps an absent gradient or Hessian absent and needs a finite scalara. -
On the
deSolvebackend,Xs()andXf()return the requested times only, warn when the solver stops early, name missing parameters, andXs()checksnamesand refuses an odemodel without sensitivities. -
[on a parframe keeps a single column as a parframe and selects columns by a single index, as for a data frame;as.parvec()of a parframe checksindex;parframe()takes as parameters the columns that are neither meta names nor objective attributes.summary()andplot()of a parlist skip aborted fits,summary()also takes an empty parlist and fits without convergence flag, andc()of parlists concatenates. -
eqnvec(),parvec()andprdlist()name unnamed entries by their position; a parvec withNAprints;wide2long()refuses duplicated names. -
plotPaths()andplotPathsMulti()build their paths in one pass;plotPathsMulti()normalizes each path to the range -1 to 1 and colors thenParslargest paths, the others in "Others".plotPars()detects steps on a frame in any order and names the columnstep. A layer added toplotCombined()uses its data.coordTransform()propagates errors only wheresigmais given.plotFluxes()takes any number of fluxes. -
wide2long(),lbind()andcombine()bind through data.table, which speeds up the data frames of predictions with many conditions. -
write.eqnlist()andas.eventlist()handle empty lists;is.eqnlist()checks that the columns of the stoichiometric matrix are the states. -
Pexpl()marks an output as fixed only when it depends on no varying input, andPimpl()keeps the rows of every solved state, so an output whose first derivative vanishes at the evaluation point keeps its second derivative.Pimpl()treats a species that no reaction changes, such as a catalyst, as a constant state, not as zero. -
In a sum of objectives,
datapointL2()finds the predictions of everynormL2()term, also behind a prior or a term thatconditionsskips.obj + NULLreturnsobj; adding anything else that is not an objfn is an error. -
constraintL2()and the other priors apply a namedsigmato the parameter it names and refuse asigmawithout an entry for a parameter ofmu.trustL1(minimize = FALSE)subtracts the penalty from the objective it maximises.evalConditionResidual()passesfixedparameters to the error model also without sensitivities. -
odemodel(deriv = FALSE)builds no derivatives, whateverderivMode.define()andinsert()refuseconditionMatchfor a trafo without a condition tree. -
eqnlist()takescompartmentswithoutcompartmentOf, with every state indefaultComp, and refusescompartmentOfwithout its compartments.c()of eqnlists keeps the layout of an eqnlist without reactions and the custom totals, andsubset()returns an empty eqnlist when no reaction matches.log10Transform()keeps names that sympy reads as constants, such asEandI.format.eqnvec()keeps powers and numbers such as1e-3whole.as.eqnvec()refuses empty and duplicated names. -
lbind()andas.data.frame()of a datalist fill a column that some conditions lack withNA.c()of datalists returns classlistas well and refuses a condition named twice.as.datalist()refuses covariate values whose condition names coincide. -
subset()of a parframe and the subsets ofplotPars(),plotValues()andplotProfile()take several expressions and see the frame of the caller.as.parframe()of a parlist fills a parameter that some fits lack withNA, andc()appends a single fit to a parlist.plotProfile(parlist = )draws the profiled parameters only,plotProfilesAndPaths()draws the total profile,plot()of a prdframe takes a data frame, andcoordTransform()keeps the order of the rows. -
modelname<-assigns a vector of names condition by condition and takesconditions;getParameters(conditions = )of a function for every condition returns its parameters. -
The code that dMod 1 contributors other than the authors wrote is reimplemented from its documented behaviour: the print methods of eqnvec, eqnlist, objlist and parvec,
summary(),plot()andas.parframe()of a parlist,fitErrorModel(),reduceReplicates(), the argument routing and fit log ofmstrust(),Xd(),define(),insert(),conservedQuantities(), the datalist methods and the prediction, data, flux, residual and array plots. With the oldsteadyStates()versions gone,DESCRIPTIONlists the two authors only. -
Xd()takes a state with a single grid point as constant;plotFluxes()labels an unnamed flux by its expression also when others are named;eqnlist()takes a plain vector as a one-column stoichiometric matrix;reduceReplicates()groups by the exact times instead of times rounded to seven digits. -
In a nested sum of objectives, a term without conditions is evaluated for any requested conditions, as in a flat sum.
-
summary()of a parlist heads its table "Best and worst fit"; the messages ofas.parvec()on an unconverged fit,as.eqnvec()on unparsable equations andsteadyStates()on namedgivenCQsare reworded.
loadDLL()also loads the shared objects of a model's auxiliary sources, such as the influx function ofPimpl(), somstrust()andprofile()with a condition axis run on PSOCK workers also for models compiled source by source.trust()andtrustL1()append the objective's error message toparinit not feasible: objfun failed.define(),insert()andrepar()ontrafo = NULLmap each left-hand side to its right-hand side, not to itself.
controls()on an objective also reaches the functions it evaluates: the prediction and error model of eachnormL2()term and the parameter transformation of anobjfn * parfn. Solver tolerances and other controls can be changed through the objective alone.conditionrestricts the lookup as for anfn.- An unknown entry of
options,optionsSensoroptionsReverseofXs()andXf()is an error, also when set later withcontrols(). Y()names a state the prediction does not return, e.g. after the controlnamesofXs()left it out.- A failed batch solve names the conditions that did not complete.
plotValues()shows unconverged fits in red with a legend entry of their own instead of a caption.
plotValues()colors converged fits by their iterations on the dMod ramp and draws unconverged ones in red, with a caption only when there are any. WithpseudoLog = TRUE, a range within one decade gets breaks at 1-2-5 or finer.- The demo
KokMSB2020shows the model scheme of the publication and the statistics of 200 fits, which ship ininst/extdata/kok/.
symmetryDetection()andsymmetryReduction()need symident 0.11.0. A direction checked without an error bound has level"checked", andinfo$assumptionslists the declared domain and the declared starts of the switches, with the error bound.symmetryReduction()names the claims of a section that are not proven underblocks[[i]]$unproven, with level"open", and listsassumptions; the level"assumed",blocks[[i]]$assumptionsand the argumentverboseare removed. A section on a face holdsfaceWhyandfaceLevel; a resting state checked only at random points of the face makes the block"checked".- In the blocks of
symmetryReduction(),newParameters,solvedDomainandchoicesreplacesurvivorMeaning,carrierDomainandadmissible. trust()takesphases, a list of phases, each with its ownhessianMethod,iterlim,sweepandqnControl, run in turn and cycled through for up toswitcheshandovers. A phase that uses up its owniterlimhands over like a soft stop.fallbackControlis removed:fallbackControl = list(hessianMethod = "sr1", limit = 3, iterlim = 200)isphases = list(list(hessianMethod = "gn"), list(hessianMethod = "sr1", iterlim = 200)), switches = 3. A run withoutphasesis unchanged.- The result of
trust()hasphaseIterations, the iterations of each phase run in run order, and withphasesthe sweeps of the phases asphaseSweepin place offallbackSweep.as.parframe()adds a columniterations.<phase>per phase. - Each run of a quasi-Newton phase has its own allowance of reseeds by
qnControl$hessianReseed = "stall";nReseedcounts over the run. mstrust(cores = c(fits = , conditions = ))runs both axes as given, each capped at the cores available, and records that layout inattr(, "run")$cores.mstrust()andprofile()stop the cluster of a conditions axis and restoreoptions(dMod.cores)on exit, on every platform.- The cores available are those of the R process,
parallelly::availableCores(), which follows the CPU affinity and a SLURM or PBS allocation. dMod2 imports parallelly. - A sum of objective functions keeps the
compileInfoof its terms, soloadDLL()on it searches the directories its sources were compiled in.
trustL1()dropsparscale, astrust()did in 0.12.2. A scaling of the parameters belongs inP().
trust()gainsfallbackControl, the settings of the fallback phase:hessianMethod(default"none"),limit,iterlim(per fallback phase),sweepandqnControl. Withsweep = NULL, the default, the fallback phase evaluates the objective backwards when the chain can.hessianFallback,fallbackLimitandparscaleare removed.trust()defaults toiterlim = 1000, counted over all phases.- The deprecated argument names and aliases are removed, among them
attr.name,opt.BLOQ,obj.attributes,sort.by,split.by,keep.covariates,rank.tol,attach.input,keep.root,optionsOde,start1stfromCenter,fterm,mterm,stepControl$theta.max, thecontrolsPTCentriesrtolandatol, thecoresofnormL2()andreadPetabYaml(),readPetabTables().repar(expr, trafo)is an error. mstrust()dropshessianfrom every fit and recordselapsed, the wall time of each fit, and the attribute"run"with the start time, total wall time, cores, optimiser, sampler, their arguments and the package versions.as.parframe()adds the columnelapsedand keeps"run".vcov()takes an evaluation of an objective function,vcov(obj(pars)), withfixedfor parameters held at a bound.vcov()on a fit without Hessian is an error instead of a matrix ofNA.- The
.Rfiles ofmstrust(cautiousMode = TRUE)restore the fits exactly.
- Needs symident 0.10.1.
- When symident needs msolve and finds none,
symmetryDetection()andsymmetryReduction()ask in an interactive session whether to build it into the symident cache, on Windows with Rtools, and then repeat the call.
- Needs cppDE 0.13.0 and symident 0.10.0.
getSymbols()andreplaceSymbols()are dMod2's own, no longer re-exported from cOde, which is now needed for the deSolve backend only.replaceSymbols()brackets a replacement that is not a single operand, replaces variables only, not function or argument names, and keeps the spacing.piis the number in equations, as in cppDE, andgetSymbols()does not list it.symmetryDetection()gainsregimes, the regimes an experiment visits in a model with switches, overall or per condition.- A direction of
symmetryDetection()has alevel("exact","verified"or"open") in place ofcertified, andstateComponents.infogainsrankStatus,proof,bound,gapCertified,limits,errorsandbranches. reconstControl()gainshomotopyDegree,degreeCapMax,exactProof,exactWorkCapandseed;NULL, the default, keeps the value of symident.symmetryReduction()stops at a direction of level"open"unlessallowOpen = TRUE. It gainsworkCap,timeoutdefaults toInf, andcheck = TRUEruns the detection again on the reduced model. The result gainslevel,limitsandchecks, and each block itslevel.plotValues()draws the objective values on a linear scale, not their distance to the best fit on a pseudo-log scale;pseudoLog = TRUEswitches to a pseudo-log10 scale.- The data set
kokand the demoKokMSB2020. The Bachmann, Boehm, bile acid and Hodgkin-Huxley scripts are demos as well,demo(package = "dMod2")lists them; the tutorials left ininst/examplesfrom dMod and the R2C templates ininst/templatesare removed, and the gradient comparisons moved tobench/gradient. bench/andbenchmarks/are one directory,bench/, by topic:harness/(the PEtab benchmark suite),gradient/,hessian/,symmetry/,multipleShooting/andnumSteady/.
- symident and python-libsbml are found in
~/.virtualenvs/r-reticulatetoo, installed on first use. summary(x, detailed = TRUE)on asymmetryReduction()result explains each section. Needs symident 0.9.3.
- New
trustL1(): the trust region oftrust()with an L1 penalty on chosen parameters around reference valuesmu, one-sided withone.sided.gatekeeps a parameter on or above its kink, so it can reach the reference exactly, andfusepenalises the pairwise differences within blocks, whose members move as one while equal. The penalised step runs in a C++ kernel.
symmetryDetection()andsymmetryReduction()install symident from PyPI on first use; underRETICULATE_PYTHON,pip install symident.
distributedComputing()runs zstd as a pipe stage of its own and moves the uploaded files one by one. macOS readstar -Ias a file list and itsmvhas no-t, so the upload and the collection of results failed there.
- Breaking:
Xs()on cppDE takesrefineandgradtolinoptionsReverse;flooris gone, and so is the weighting of the backward grid by the previous evaluation's adjoint.refine = TRUEholds each step of the backward sweep to the error test of the CVODES backward problem underabstolandreltol, throughcppDE::adjointControl(refine = TRUE), andgradtolis then the absolute tolerance on each step's share of the gradient. The reverse value pass takes the grid of a value solve. - New
symmetryDetection(),reconstControl()andsymmetryReduction(): structural non-identifiabilities of an ODE model, as generators of its symmetries, and the parameter transformation that removes them. The computation runs in the Python package symident through reticulate, which installs it on first use from https://github.com/dModverse/symident; R arguments and fields are camelCase,print()andsummary()are the reports of symident. symident's switches are read fromSYMIDENT_*and from the optionsdMod.sym.*. msolve, for coupled steady states, comes fromsymident.install_msolve()orSYMIDENT_MSOLVE. - reticulate no longer provisions symengine.
Xs()on cppDE and Sundials stops on a failed solve (onFailure = "stop"), asXf()does; a fit takes the error as a rejected step.Xs()on cppDE:optionsSens = list(sensErrCon = FALSE)takes the sensitivities out of the error test.odemodel(backend = "Sundials", derivMode = "reverse")takes events: the reverse objective takes the adjoint and the gradient through each jump.- Reverse objectives run their solves through prepared batch handles, evaluate the error model batched and keep the data-to-prediction indices between evaluations. Under second order a transformation that no tangent reaches skips its second derivatives.
normL2()with one condition evaluates the error model without derivatives when none are requested.importPEtab()gainssparse, which pins the sparse (KLU) or dense linear solver of the cppDE and Sundials backends.benchmarks/, outside the built package: value, forward and reverse gradients of PEtab benchmark problems by tier (run-benchmarks.R --tier tiny|medium|full), with the scripts ofinst/benchmarks/underbenchmarks/scripts/.- A sum of objectives keeps the
sweepattribute its terms agree on, so a reverse objective with a prior says it was evaluated in reverse. constraintL1()andconstraintL2()are S3 generics; the default methods are the priors on named parameters.- The option
dMod.outdirsets whereodemodel(),P(),Pimpl(),Y()and the PEtab import write generated sources and shared objects; unset, it is the working directory as before. getSymbols()andreplaceSymbols()are reexported from cOde.- Breaking:
Xs()takes its solver options asoptions, which applies to every solve, andoptionsSens, which now holds only the entries that overrideoptionsfor the solves with sensitivities.Xf()andimportPEtab()takeoptionsas well, and the controls are namedoptions,optionsSensandoptionsReverse.optionsOdeis still accepted as an argument and bycontrols(), with a deprecation warning. On deSolve the default methods stay"lsoda"and, for the solves with sensitivities,"lsodes". Xs()on cppDE warns about an unknown entry ofoptionsReverse, given to the constructor or set throughcontrols().- The solver option for the limit of consecutive rejected steps is
maxattempts, as in cppDE. - Fix:
constraintL2() * P()takesfixed, and with ittrust(fixed = )andprofile()on such priors. - Fix:
Xt()returns its sensitivities by the outer parameters, sog * Xt() * phas derivatives andnormL2()on it takesfixed. - Fix:
Y()differentiates by its own parameters when no transformation precedes the prediction:g * xlacked the derivatives by observation and error parameters that are not parameters ofx. - Fix:
normL2()takes an error model whose parameters are all fixed. - Fix:
as.datalist()stores asigmaorlloqcolumn ofNAalone as numeric. - Fix:
mstrust(cautiousMode = TRUE)removes the.Rdafits of its own run rather than files of that name in the working directory. - Breaking: arguments in dot case are camelCase.
Y()andPexpl()takeattachInput(wasattach.input),Pimpl()keepRoot(waskeep.root); both are also the names of the controls. - Breaking:
normL2(),datapointL2()and theconstraint*()priors takeattrName(wasattr.name), also as control ofdatapointL2();normL2()andevalConditionResidual()takeoptBLOQ(wasopt.BLOQ). - Breaking:
datapointL2()takesparameter(wasvalue), the name of the parameter holding the data value. - Breaking:
as.datalist()takessplitByandkeepCovariates(weresplit.byandkeep.covariates),as.parframe()sortBy(wassort.by),parframe()objAttributes(wasobj.attributes; the attribute keeps its name). - Breaking:
reml()andremlLeverage()takerankTol(wasrank.tol). - Breaking:
trust(stepControl = )takesthetaMax(wastheta.max);mstrust()takesstartFromCenter(wasstart1stfromCenter). - Breaking:
distributedComputing()takesmemPerCore,sshPasswd,varValues,nRep,purgeLocalandcustomFolders(weremem_per_core,ssh_passwd,var_values,no_rep,purge_localandcustom_folders); itspurge()takespurgeLocal. - Breaking:
plotPathsMulti()andplotProfilesAndPaths()takewhichParandnPars(werewhichparsandnpars),plotProfilesAndPaths()ncol(wasncols),plotArray()nSim(wasnsimus),plotFluxes()legendTitle(wasnameFlux),attrs()which(wasatr) andprofileParsPerNode()parsPerNode(wasfits_per_node). - Breaking:
reconstControl()takesminSupportCandCap,perPrimeCapandperPrimeMinPrimes(wereminsupportCandCap,perprimeCapandperprimeMinPrimes). - Breaking:
readPEtabYaml()andreadPEtabTables()replacereadPetabYaml()andreadPetabTables(). An imported problem prints as<PEtab problem ...>. - Breaking:
repar(trafo, expr)takes the transformation first, asdefine()andinsert()do. The orderrepar(expr, trafo)is recognised by an expression in first place. - Breaking:
profile()andvcov()are methods of the generics of stats, which dMod2 no longer masks.profile()dispatches on objective functions and plain R functions, first argumentfitted(wasobjfun).trust()returns a list of classtrustfit, kept bymstrust(), on whichvcov()dispatches, first argumentobject(wasfit); a fit stored before needsclass(fit) <- c("trustfit", "list"). - Breaking: the operands of
+and*on dMod functions, objectives, objlists and datalists are namede1ande2. - The old names above are accepted with a deprecation warning; giving the old and the new name together is an error.
- Breaking:
distributedComputing()submits by default (recover = FALSE). Theget()ofdistributedComputing()andrunbg()returns the results instead of assigningcluster_resultor.runbgOutputin the global environment;runbg(wait = TRUE)returns them too. The former behaviour iscluster_result <- job$get(). - Breaking:
steadyStates()writes nothing unlessfileis given; the model for the solver goes to a temporary file. Its unused argumentratesis deprecated and ignored, andverbosedefaults toFALSE. - Breaking:
fitErrorModel()plots only withplotting = TRUE;resetWarmStarts()andsymmetryDetection()print only withverbose = TRUE. - Breaking:
symmetryDetection()requiresfandg. ggopen()opens the PDF with the viewer of the platform by default:openon macOS,xdg-openon Linux,shell.exec()on Windows.trust(parscale = )is deprecated and warns; it still scales the trust region. Fit on log scale instead.msParframe()leaves the global random number generator as it was.mstrust()takessamplefunas a function as well as its name.rref()returns a list namedrrefandpivots; its unused argumentfractionsis deprecated and ignored.- Breaking:
Pimpl(controlsPTC)takes the tolerances asreltolandabstol, ascppDE::ptc()does;rtolandatolare deprecated. compile()restoresPKG_*flags set before the call without an error, and leaves nosymbols.rdsin the working directory underR CMD check.reconstControl(homotopy = )switches the last reconstruction route of symident, which fits directions still open underequilibrate = TRUEalong lines in parameter space; on by default.- Needs cppDE 0.12.0.
- The help pages of the constructors (
odemodel(),Xs(),Xf(),Xd(),Xt(),Y(),P(),Pexpl(),Pimpl(),normL2(),datapointL2(), theconstraint*()priors) and of the function classes are rewritten, with examples forXs(),P(),Pexpl()andPimpl(). Xs()andXf()on cppDE and Sundials return the forcings after the states, so observables can contain them, and take forcing names as factor. A forcing holds its first and last value outside its points and may be a single point.Xf()on cppDE and Sundials starts states missing fromparsat 0.datapointL2()returns anobjlist.constraintL2()stops on asigmathat mixes numbers and parameter names.normL2()names the data conditionsxdoes not have.plotValues()puts its breaks at 0 and at the decades from 10 on and no longer fixes the limits of the y axis.- Needs cppDE 0.11.3.
runbg()anddistributedComputing()withcompile = TRUEtake SUNDIALS' include path, and MPI's if that SUNDIALS was built with it, from the cppDE installed on the remote machine instead of the submitting one. A model using SUNDIALS' LAPACK dense solver links it there too, and the build stops with a message if the remote SUNDIALS lacks it.
Pimpl()draws its random starts uniformly overstartRange = c(-5, 5)onstartScale, log10 by default and linear ifpositive = FALSE. Each warm-start cache, i.e. each condition in each fit ofmstrust(), draws from its own stream; the global RNG is read, not advanced.maxitdefaults toceiling(70 * log(n + 1))fornstates (was 400).- Needs cppDE 0.11.2:
ptc()scales a row by its largest partial derivative, so a row whose own rate is tiny next to its partner's no longer stalls it, and it builds on macOS.
normL2()losest0.timesis a vector for all conditions or a list named by condition, and a condition's prediction starts at the first time of its grid, its data times and itstimes:times = list(C1 = -10, C2 = 5)is whatt0 = c(C1 = -10, C2 = 5)was. Withouttimesthe grid starts at the first data point, which matters only for a model built withodemodel(..., includeTimeZero = FALSE); with the default the grid holds 0. PEtab start times go throughtimes, in one objective for all conditions.- Multiple shooting.
normL2(data, g * x * p, multipleShootingControl = TRUE), or a list of settings, makestrust()andmstrust()fit by Bock's generalised Gauss-Newton method. The time axis is cut into segments with node values of their own; continuity is a constraint, linearised and eliminated at every iterate (condensing), so the subproblem keeps the single-shooting size and box bounds work as before. At convergence the result is a stationary point of the single-shooting objective, which the objective still is when called as a function. The ODE model needsodemodel(..., includeTimeZero = FALSE);profile()re-optimises by single shooting. - The trust region bounds parameters and node values together. Steps are
accepted by a filter or an l2 merit function (
stepControl$acceptance), with a second-order correction (stepControl$soc) and a gap tolerance (tolControl$ctol). stepControl$anneal(default) first penalises the gaps,f + w ||c||^2withwrising tenfold per stage, then hands over to the exact method. A stage ends once a step gains less than a thousandth; the annealing ends early where the penalty dominates and a stage no longer closes the gaps.nodes = "auto"(default) starts from ten segments per condition, or one per half oscillation of the data where that is more, and cuts a segment where its propagation matrix grows by more thangrowth, where it misses its data by more thanmisfitand by twice the median of its condition, or where its solve fails. A misfit everywhere is left to the parameters, and each half of a cut keepsminPointsdata points, by default one more than the states. Unobserved node values start from a run synchronised to the data.chartsputs node values on a log10, linear or angle scale per state; the gap of an angle is taken modulo 2 pi. The scale of an unobserved linear state is the largest range it covers within one segment.hessianMethod = "bfgs"or"sr1"keeps one quasi-Newton block per segment.- The cost of an evaluation is linear in the number of segments: segment
parameters are built per condition, the annealing model is assembled in one
pass, and the values at new nodes come from one batched solve.
fit$multipleShooting$cutscounts the cuts by growth, misfit and failed solve. c()of one parameter vector, alone or besideNULLor an empty one, as every prediction builds itsparameters, skips the list machinery.- Methods of the Freiburg group (Bock; Horbelt, Timmer and Voss 2002; Peifer
and Timmer 2007; Voss, Timmer and Kurths 2004), off by default:
stepControl$acceptance = "natural"(damped Gauss-Newton on the natural level function, needs a fixed sigma),stepControl$twoPhase(parameters first, then nodes),stepControl$regularise,stepControl$restore(restoration phase of the filter),init = "spline",minPoints,nodes = "transitions"(a node just before each fast change of the data) andbreaks(nodes without continuity). - Spiking data need tight integrator tolerances (1e-10); at 1e-8 the objective is too noisy far from the optimum and the trust region collapses.
- The objective returns the weighted residuals of each segment and their
Jacobian on request (
residuals = TRUE). A trial point after a rejected step is evaluated without sensitivities. - The reverse walk takes several first-order seeds at once, which
Xs()answers in one backward solve. Checkpoints serve any number of backward sweeps. - New benches
bench/multipleShooting_preBotC.Randbench/multipleShooting_macrospin.R.
Pimpl()runs a multistart when the warm starts and the initial guess fail:controlsPTC(nStarts = 20, startSd = 2, seed = 1), log-normal around the guess, from a fixed seed and without touching the global RNG. If every start fails it is an error.plotValues()draws the objective values minus the best one on a pseudo-log10 scale, linear below 1;"data"holds the distance indelta.
exportPEtab()reads the parameter scale offp: a parameter entering only asexp(X)is written withparameterScale = log, one entering only as10^(X)orexp10(X)withlog10, with linear nominal values and bounds five decades around them.parameterScale = NULLis the new default. v2 writes the linear values.exportPEtab()no longer writes self-referencing files for a trafo likeA ~ exp(A). An outer parameter named like a state goes out asinit_<state>and the condition table sets the species to it; one named like an inner parameter keeps its name under a pure scale wrap and becomes<name>_outerotherwise. Symbolic initial values go to the condition table instead of SBML<initialAssignment>. An export that would still refer an assignment to its own target stops with an error.exportPEtabObject()to v1 keeps preequilibration conditions in the condition table, leaves condition targets and compartment sizes out of the parameter table, writes observables in a v2 noise formula as their formula, and refuses condition switches during a simulation, which v1 cannot hold. A v2 experiment names no condition where the table has none.exportPEtabObject()to v1 moves an expression that sets a species in every condition into an SBML initial assignment and warns about per-condition expressions, which v1 does not allow; start times after 0 stop the export. Measurements name their conditions by id, as the condition table does.exportPEtabObject()to v1 keeps the SBML constants of an imported problem in the SBML, as v2 does. Listed in parameters.tsv they lost their role in a preequilibration, where a constant at 0 idles its reaction, and the reimport of Isensee_JCB2018 equilibrated to another state.exportPEtabObject()declares condition targets and compartment sizes that nothing else states, with the SBML default recorded at import (a surface compartment of Lang_PLOSComputBiol2024 went out sized by itself).exportPEtabObject()writes a single measured sigma per observable as its noise formula instead of per-row noise parameters.exportPEtabObject()translates priors between the versions: v1 spells the log familieslogNormal/logLaplace, and a v1 prior on the parameter scale goes to v2 as the matching prior of the linear parameter.importPEtab()reads the v1 spellings and truncates a prior on the parameter scale at the bounds on that scale; it used the linear bounds (Schwen_PONE2014).importSbml()substitutes the initial assignment of a constant parameter wherever the parameter appears, with species at their initial values. It kept the SBML default, so a condition that set the source parameter did not reach the rates (Laske_PLOSComputBiol2019:ModelValue_80 := k_syn_P; Bertozzi_PNAS2020:beta_N := R0_*gamma_/N_).repar(),insert()anddefine()keep a trailing underscore in an identifier.gamma_becamegamma, so a PEtab condition that setgamma_,N_,I0_orR0_in Bertozzi_PNAS2020 changed nothing.Y()returns a NaN observable instead of stopping, andnormL2()stops only on a NaN at a data point. A ratio of states that all start at 0 is undefined at t0 alone, where Laske_PLOSComputBiol2019 has no data, and the objective could not be evaluated.importPEtab(derivMode = c("forward", "reverse"))builds the observation, error model and parameter transformation for the reverse sweep as well; only the ODE model had it, soobj(sweep = "reverse")failed.exportSbml()declares species that enter a rate without being consumed or produced as modifiers, as SBML requires, and stops on a compartment sized by its own symbol without a value.
Pimpl()solves steady states by pseudo-transient continuation: implicit Euler steps under local error control that turn into Newton steps near the root. It follows the flow to a stable steady state and checks the eigenvalues on the conservation manifold;flow = FALSEreturns any regular root of plain equations. The iteration runs in C++,cppDE::ptc().Pimpl()keeps every species and adds the conservation rowsC x = Tinstead of eliminating pivot species. Iterates stay positive and on the manifold. States without influx at the given parameter values are 0.Pimpl()starts deterministically: the nearest kept root of the condition corrected by its sensitivities, then the initial guess. Conditions with identical inputs are solved once. A failed solve is an error.Pimpl():controlsPTCreplacescontrolsMSandcontrolsNleqslv,expressInTotalsis gone. nleqslv is no longer a dependency.Pequil()is removed.P()builds aneqnlistwithPimpl();method = "equilibrate"is gone.importPEtab()computes preequilibration withPimpl(). Events do not act during preequilibration (PEtab v2 case 0023).importPEtab(): species a preequilibration cannot move keep their initial values. They take part only in reactions idled by SBML constants at 0 or by species that start and stay at 0, and the steady-state equations leave them undetermined (Isensee_JCB2018: AC/pAC and PDE/pPDE).
steadyStates()1.3 and 1.4 detect structurally zero clusters that span compartments, such as a ligand in the medium bound by a receptor in the cell. The sink-cluster test runs in amounts, over reactions rather than the rows per volume ratio the backend reads. As in 0.8.3, a conserved moiety in the support of the test stays nonzero.
- Structurally zero states: a conserved moiety listed before a leaking cluster was declared zero together with it, since the sink-cluster LP returns their union. Only species that drain into a leaking reaction are zero now.
steadyStates(version = "1.4")prints the resolved expressions through sympy again, as 1.3 does: no redundant brackets from the substitution, nothing expanded.steadyStates(version = "1.4", solveQuadratic = TRUE): a state whose own balance, denominators cleared, is quadratic with positive production and quadratic consumption takes the unique positive root, tried before any rate constant. A root that would take a balance another unknown needs is dropped for the next attempt.steadyStates(version = "1.4", givenCQs = ): each conserved quantity keeps one of its states free, the first one whose balance can be spent elsewhere.customTotals()reach it as before.steadyStates(version = "1.4")never solves for a rate constant whose fluxes all contain aneglected symbol, which would divide by that symbol.steadyStates(version = "1.4")prints one line per attempt instead of the intermediate expressions, notes after its summary, and on failure the balances left.steadyStates()handles species named like sympy objects (Ci,E,S,Q,gamma).
steadyStates(version = "1.4"): a new core on the same interface. Balances are linear forms over the flux terms, and each unknown comes from a combination of balances a linear program finds, so every solution is a ratio of positive sums. Where no single unknown is left, one side of a balance shares its sum by new flux ratiosr_*. Needspositive = TRUE; volumes are never unknowns, andsolveQuadratic,branchesandgivenCQsdo not apply.steadyStates(verbose = ):TRUE(default) prints a few progress lines and the result,FALSEthe result only,"full"every step.steadyStates()1.3 solves linear balances directly instead of bysolve(), reads signs off the expression tree after a numeric screen, and searches a block balance by balance when no pivot combination is positive.steadyStates()solves a strongly connected block of states jointly where no single balance stays positive, e.g. receptors recycling through endosomes.steadyStates()solves every rate constant at most once. A symbol defined twice, by itself or in a cycle is an error.steadyStates()factors out common terms instead of cancelling the resolved expressions, which stalled large models.steadyStates(testSteady = "fast")checks mod p only. A residual it cannot evaluate is an error with its reason, not a fallback to the symbolic test.normL2(t0 = )takes one start per condition, named by condition.Pexpl()andY()buildderivMode = "forward"by default and accept"forward-reverse", which the reverse sweep withderiv2 = TRUEneeds.controls()reaches every option an object uses at run time, andcontrols<-changes it for the next call.Pimpl()(keep.root,controlsMS,controlsNleqslv) andPexpl()(attach.input) have controls now;Pequil()readskeep.rootandattach.inputfrom them, and its memo no longer answers a call made with other controls.Xs()andXf()merge the solver defaults at each call, so replacingoptionsOdekeeps the options it does not name. The accessor walksg * x * p, the summands of a sum (attributeterms) and the objective inside%.*%andobjfn * parfn;condition = NULLsets every condition, and an unknown name or condition is an error.mstrust()routes its arguments by the formals ofoptmethodinstead of those oftrust(), and no longer passescautiousModeto the objective.- Bug fixes.
normL2()no longer reads past a prediction the solver cut short, which gave a wrong value or a segfault; it is an error naming the time the solver reached, whichtrust()takes as a failed step.Pexpl(attach.input = TRUE)gives the inputs it passes through their forward derivatives, so the forward gradient and Hessian along them are no longer zero; inputs infixedare handed on as constants. The name clash check ofmstrust()measures the intersection. - Requires cppDE 0.10.2.
- Reverse mode.
obj(pars, sweep = "reverse")walksnormL2 -> Y -> Xs -> Pbackwards in one sweep, at a cost that does not grow with the number of parameters; withderiv2 = TRUEit returns the exact Hessian at a cost linear in it. Needsodemodel(..., derivMode = c("forward", "reverse")), or"forward-reverse"for second order;backend = "Sundials"uses the CVODES adjoint, andimportPEtab(backend = "Sundials")imports onto it. The reverse gradient belongs to the trajectory of the value solve and reuses its checkpoints;Xs(optionsReverse = list(gradtol = ))weights the backward step size. - Derivative arguments.
deriv,hessian,deriv2andsweepeach ask for one thing, andcurvatureis gone.derivModenames the direction at build time inodemodel(),Y(),Pexpl()andimportPEtab();reverse = TRUEand"dual"are errors."symbolic"is gone, so transformation and observation functions evaluate only aftercompile(). - Optimisation.
trust()asks per evaluation for a value, a gradient, a Gauss-Newton or an exact Hessian.qnControl$hessianInit = "exact"seeds a quasi-Newton run,hessianReseed = "stall"fetches a fresh Hessian when it stalls, andhessianMethod = "exact"is a Newton run. - Requires cppDE 0.10.0 and follows its names:
tangent,hessian,cotangentandcurvature, andcppFUN()forfunCpp(). distributedComputing()andrunbg()gainlibs, library paths put in front on the remote side.compile()reports the toolchain and reused objects throughmessage()rather thancat(), sosuppressMessages()silences it.- Bug fixes. A summed objective keeps the curvature of every term. An
objective that declines a Hessian no longer crashes
trust(). A backward solve without an answer is an error instead of a zero gradient. A PEtab prior honourshessian = FALSE. The batched backward path propagates more than one direction.compile()stays within the 8191-character command line of Windows. - New examples
inst/examples/example_ReverseAD.R,example_AdjointComparison.Randexample_BachmannReverse.R; the tests compile the models of each file jointly.
runbg(compile = TRUE)builds on Debian and Ubuntu machines again. The generated build script looked for R's headers in$(R RHOME)/include, which Debian does not have, so every source failed withR.h: No such file or directoryand the job never started. The script now asks the remote R forR.home("include"), ascompile()does locally.runbg(compile = TRUE)compiles in parallel again. It exportsOMP_NUM_THREADS=1before the build, andnprochonours that, so the build counted one core and compiled every source in turn. The count now ignores the OpenMP variables.runbg()also gainsbuildCoresandbuildBundle, asdistributedComputing()has.
- A prepared ODE batch handle no longer outlives the shared object it was
resolved from.
Xs()keyed its cache on shapes and labels only, so a workspace shipped to a cluster node called an entry point that was never built there and everymstrust()start failed withparinit not feasible. The cache now checks that the object is loaded,modelname<-drops it, and the rename loop inrunbg()anddistributedComputing()covers objectives. - Needs cppDE 0.9.5, where a Windows install detects OpenMP. Without it a batched solve stays serial there.
- SBML and PEtab import work again in a fresh session.
.dmod_libsbml_python()readreticulate::py_exe()before Python was initialised, where it falls back to discovery, misses the reticulate-managed environment and returns an empty string. Thesystem2()that followed then failed with status 127, which the resolver reported as a missingpython-libsbml, so a working environment looked like a broken one. It now initialises Python first and rejects an empty path. - The libsbml integration tests did not catch that, because
.libsbml_works()turns any error intoFALSEand every test behind it into a skip. A test now covers the interpreter resolution itself, outside that guard.
trust()reports the iterate the run actually reached. The best-iterate bookkeeping introduced with non-monotone acceptance took its snapshot at the head of an iteration, so a run ending onfvalue,preddiff,stepor the iteration limit reported the second to last accepted point.argument,value,gradient,hessianandatBoundwere affected, and with themvcov(),confint()and profiles. The objective value moved by less thanftolon the three soft stops, but the gradient could be orders of magnitude too large, and an iteration-limit stop lost a full step. Runs ending ongradient,stagnation,radiusorobjfunwere never affected.trust()takeshessianFallbackandfallbackLimit. The fallback takes over at the first value, model, step or stagnation stop instead of ending the run;fallbackLimitabove one alternates back to the primary source at the next such stop, re-seeding a Gauss-Newton phase from a fresh Hessian at the cost of one evaluation.hessianMethod = "hybrid"is gone, not deprecated: a run is a method and, optionally, a fallback, and no pair has a name of its own. It was"gn"withhessianFallback = "bfgs", andtrust()says so when it is passed.as.parframe()includesnSwitch, the number of Hessian source handovers, so a multi-start can be scored on them.trustL1()is gone. It had its own trust-region driver, shared no acceptance semantics withtrust()after the changes above, and had no caller left here: the L1 penalty, its clustering and the EM layer that use it live ondevel-EM, and it returns with them once that layer lands, aligned with the interface above.constraintL1()stays as the Laplace prior it always was.createExample(),exmpextr()andextractExamples()are gone. They served the pre-testthatinst/testslayout, andcreateExample()opened an editor on a template to write a unit test, which is maintainer tooling rather than package API.distributedComputing()andrunbg()attachdMod2on the remote machine and nothing else. They used to replicate every package attached in the submitting session, plus a hard-codedtidyverse, which made a job depend on the session that sent it and filled the node logs with failed loads. A job that needs another package should attach it in its own expression.trust()groups its arguments intotolControl,qnControlandstepControl. The flatftol,mtol,gtol,xtol,rmin,boundary,theta.max,hessianInit,qnMemoryandqnCautiousare gone, not deprecated;trust()andmstrust()name the control member a moved argument became.trustL1keeps its flat tolerances.trust()takesqnControl$qnRejected(SR1 also updates from a rejected trial point, Nocedal and Wright sec. 6.2; defaultTRUEwherever an SR1 phase can occur) andstepControl$nonmonotone(Zhang-Hager acceptance, default0, the monotone rule).
trust()ends a quasi-Newton run on the gradient.fvalue,preddiffandstepread the prediction of the quadratic model, which during a quasi-Newton phase describes the approximation being built rather than the iterate, so they no longer terminate such a run;gradient,stagnationand the radius still do. This follows the usual recommendation for quasi-Newton methods. On the Boehm model it raises the starts that reach the published optimum from 10 to 13 for"bfgs"and from 12 to 18 for"sr1", and lowers the evaluations spent per such start in both.trust()takes an interchangeable Hessian source throughhessianMethod."gn"(default) is the Gauss-NewtonJ^T Jas before;"bfgs"and"sr1"maintain a dense quasi-Newton update seeded from it;"hybrid"runs"gn"until it stagnates, then switches to"bfgs"once. The quasi-Newton phase consumes only the gradient. Reflective boundary only.- Objective functions take a call-time
hessianargument. Withhessian = FALSEtheJ^T Jcontraction never runs and the result has aNULLhessian.trust()uses this in the quasi-Newton phase, and it propagates through objective composition. trust()reportsnevaland, underblather, thehessianSourceper iteration, so a multi-start can be scored on gradient evaluations. These reachas.parframe()as columns.mstrust()drops the deprecatedstudynameargument. Usename.compile()builds model shared objects correctly on Windows. The per-source link and the combined-output object compile ransystem()with a2>&1token, but R'ssystem()on Windows has no shell: the token reachedR CMD SHLIBas the overridePKG_LIBS=2>&1and the compiler as an input file, so the object never built and the.dllnever linked against BLAS, LAPACK or the Sundials solver. Both now go throughsystem2().- Adds
inst/examples/example_Boehm_JProteomeRes2014.R, which builds the Boehm et al. (2014) STAT5 model and compares the Hessian sources over a multi-start. define()andinsert()resolve the values passed through...in the frame they were called from. Resolution reached the global environment only, so the same call that worked at top level failed with an "object not found" inside a function or a knitr chunk. Condition columns and.currentSymbolskeep precedence over the calling frame. The internal names ofinsert()are no longer visible to those values either;.currentTrafoand.currentSymbolsare the documented way to reach the branch being rewritten.plotValues()marks a converged fit with a circle and an unconverged one with a triangle, in every plot, and keeps both in the legend. The mapping followed the levels present in the data, so it could differ between two plots of the same kind.subset()on aparframeevaluates its condition against the columns of the frame and then in the frame it was written in. The second step reached the package namespace instead, so a condition naming a local variable failed from inside a function.- Adds the
Optimisationvignette on Hessian sources, seeding and the economics of a multi-start, shipped pre-rendered, withinst/benchmarks/bench_hessianSource.Rbehind its numbers. - Attaching the package reports that BLAS is pinned to one thread inside the
forked workers of
mstrust()andprofile(), so a serial BLAS there is not a surprise, and says so when no thread-control entry point was found and the deadlock is still reachable.options(dMod.quiet = TRUE)suppresses the line. The pin itself is cppDE's and needs no cooperation here; requires cppDE 0.9.4.
- PEtab import and export.
importPEtab()reads a v1 or v2 problem and returns the composedg * x * p, the error model, the objective with its priors and the published parameter vector.exportPEtabObject()writes one back, andexportPEtab()builds a problem from a hand-written dMod model. - SBML import and export through libsbml, reached over reticulate, so a problem
that needs neither stays Python-free. Identifiers that R or C++ reserve are
renamed throughout, including in the PEtab tables, and the MathML constants
pi,exponentialeandavogadroare folded to their value. normL2()returns gradient and Hessian in the order of the parameter vector it was called with. They used to follow the union of the per-condition sensitivity blocks, whiletrust()andoptim()read them positionally, so every fit optimised a permuted model.normL2()reports its sum of squares as achi2attribute. Adding objectives pools the terms sharing anattr.nameand splits the rest intochi2_<attr.name>.compile(output = )points the cOde models inside a function object at the batched shared object. Without that the deSolve backend looked for an entry point in a library that was never built.- Pre-equilibration takes its conserved moieties from the same reduced network
Pequil()uses, freezes a time dependent input at the start time, and leaves symbols no outer parameter reaches out of the sensitivity system. A failure now reports the solver error and the non-finite parameters. Y()accepts one observation function per condition and takescores.importPEtab()takescores,deriv,outdir,optionsOdeandoptionsSens, and derives a model name that steps aside for what is already loaded.as.data.frame()on a prediction joins the error model by observable rather than by row order, which used to hand each observable the wrong sigma.wide2long()on an unnamed list numbers its conditions instead of dropping theconditioncolumn, and forwardskeepandna.rm.print()on an objective result from aderiv = FALSEcall shows the value alone, where it used to fail on the absent Hessian.- New data set
bachmannwithinst/examples/example_BachmannMSB2011.Rbuilding the model: 25 states, 36 conditions, 113 estimated parameters. - All 31 cases of the PEtab v2 test suite reproduce the published likelihood and survive the export round trip. Of the 35 problems in the Benchmark-Models collection, 33 import, evaluate with derivatives, export and reimport.
compile(output = )no longer links into a shared object the process already holds. Overwriting one is not portable: Windows may keep the file handle and macOS may keep the image resident, so the reload would serve the old code. A loaded name is replaced by<name>_2with a warning, the rule cppDE's model constructors already apply to their own names.
- New soft constraints next to
constraintL2():constraintL1()(Laplace),constraintCauchy(),constraintGamma(),constraintExponential(),constraintChisq()andconstraintRayleigh(). Each is the-2 logdensity including its normalisation, unlikeconstraintL2(), which is the penalty form, and each applies the chain rule soconstraint * P()is exact. - The multivariate-normal path of
constraintL2()is gone, together with thepenaltySpecplumbing and theplotIndivs()/plotHistIndivs()generics. All of them need anomegaSpecor apenaltyspec, which only the NLME layer builds, so here they were unreachable. constraintExp2(), an unused box prior inherited from dMod 1.x, is removed.
normL2()takest0, the time at which initial values take effect. The time grid starts there rather than at 0, so a prediction may begin after 0.- The error model keeps its fixed parameters under
deriv = FALSE. They were derived from the sensitivity rows, which are absent then. Xt()reports its fixed parameters, asXs()does, so an error model can read them off the prediction.eqnlist()accepts a model without species, which is a parameter-only problem rather than an error.- The steady-state heuristic evaluates SBML's
piecewise()and declares a state structurally zero only when its rate vanishes at zero alone.
- Reloading a model's shared object no longer invalidates anything. cppDE
(>= 0.9.2) resolves entry points by name, so the batch handle cached by
Xs(), and every objective built on it, survives an unload and reload. This replaces the symbol-cache flush counter of 0.6.1, which only caught reloads that went throughloadDLL()orcompile(). loadDLL()also searches the directories the sources were generated in, so a model compiled into a temporary folder is found rather than silently skipped.normL2()attachescompileInfo, as*and+already did, soloadDLL()andcompile()reach the shared objects of a composed objective.
loadDLL()skips shared objects that are already loaded in the current process. Unloading them nulled the native symbol pointers held by the prediction, observation and parameter functions built from them, leaving every later call without a way to resolve them again.Xs()caches a prepared batch handle that holds such a pointer of its own, which no symbol-cache flush reached. It is now rebuilt whenever the cache is flushed.
- The core is complete: equations, compiled prediction, observation and parameter transformations, objectives, trust region optimisation and profile likelihood. The PEtab, symmetry, mixed effects and Bayesian layers live on their own development branches.
- Test suite over equations, compilation, prediction, objectives, the trust region and profiles.
plotProfile(),plotPaths()andplotValues()build on one long-format frame; theparlistoverlay is measured against the profile optimum.- dMod palettes,
theme_dMod()and the colour scales.
runbg()runs jobs in the background on a local or remote host.distributedComputing()spreads fits and profiles over a cluster.
reml()estimates the error model from the restricted likelihood, charging every data point its own leverage.remlLeverage()reports the hat values, their rank and the effective degrees of freedom per observable.
mstrust(cores = c(fits = , conditions = ))andprofile(cores = c(pars = , conditions = ))split the two parallel axes.profileThreshold()computes the chi-square or finite-sample F threshold thatprofile(),plotProfile()andconfint()share.
steadyStates()solves for symbolic steady states through sympy.
Pexpl(),Pimpl()andPequil()build explicit, implicit and equilibrated parameter transformations.define(),insert(),branch()andrepar()assemble transformations.- Warm starts are kept per condition, so a condition restarts from its own previous root.
Xs()hands every condition tocppDE::solveODEBatch()in one call.Y()andP(method = "explicit")evaluate all conditions in one call.coresis a call-time argument on every function object and defaults togetOption("dMod.cores", 1L).
- Fixed: composing a condition-less function with a multi-condition one kept only the first condition and returned an unnamed result.
- The loop over experimental conditions moved from outside a composed chain into its leaves, so a leaf sees every condition at once and can batch them.
- Predictions and objectives are identical across thread counts.
- Constraint and parameter-vector kernels in C++.
trust()andtrustL1()run on a C++ trust region kernel.
normL2()evaluates its data term in a C++ kernel.
- Residual kernel in C++, with error models and BLOQ handling (
M1,M3,M4NM,M4BEAL).
compile()links model sources into one shared object, with compile and link flags taken per file from the backend that produced it.- Objects are reused when source and command are unchanged, recorded in a
.dMod_objectsindex.
odemodel()generates C++ with first and second order sensitivities. BackendscppDE,SundialsanddeSolve.
- An
eventlistis declared onodemodel(), so the sensitivity equations are extended consistently.
eqnliststores compartments and volumes;assignCompartment()andsetCompartmentVolume()set them.conservedQuantities()andgetTotals()report the conserved moieties.
- Package entry point renamed to
dMod2.
- Symmetry detection and SBML import moved out of the core.
- Removed
dModFrame,dfoptimandnormIndiv.
- The per-author
tools*.Rfiles are merged intoutils.Randaccessors.R.
- Fork of dMod 0.5.