COPASI API  4.16.103
CNewtonMethod.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <algorithm>
16 #include <cmath>
17 
18 #include "copasi.h"
19 
20 #include "CSteadyStateMethod.h"
21 #include "CSteadyStateProblem.h"
22 #include "CSteadyStateTask.h"
23 
26 #include "model/CState.h"
27 #include "model/CModel.h"
28 #include "model/CCompartment.h"
31 //#include "trajectory/CTrajectoryMethod.h"
33 #include "utilities/utility.h"
35 
36 #include "lapack/lapackwrap.h"
37 #include "lapack/blaswrap.h"
38 
40  CSteadyStateMethod(CCopasiMethod::Newton, pParent),
41  mIpiv(NULL),
42  mpTrajectory(NULL),
43  mStartState()
44 {
46 }
47 
49  const CCopasiContainer * pParent):
50  CSteadyStateMethod(src, pParent),
51  mIpiv(NULL),
52  mpTrajectory(NULL),
53  mStartState()
54 {
56 }
57 
59 {cleanup();}
60 
62 {
63  CCopasiParameter *pParm;
64 
65  assertParameter("Use Newton", CCopasiParameter::BOOL, true);
66  assertParameter("Use Integration", CCopasiParameter::BOOL, true);
67  assertParameter("Use Back Integration", CCopasiParameter::BOOL, true);
68  assertParameter("Accept Negative Concentrations", CCopasiParameter::BOOL, false);
69  assertParameter("Iteration Limit", CCopasiParameter::UINT, (unsigned C_INT32) 50);
70  assertParameter("Maximum duration for forward integration", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1e9);
71  assertParameter("Maximum duration for backward integration", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1e6);
72  //assertParameter("Force additional Newton step", CCopasiParameter::BOOL, true);
73  //assertParameter("Keep Protocol", CCopasiParameter::BOOL, true);
74 
75  // Check whether we have a method with the old parameter names
76  if ((pParm = getParameter("Newton.UseNewton")) != NULL)
77  {
78  setValue("Use Newton", *pParm->getValue().pBOOL);
79  removeParameter("Newton.UseNewton");
80 
81  if ((pParm = getParameter("Newton.UseIntegration")) != NULL)
82  {
83  setValue("Use Integration", *pParm->getValue().pBOOL);
84  removeParameter("Newton.UseIntegration");
85  }
86 
87  if ((pParm = getParameter("Newton.UseBackIntegration")) != NULL)
88  {
89  setValue("Use Back Integration", *pParm->getValue().pBOOL);
90  removeParameter("Newton.UseBackIntegration");
91  }
92 
93  if ((pParm = getParameter("Newton.acceptNegativeConcentrations")) != NULL)
94  {
95  setValue("Accept Negative Concentrations", *pParm->getValue().pBOOL);
96  removeParameter("Newton.acceptNegativeConcentrations");
97  }
98 
99  if ((pParm = getParameter("Newton.IterationLimit")) != NULL)
100  {
101  setValue("Iteration Limit", *pParm->getValue().pUINT);
102  removeParameter("Newton.IterationLimit");
103  }
104 
105  removeParameter("Newton.LSODA.RelativeTolerance");
106  removeParameter("Newton.LSODA.AbsoluteTolerance");
107  removeParameter("Newton.LSODA.AdamsMaxOrder");
108  removeParameter("Newton.LSODA.BDFMaxOrder");
109  removeParameter("Newton.LSODA.MaxStepsInternal");
110  }
111 }
112 
114 {
117  return true;
118 }
119 
121 {
122  if (mIpiv) delete [] mIpiv; mIpiv = NULL;
123 
125 }
126 
127 void CNewtonMethod::load(CReadConfig & configBuffer,
129 {
130  if (configBuffer.getVersion() < "4.0")
131  {
132  C_FLOAT64 Dbl;
133  C_INT32 Int;
134  bool Bool;
135 
136  configBuffer.getVariable("SSStrategy", "C_INT32", &Int, CReadConfig::LOOP);
137 
138  switch (Int)
139  {
140  case 0:
141  setValue("Use Newton", true);
142  setValue("Use Integration", true);
143  setValue("Use Back Integration", false);
144  break;
145 
146  case 1:
147  setValue("Use Newton", false);
148  setValue("Use Integration", true);
149  setValue("Use Back Integration", false);
150  break;
151 
152  case 2:
153  setValue("Use Newton", true);
154  setValue("Use Integration", false);
155  setValue("Use Back Integration", false);
156  break;
157 
158  case 3:
159  setValue("Use Newton", false);
160  setValue("Use Integration", false);
161  setValue("Use Back Integration", true);
162  break;
163 
164  default:
165  fatalError();
166  break;
167  }
168 
169  configBuffer.getVariable("SSBackIntegration", "bool", &Bool);
170  setValue("Use Back Integration", Bool);
171 
172  configBuffer.getVariable("NewtonLimit", "C_INT32", &Int,
174  setValue("Iteration Limit", Int);
175 
176  configBuffer.getVariable("SSResoltion", "C_FLOAT64", &Dbl); //typo is necessary!!
177  //setValue("Steady State Resolution", Dbl);
178  //setValue("Derivation Resolution", Dbl);
179  //setValue("Stability Resolution", Dbl);
180  setValue("Resolution", Dbl);
181  }
182 }
183 
184 //**********************************************************************************
185 
187 {
188  C_FLOAT64 iterationFactor = forward ? 10.0 : 2.0;
189  C_FLOAT64 maxDuration = forward ? mMaxDurationForward : -mMaxDurationBackward;
190  //minimum duration is either hardcoded or equal to maximum duration, whichever is smaller.
191  C_FLOAT64 minDuration = forward ? (mMaxDurationForward < 1e-1 ? mMaxDurationForward : 1e-1)
192  : -(mMaxDurationBackward < 1e-2 ? mMaxDurationBackward : 1e-2);
193 
194  //progress bar
195  size_t hProcess;
196  unsigned C_INT32 Step = 0;
197  unsigned C_INT32 MaxSteps;
198  MaxSteps = (unsigned C_INT32) ceil(log(maxDuration / minDuration) / log(iterationFactor));
199 
200  std::string tmpstring = forward ? "forward integrating..." : "backward integrating...";
201 
202  if (mpCallBack)
203  hProcess = mpCallBack->addItem(tmpstring,
204  Step,
205  & MaxSteps);
206 
207  //setup trajectory
208  CTrajectoryProblem * pTrajectoryProblem = NULL;
209 
210  if (mpTrajectory)
211  {
213 
214  pTrajectoryProblem =
215  dynamic_cast<CTrajectoryProblem *>(mpTrajectory->getProblem());
216  assert(pTrajectoryProblem);
217  pTrajectoryProblem->setStepNumber(1);
218  }
219 
220  bool stepLimitReached = false;
221  C_FLOAT64 duration;
222 
223  for (duration = minDuration; fabs(duration) <= fabs(maxDuration); duration *= iterationFactor, Step++)
224  {
225  if (mpCallBack && !mpCallBack->progressItem(hProcess)) break;
226 
227  pTrajectoryProblem->setDuration(duration);
228 
229  try
230  {
231  // We must not use useInitialValues = true here as this will interfere with scan continuation.
233  stepLimitReached = !mpTrajectory->process(false); //single step
234  }
235 
236  catch (CCopasiException & /*Exception*/)
237  {
239 
240  if (mKeepProtocol)
241  mMethodLog << " Integration with duration " << duration << " failed (Exception).\n\n";
242 
243  break;
244  }
245 
247 
248  if (containsNaN())
249  {
250  if (mKeepProtocol)
251  mMethodLog << " Integration with duration " << duration << " failed (NaN).\n\n";
252 
253  break;
254  }
255 
256  if (!(mAcceptNegative || allPositive()))
257  {
258  if (mKeepProtocol)
259  mMethodLog << " Integration with duration " << duration
260  << " resulted in negative concentrations.\n\n";
261 
262  break;
263  }
264 
266  C_FLOAT64 value = targetFunction(mdxdt);
267 
268  if (isSteadyState(value))
269  {
270  if (mpCallBack) mpCallBack->finishItem(hProcess);
271 
272  if (mKeepProtocol)
273  mMethodLog << " Integration with duration " << duration
274  << ". Criterium matched by " << value << ".\n\n";
275 
276  return CNewtonMethod::found;
277  }
278  else
279  {
280  if (mKeepProtocol)
281  mMethodLog << " Integration with duration " << duration
282  << ". Criterium not matched by " << value << ".\n\n";
283  }
284 
285  if (mUseNewton)
286  {
287  if (mKeepProtocol) mMethodLog << " Try Newton's method from this starting point. \n";
288 
289  NewtonResultCode returnCode = processNewton();
290 
291  if (mKeepProtocol) mMethodLog << "\n";
292 
293  // mpParentTask->separate(COutputInterface::DURING);
294 
295  if (returnCode == CNewtonMethod::found)
296  {
297  if (mpCallBack) mpCallBack->finishItem(hProcess);
298 
299  return CNewtonMethod::found;
300  }
301  }
302 
303  if (stepLimitReached)
304  {
305  if (mKeepProtocol)
306  mMethodLog << " Integration with duration " << duration
307  << " reached internal step limit.\n";
308 
309  break;
310  }
311  }
312 
313  if (mpCallBack) mpCallBack->finishItem(hProcess);
314 
316 }
317 
318 //**********************************************************************************
319 
321 {
322  //clear log
323  mMethodLog.str("");
324 
325  if (mpCallBack)
326  mpCallBack->setName("performing steady state calculation...");
327 
330 
331  NewtonResultCode returnCode;
332 
333  // Newton
334  if (mUseNewton)
335  {
336  if (mKeepProtocol) mMethodLog << "Try Newton's method. \n";
337 
338  returnCode = processNewton();
339 
340  if (returnCode == CNewtonMethod::found)
341  return returnProcess(true);
342  }
343 
344  // forward integration
345  if (mUseIntegration)
346  {
347  if (mKeepProtocol) mMethodLog << "\nTry forward integration. \n";
348 
349  returnCode = doIntegration(true); //true means forward
350 
351  if (returnCode == CNewtonMethod::found)
352  return returnProcess(true);
353  }
354 
355  // backward integration
357  {
358  if (mKeepProtocol) mMethodLog << "\nTry backward integration. \n";
359 
360  returnCode = doIntegration(false); //false means backwards
361 
362  if (returnCode == CNewtonMethod::found)
363  return returnProcess(true);
364  }
365 
366  return returnProcess(false);
367 }
368 
369 //**************************************************************
370 
372 {
373  memcpy(mXold.array(), mpX, mDimension * sizeof(C_FLOAT64));
374 
375  // DebugFile << "Iteration: " << k << std::endl;
376 
377  calculateJacobianX(currentValue);
378 
379  if (solveJacobianXeqB(mH, mdxdt) != 0.0)
380  {
381  // We need to check that mH != 0
382  C_FLOAT64 * pH = mH.array();
383  C_FLOAT64 * pHEnd = pH + mH.size();
384 
385  for (; pH != pHEnd; ++pH)
386  {
387  if (fabs(*pH) > 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon())
388  break;
389  }
390 
391  if (pH == pHEnd)
392  {
393  if (mKeepProtocol)
394  mMethodLog << " Newton step failed. Jacobian could not be inverted.\n\n";
395 
397  }
398  }
399 
400  C_FLOAT64 newValue = currentValue * 1.001;
401 
402  //repeat till the new max rate is smaller than the old.
403  //max 32 times
404  size_t i;
405 
406  for (i = 0; (i < 32) && !((newValue < currentValue)); i++)
407  {
408  C_FLOAT64 * pXit = mpX;
409  C_FLOAT64 * pXoldIt = mXold.array();
410  C_FLOAT64 * pHit = mH.array();
411  C_FLOAT64 * pEnd = pHit + mDimension;
412 
413  for (; pHit != pEnd; ++pHit, ++pXit, ++pXoldIt)
414  {
415  *pXit = *pXoldIt - *pHit;
416  (*pHit) *= 0.5;
417  }
418 
420  newValue = targetFunction(mdxdt);
421 
422  // mpParentTask->output(COutputInterface::DURING);
423  }
424 
425  if (i == 32)
426  {
427  //discard the step
428  memcpy(mpX, mXold.array(), mDimension * sizeof(C_FLOAT64));
429 
431  currentValue = targetFunction(mdxdt);
432 
433  if (mKeepProtocol) mMethodLog << " Newton step failed. Damping limit exceeded.\n";
434 
436 
437  // if (isSteadyState(oldMaxRate) && (mAcceptNegative || allPositive()))
438  // ReturnCode = CNewtonMethod::found;
439  // else if (oldMaxRate < *mpSSResolution)
440  // ReturnCode = CNewtonMethod::notFound;
441  // else
442  // ReturnCode = CNewtonMethod::dampingLimitExceeded;
443 
444  //if (mpCallBack) mpCallBack->finish(hProcess);
445  }
446 
447  if (!(mAcceptNegative || allPositive()))
448  {
449  if (mKeepProtocol)
450  mMethodLog << " Newton step failed. Negative volume or concentration found.\n\n";
451 
453  }
454 
455  currentValue = newValue; //return the new target value
456 
457  if (mKeepProtocol)
458  {
459  if (i <= 1)
460  mMethodLog << " Regular Newton step. New value: " << currentValue << "\n";
461  else
462  mMethodLog << " Newton step with damping. New value: " << currentValue
463  << " (" << i - 1 << " damping iteration(s))\n";
464  }
465 
467 }
468 
469 //************************************************************
470 
472 {
474  unsigned C_INT32 k;
475 
476  k = 0;
477  //start progress bar
478  size_t hProcess;
479 
480  if (mpCallBack)
481  hProcess = mpCallBack->addItem("Newton method...",
482  k,
483  & mIterationLimit);
484 
485  C_FLOAT64 targetValue;
486 
488  targetValue = targetFunction(mdxdt);
489 
490  {
491  if (mKeepProtocol) mMethodLog << " Starting Newton Iterations...\n";
492 
493  for (k = 0; k < mIterationLimit && !isSteadyState(targetValue); k++)
494  {
495  if (mpCallBack && !mpCallBack->progressItem(hProcess)) break;
496 
497  result = doNewtonStep(targetValue);
498 
499  if (singularJacobian == result) break;
500 
501  if (dampingLimitExceeded == result) break;
502 
503  if (negativeValueFound == result) break;
504  }
505  }
506 
507  //check if ss was found. If not make sure the correct return value is set
508  if (isSteadyState(targetValue))
509  result = CNewtonMethod::found;
510  else if (CNewtonMethod::stepSuccesful == result)
512 
513  //log
514  if (mKeepProtocol)
515  {
516  if (CNewtonMethod::found == result)
517  mMethodLog << " Success: Target criterium matched by " << targetValue << ".\n";
518  else if (CNewtonMethod::dampingLimitExceeded == result)
519  mMethodLog << " Failed: Target criterium not matched after reaching iteration limit. " << targetValue << "\n";
520  }
521 
522  //do an additional Newton step to refine the result
523  if ((CNewtonMethod::found == result) && mForceNewton && targetValue > 0.0)
524  {
525  bool tmp = true;
526 
527  ++k; if (mpCallBack && !mpCallBack->progressItem(hProcess)) tmp = false;
528 
529  if (tmp)
530  {
531  if (mKeepProtocol) mMethodLog << " Do additional step to refine result...\n";
532 
533  result = doNewtonStep(targetValue);
534 
535  if (CNewtonMethod::stepSuccesful == result)
536  result = CNewtonMethod::found;
537 
538  if (CNewtonMethod::singularJacobian == result)
539  {
540  if (mKeepProtocol) mMethodLog << " Additional step failed. Old values restored.\n";
541 
542  result = CNewtonMethod::found;
543  }
544 
546  {
547  if (mKeepProtocol) mMethodLog << " Additional step failed. Old values restored.\n";
548 
549  result = CNewtonMethod::found;
550  }
551 
552  if (CNewtonMethod::negativeValueFound == result)
553  {
554  if (mKeepProtocol) mMethodLog << " Additional step failed. Old values restored.\n";
555 
556  result = CNewtonMethod::found;
557  }
558  }
559  }
560 
561  //end progress bar
562  if (mpCallBack) mpCallBack->finishItem(hProcess);
563 
564  return result;
565 }
566 
568 {
572 }
573 
575 {
576  //checks for NaNs
577  const C_FLOAT64 * pIt = mpSteadyState->beginIndependent();
578  const C_FLOAT64 * pEnd = mpSteadyState->endIndependent();
579 
580  for (; pIt != pEnd; ++pIt)
581  if (isnan(*pIt))
582  return true;
583 
584  return false;
585 }
586 
588 {
589  if (value > *mpSSResolution)
590  return false;
591 
592  if (containsNaN())
593  return false;
594 
595  return true;
596 }
597 
599 {
600  // New criterion: We solve Jacobian * x = current rates and compare x with the current state
601  // Calculate the Jacobian
603 
604  CVector< C_FLOAT64 > Distance;
605 
606  C_FLOAT64 Error = solveJacobianXeqB(Distance, particlefluxes);
607 
608  // We look at all ODE determined entity and dependent species rates.
609  C_FLOAT64 * pDistance = Distance.array();
610  C_FLOAT64 * pDistanceEnd = pDistance + Distance.size();
611  C_FLOAT64 * pCurrentState = mpSteadyState->beginIndependent();
612  const C_FLOAT64 * pAtol = mAtol.array();
613 
614  // Assure that all values are updated.
616 
617  CModelEntity *const* ppEntity = mpModel->getStateTemplate().beginIndependent();
618  const CMetab * pMetab = NULL;
619  C_FLOAT64 Number2Quantity = mpModel->getNumber2QuantityFactor();
620 
621  C_FLOAT64 AbsoluteDistance = 0.0; // Largest relative distance
622  C_FLOAT64 RelativeDistance = 0.0; // Total relative distance
623 
624  C_FLOAT64 tmp;
625 
626  for (; pDistance != pDistanceEnd; ++pDistance, ++pCurrentState, ++pAtol, ++ppEntity)
627  {
628  // Prevent division by 0
629  tmp = fabs(*pDistance) / std::max(fabs(*pCurrentState), *pAtol);
630  RelativeDistance += tmp * tmp;
631 
632  tmp = fabs(*pDistance);
633 
634  if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
635  {
636  tmp *= Number2Quantity / fabs(pMetab->getCompartment()->getValue());
637  }
638 
639  AbsoluteDistance += tmp * tmp;
640  }
641 
642  RelativeDistance =
643  isnan(RelativeDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(RelativeDistance);
644  AbsoluteDistance =
645  isnan(AbsoluteDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(AbsoluteDistance);
646 
647  C_FLOAT64 TargetValue = std::max(RelativeDistance, AbsoluteDistance);
648 
649  if (Error < TargetValue)
650  return TargetValue * (1.0 + Error);
651  else
652  return Error;
653 }
654 
655 //virtual
657 {
658  if (!CSteadyStateMethod::isValidProblem(pProblem)) return false;
659 
660  const CModel * pModel = pProblem->getModel();
661 
662  if (!pModel->isAutonomous() &&
663  *getValue("Use Newton").pBOOL)
665 
666  //const CSteadyStateProblem * pP = dynamic_cast<const CSteadyStateProblem *>(pProblem);
667 
668  if (!((* getValue("Use Newton").pBOOL)
669  || (* getValue("Use Integration").pBOOL)
670  || (* getValue("Use Back Integration").pBOOL)))
671  {
672  //would do nothing
673  CCopasiMessage(CCopasiMessage::ERROR, "At least one of the features \n - UseNewton\n - UseIntegration\n - UseBackIntegration\nmust be activated.");
674  return false;
675  }
676 
677  if (*getValue("Maximum duration for forward integration").pUDOUBLE <= 0)
678  {
679  CCopasiMessage(CCopasiMessage::ERROR, "Maximum duration for forward integration needs to be positive.");
680  return false;
681  }
682 
683  if (*getValue("Maximum duration for backward integration").pUDOUBLE <= 0)
684  {
685  CCopasiMessage(CCopasiMessage::ERROR, "Maximum duration for backward integration needs to be positive.");
686  return false;
687  }
688 
689  return true;
690 }
691 
693 {
694  if (!CSteadyStateMethod::initialize(pProblem)) return false;
695 
696  CTrajectoryProblem * pTrajectoryProblem = NULL;
697  CTrajectoryMethod * pTrajectoryMethod = NULL;
698 
699  cleanup();
700 
701  /* Configure Newton */
703  = mForceNewton = mKeepProtocol = false;
704 
705  if (* getValue("Use Newton").pBOOL)
706  mUseNewton = true;
707 
708  if (* getValue("Use Integration").pBOOL)
709  mUseIntegration = true;
710 
711  if (* getValue("Use Back Integration").pBOOL)
712  mUseBackIntegration = true;
713 
714  if (* getValue("Accept Negative Concentrations").pBOOL)
715  mAcceptNegative = true;
716 
717  //if (* getValue("Force additional Newton step").pBOOL)
718  mForceNewton = true;
719  //if (* getValue("Keep Protocol").pBOOL)
720  mKeepProtocol = true;
721 
722  mIterationLimit = * getValue("Iteration Limit").pUINT;
723 
724  mMaxDurationForward = *getValue("Maximum duration for forward integration").pUDOUBLE;
725  mMaxDurationBackward = *getValue("Maximum duration for backward integration").pUDOUBLE;
726 
727  //mFactor = * getValue("Derivation Factor").pUDOUBLE;
728  //mSSResolution = * getValue("Steady State Resolution").pUDOUBLE;
729  //mScaledResolution =
730  // mSSResolution; // * initialState.getModel()->getQuantity2NumberFactor();
731  // :TODO: discuss scaling
732 
733  // convert CState to CStateX
734  //mInitialStateX = mpProblem->getInitialState();
735  //*mpSteadyStateX = mInitialStateX; //not strictly necessary
736 
738 
743  mIpiv = new C_INT [mDimension];
744 
746  {
747  // create an appropriate trajectory task
748  CCopasiDataModel* pDataModel = getObjectDataModel();
749  assert(pDataModel != NULL);
750  CTrajectoryTask * pSrc =
751  dynamic_cast< CTrajectoryTask * >((*pDataModel->getTaskList())["Time-Course"]);
752 
753  if (pSrc)
754  mpTrajectory = new CTrajectoryTask(*pSrc, this);
755  else
756  mpTrajectory = new CTrajectoryTask(this);
757 
758  pTrajectoryProblem =
759  dynamic_cast<CTrajectoryProblem *>(mpTrajectory->getProblem());
760  assert(pTrajectoryProblem);
761 
764 
765  pTrajectoryMethod =
766  dynamic_cast<CTrajectoryMethod *>(mpTrajectory->getMethod());
767  assert(pTrajectoryMethod);
768 
769  pTrajectoryProblem->setModel(mpProblem->getModel());
770  pTrajectoryProblem->setStepNumber(1);
771 
773  }
774 
775  return true;
776 }
777 
779 {
780  X = B;
781 
782  C_INT M = (C_INT) mpJacobianX->numCols();
783  C_INT N = (C_INT) mpJacobianX->numRows();
784 
785  if (M == 0 || N == 0 || M != N)
786  {
787  return std::numeric_limits< C_FLOAT64 >::infinity();
788  }
789 
790  C_INT LDA = std::max< C_INT >(1, M);
791  C_INT NRHS = 1;
792 
793  // We need the transpose of the Jacobian;
794  CMatrix< C_FLOAT64 > JT(M, N);
795  C_FLOAT64 * mpJ = mpJacobianX->array();
796  C_FLOAT64 * mpJTcolumn = JT.array();
797  C_FLOAT64 * mpJTcolumnEnd = mpJTcolumn + M;
798  C_FLOAT64 * mpJT = JT.array();
799  C_FLOAT64 * mpJTEnd = mpJT + JT.size();
800 
801  for (; mpJTcolumn != mpJTcolumnEnd; ++mpJTcolumn)
802  {
803  mpJT = mpJTcolumn;
804 
805  for (; mpJT < mpJTEnd; mpJT += M, ++mpJ)
806  {
807  *mpJT = *mpJ;
808  }
809  }
810 
811  CVector< C_INT > JPVT(M);
812  JPVT = 0;
813 
814  C_FLOAT64 RCOND = 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon();
815 
816  C_INT RANK = 0;
817 
818  CVector< C_FLOAT64 > WORK(1);
819  C_INT LWORK = -1;
820  C_INT INFO;
821 
822  /*
823  SUBROUTINE DGELSY(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
824  $ WORK, LWORK, INFO)
825  *
826  * -- LAPACK driver routine (version 3.2) --
827  * -- LAPACK is a software package provided by Univ. of Tennessee, --
828  * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
829  * November 2006
830  *
831  * .. Scalar Arguments ..
832  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
833  DOUBLE PRECISION RCOND
834  * ..
835  * .. Array Arguments ..
836  INTEGER JPVT(*)
837  DOUBLE PRECISION A(LDA, * ), B(LDB, * ), WORK(*)
838  * ..
839  *
840  * Purpose
841  * =======
842  *
843  * DGELSY computes the minimum-norm solution to a real linear least
844  * squares problem:
845  * minimize || A * X - B ||
846  * using a complete orthogonal factorization of A. A is an M-by-N
847  * matrix which may be rank-deficient.
848  *
849  * Several right hand side vectors b and solution vectors x can be
850  * handled in a single call; they are stored as the columns of the
851  * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
852  * matrix X.
853  *
854  * The routine first computes a QR factorization with column pivoting:
855  * A * P = Q * [ R11 R12 ]
856  * [ 0 R22 ]
857  * with R11 defined as the largest leading submatrix whose estimated
858  * condition number is less than 1/RCOND. The order of R11, RANK,
859  * is the effective rank of A.
860  *
861  * Then, R22 is considered to be negligible, and R12 is annihilated
862  * by orthogonal transformations from the right, arriving at the
863  * complete orthogonal factorization:
864  * A * P = Q * [ T11 0 ] * Z
865  * [ 0 0 ]
866  * The minimum-norm solution is then
867  * X = P * Z' [ inv(T11)*Q1'*B ]
868  * [ 0 ]
869  * where Q1 consists of the first RANK columns of Q.
870  *
871  * This routine is basically identical to the original xGELSX except
872  * three differences:
873  * o The call to the subroutine xGEQPF has been substituted by the
874  * the call to the subroutine xGEQP3. This subroutine is a Blas-3
875  * version of the QR factorization with column pivoting.
876  * o Matrix B (the right hand side) is updated with Blas-3.
877  * o The permutation of matrix B (the right hand side) is faster and
878  * more simple.
879  *
880  * Arguments
881  * =========
882  *
883  * M (input) INTEGER
884  * The number of rows of the matrix A. M >= 0.
885  *
886  * N (input) INTEGER
887  * The number of columns of the matrix A. N >= 0.
888  *
889  * NRHS (input) INTEGER
890  * The number of right hand sides, i.e., the number of
891  * columns of matrices B and X. NRHS >= 0.
892  *
893  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
894  * On entry, the M-by-N matrix A.
895  * On exit, A has been overwritten by details of its
896  * complete orthogonal factorization.
897  *
898  * LDA (input) INTEGER
899  * The leading dimension of the array A. LDA >= max(1,M).
900  *
901  * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
902  * On entry, the M-by-NRHS right hand side matrix B.
903  * On exit, the N-by-NRHS solution matrix X.
904  *
905  * LDB (input) INTEGER
906  * The leading dimension of the array B. LDB >= max(1,M,N).
907  *
908  * JPVT (input/output) INTEGER array, dimension (N)
909  * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
910  * to the front of AP, otherwise column i is a free column.
911  * On exit, if JPVT(i) = k, then the i-th column of AP
912  * was the k-th column of A.
913  *
914  * RCOND (input) DOUBLE PRECISION
915  * RCOND is used to determine the effective rank of A, which
916  * is defined as the order of the largest leading triangular
917  * submatrix R11 in the QR factorization with pivoting of A,
918  * whose estimated condition number < 1/RCOND.
919  *
920  * RANK (output) INTEGER
921  * The effective rank of A, i.e., the order of the submatrix
922  * R11. This is the same as the order of the submatrix T11
923  * in the complete orthogonal factorization of A.
924  *
925  * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
926  * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
927  *
928  * LWORK (input) INTEGER
929  * The dimension of the array WORK.
930  * The unblocked strategy requires that:
931  * LWORK >= MAX(MN+3*N+1, 2*MN+NRHS ),
932  * where MN = min(M, N ).
933  * The block algorithm requires that:
934  * LWORK >= MAX(MN+2*N+N*(N+1), 2*MN+NB*NRHS ),
935  * where NB is an upper bound on the blocksize returned
936  * by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
937  * and DORMRZ.
938  *
939  * If LWORK = -1, then a workspace query is assumed; the routine
940  * only calculates the optimal size of the WORK array, returns
941  * this value as the first entry of the WORK array, and no error
942  * message related to LWORK is issued by XERBLA.
943  *
944  * INFO (output) INTEGER
945  * = 0: successful exit
946  * < 0: If INFO = -i, the i-th argument had an illegal value.
947  *
948  * Further Details
949  * ===============
950  *
951  * Based on contributions by
952  * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
953  * E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
954  * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
955  *
956  * =====================================================================
957  */
958 
959  dgelsy_(&M, &N, &NRHS, JT.array(), &LDA, X.array(), &LDA, JPVT.array(), &RCOND, &RANK,
960  WORK.array(), &LWORK, &INFO);
961 
962  if (INFO < 0)
963  {
964  return std::numeric_limits< C_FLOAT64 >::infinity();
965  }
966 
967  LWORK = (C_INT) WORK[0];
968  WORK.resize(LWORK);
969 
970  dgelsy_(&M, &N, &NRHS, JT.array(), &LDA, X.array(), &LDA, JPVT.array(), &RCOND, &RANK,
971  WORK.array(), &LWORK, &INFO);
972 
973  if (INFO < 0)
974  {
975  return std::numeric_limits< C_FLOAT64 >::infinity();
976  }
977 
978  C_FLOAT64 Error = 0;
979 
980  if (RANK != M)
981  {
982  // We need to check whether the || Ax - b || is sufficiently small.
983  // Calculate Ax
984  char T = 'N';
985  M = 1;
986  C_FLOAT64 Alpha = 1.0;
987  C_FLOAT64 Beta = 0.0;
988 
989  CVector< C_FLOAT64 > Ax = B;
990 
991  dgemm_(&T, &T, &M, &N, &N, &Alpha, X.array(), &M,
992  mpJacobianX->array(), &N, &Beta, Ax.array(), &M);
993 
994  // Calculate absolute and relative error
995  C_FLOAT64 *pAx = Ax.array();
996  C_FLOAT64 *pAxEnd = pAx + Ax.size();
997  const C_FLOAT64 *pB = B.array();
998  C_FLOAT64 * pCurrentState = mpSteadyState->beginIndependent();
999  const C_FLOAT64 * pAtol = mAtol.array();
1000 
1001  // Assure that all values are updated.
1003 
1004  CModelEntity *const* ppEntity = mpModel->getStateTemplate().beginIndependent();
1005  const CMetab * pMetab = NULL;
1006  C_FLOAT64 Number2Quantity = mpModel->getNumber2QuantityFactor();
1007 
1008  C_FLOAT64 AbsoluteDistance = 0.0; // Largest relative distance
1009  C_FLOAT64 RelativeDistance = 0.0; // Total relative distance
1010 
1011  C_FLOAT64 tmp;
1012 
1013  for (; pAx != pAxEnd; ++pAx, ++pB, ++pCurrentState, ++pAtol, ++ppEntity)
1014  {
1015  // Prevent division by 0
1016  tmp = fabs(*pAx - *pB) / std::max(fabs(*pCurrentState), *pAtol);
1017  RelativeDistance += tmp * tmp;
1018 
1019  tmp = fabs(*pAx - *pB);
1020 
1021  if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
1022  {
1023  tmp *= Number2Quantity / fabs(pMetab->getCompartment()->getValue());
1024  }
1025 
1026  AbsoluteDistance += tmp * tmp;
1027  }
1028 
1029  RelativeDistance =
1030  isnan(RelativeDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(RelativeDistance);
1031  AbsoluteDistance =
1032  isnan(AbsoluteDistance) ? std::numeric_limits< C_FLOAT64 >::infinity() : sqrt(AbsoluteDistance);
1033 
1034  Error = std::max(RelativeDistance, AbsoluteDistance);
1035  }
1036 
1037  return Error;
1038 }
CCopasiDataModel * getObjectDataModel()
CVector< C_FLOAT64 > mH
Definition: CNewtonMethod.h:65
#define C_INT
Definition: copasi.h:115
virtual bool setCallBack(CProcessReport *pCallBack)
#define pdelete(p)
Definition: copasi.h:215
virtual bool initialize(const OutputFlag &of, COutputHandler *pOutputHandler, std::ostream *pOstream)
CMatrix< C_FLOAT64 > * mpJacobianX
CCopasiProblem * getProblem()
C_FLOAT64 mMaxDurationBackward
Definition: CNewtonMethod.h:60
int dgelsy_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *jpvt, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, integer *info)
size_t getNumIndependent() const
Definition: CState.cpp:222
virtual bool isValidProblem(const CCopasiProblem *pProblem)
virtual bool elevateChildren()
virtual bool setName(const std::string &name)
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
virtual size_t numRows() const
Definition: CMatrix.h:138
#define fatalError()
void initializeParameter()
virtual bool process(const bool &useInitialValues)
std::ostringstream mMethodLog
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CState mStartState
Definition: CNewtonMethod.h:72
void setDuration(const C_FLOAT64 &duration)
std::string getVersion()
Definition: CReadConfig.cpp:76
CVector< C_FLOAT64 > mdxdt
Definition: CNewtonMethod.h:67
const CCopasiMethod::SubType & getSubType() const
C_FLOAT64 targetFunction(const CVector< C_FLOAT64 > &particleFluxes)
#define C_UNUSED(p)
Definition: copasi.h:220
#define C_INT32
Definition: copasi.h:90
Definition: CMetab.h:178
bool mUseBackIntegration
Definition: CNewtonMethod.h:54
bool isSteadyState(C_FLOAT64 value)
virtual bool progressItem(const size_t &handle)
virtual bool setModel(CModel *pModel)
C_FLOAT64 * endIndependent()
Definition: CState.cpp:329
bool removeParameter(const std::string &name)
CVector< C_FLOAT64 > mXold
Definition: CNewtonMethod.h:66
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
void setState(const CState &state)
Definition: CModel.cpp:1785
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
void calculateJacobianX(const C_FLOAT64 &oldMaxRate)
virtual CSteadyStateMethod::ReturnCode processInternal()
CNewtonMethod::NewtonResultCode processNewton()
const bool & isAutonomous() const
Definition: CModel.cpp:3956
const Value & getValue() const
CCopasiVectorN< CCopasiTask > * getTaskList()
void setStepNumber(const unsigned C_INT32 &stepNumber)
bool setValue(const std::string &name, const CType &value)
virtual CSteadyStateMethod::ReturnCode returnProcess(bool steadyStateFound)
unsigned C_INT32 * pUINT
CVector< C_FLOAT64 > mAtol
Definition: CNewtonMethod.h:64
virtual bool finishItem(const size_t &handle)
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
virtual bool isValidProblem(const CCopasiProblem *pProblem)
CCopasiParameter * getParameter(const std::string &name)
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
CTrajectoryTask * mpTrajectory
Definition: CNewtonMethod.h:70
virtual bool initialize(const CSteadyStateProblem *pProblem)
const CSteadyStateProblem * mpProblem
C_FLOAT64 * mpX
Definition: CNewtonMethod.h:63
size_t size() const
Definition: CVector.h:100
unsigned C_INT32 mIterationLimit
Definition: CNewtonMethod.h:58
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
const C_FLOAT64 & getValue() const
virtual void load(CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)
CNewtonMethod(const CCopasiContainer *pParent=NULL)
CCopasiMethod * getMethod()
void calculateDerivativesX()
#define C_FLOAT64
Definition: copasi.h:92
CNewtonMethod::NewtonResultCode doNewtonStep(C_FLOAT64 &currentValue)
virtual bool initialize(const CSteadyStateProblem *pProblem)
CType * array()
Definition: CVector.h:139
C_FLOAT64 solveJacobianXeqB(CVector< C_FLOAT64 > &X, const CVector< C_FLOAT64 > &B) const
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
#define MCSteadyState
virtual size_t size() const
Definition: CMatrix.h:132
C_FLOAT64 * mpSSResolution
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
size_t mDimension
Definition: CNewtonMethod.h:62
CVector< C_FLOAT64 > initializeAtolVector(const C_FLOAT64 &baseTolerance, const bool &reducedModel) const
Definition: CModel.cpp:4368
Definition: CModel.h:50
C_INT32 getVariable(const std::string &name, const std::string &type, void *pout, CReadConfig::Mode mode=CReadConfig::NEXT)
Definition: CReadConfig.cpp:81
bool mUseIntegration
Definition: CNewtonMethod.h:53
CNewtonMethod::NewtonResultCode doIntegration(bool forward)
virtual size_t numCols() const
Definition: CMatrix.h:144
bool containsNaN() const
virtual bool setMethodType(const int &type)
const CCompartment * getCompartment() const
Definition: CMetab.cpp:222
CProcessReport * mpCallBack
virtual CType * array()
Definition: CMatrix.h:337
CModel * getModel() const
C_FLOAT64 mMaxDurationForward
Definition: CNewtonMethod.h:59
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
bool mAcceptNegative
Definition: CNewtonMethod.h:55
#define max(a, b)
Definition: f2c.h:176