Merge branch 'develop' into stable

fabian.froehlich · fabian.froehlich · commit 8e3f5d50faea · 2016-09-02T15:22:57.000+02:00
# Conflicts:
#	@amimodel/makeEvents.m
diff --git a/@amimodel/amimodel.m b/@amimodel/amimodel.m
@@ -92,7 +92,7 @@
 
         % counter that allows enforcing of recompilation of models after
         % code changes
-        compver = 15;
+        compver = 16;
     end
     
     properties ( GetAccess = 'public', SetAccess = 'public' )
diff --git a/@amimodel/generateC.m b/@amimodel/generateC.m
@@ -458,6 +458,7 @@ function generateC(this)
             % they should act as placeholders
             fprintf(fid,'                    return(0);\n');
         else
+            fprintf(fid,['                    mexWarnMsgIdAndTxt("AMICI:mex:' iffun{1} ':NotAvailable","ERROR: The function ' iffun{1} ' was called but not compiled for this model.");\n']);
             fprintf(fid,'                    return(-1);\n');
         end
     end
diff --git a/@amimodel/generateM.m b/@amimodel/generateM.m
@@ -377,13 +377,15 @@ function generateM(this, amimodelo2)
         fprintf(fid,'    sol.sx = bsxfun(@times,sol.sx,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
         fprintf(fid,'    sol.sy = bsxfun(@times,sol.sy,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
         fprintf(fid,'    sol.sz = bsxfun(@times,sol.sz,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
+        fprintf(fid,'    sol.srz = bsxfun(@times,sol.srz,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
         fprintf(fid,'    sol.ssigmay = bsxfun(@times,sol.ssigmay,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
         fprintf(fid,'    sol.ssigmaz = bsxfun(@times,sol.ssigmaz,permute(theta(options_ami.sens_ind),[3,2,1]));\n');
     case 'log10'
         fprintf(fid,'    sol.sllh = sol.sllh.*theta(options_ami.sens_ind)*log(10);\n');
         fprintf(fid,'    sol.sx = bsxfun(@times,sol.sx,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
         fprintf(fid,'    sol.sy = bsxfun(@times,sol.sy,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
         fprintf(fid,'    sol.sz = bsxfun(@times,sol.sz,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
+        fprintf(fid,'    sol.srz = bsxfun(@times,sol.srz,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
         fprintf(fid,'    sol.ssigmay = bsxfun(@times,sol.ssigmay,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
         fprintf(fid,'    sol.ssigmayz = bsxfun(@times,sol.ssigmaz,permute(theta(options_ami.sens_ind),[3,2,1])*log(10));\n');
     otherwise
@@ -581,6 +583,7 @@ function generateM(this, amimodelo2)
 
 for fun = this.mfuns
     if(isfield(this.fun,fun{1}))
+        fprintf([fun{1} ' | ']);
         this.fun.(fun{1}).writeMcode(this);
     end
 end
diff --git a/@amimodel/makeEvents.m b/@amimodel/makeEvents.m
diff --git a/@amimodel/parseModel.m b/@amimodel/parseModel.m
@@ -82,13 +82,13 @@ function parseModel(this)
     % compute functions
     
     % do not change the ordering, it is essential for correct dependencies
-    funs = {'xdot','w','dwdx','dwdp','J','x0','Jv','JBand','JSparse','y','z','deltax','dydp','dxdotdp','root','Jy','dJydx','dJydp','sJy','Jz','dJzdx','dJzdp','sJz'};
+    funs = {'xdot','w','dwdx','dwdp','J','x0','Jv','JBand','JSparse','y','z','deltax','dydp','dxdotdp','root','Jy','Jz','sigma_y','sigma_z'};
     
     if(this.forward)
-        funs = {funs{:},'sxdot','sx0','sy','sz','sz_tf','deltasx','stau','sroot'};
+        funs = {funs{:},'sxdot','sx0','sy','sz','sz_tf','deltasx','stau','sroot','sJy','sJz'};
     end
     if(this.adjoint)
-        funs = {funs{:},'xBdot','qBdot','JB','JvB','JBandB','JSparseB','dydx','dzdx','dzdp','deltaxB','deltaqB','sigma_y','sigma_z','dsigma_ydp','dsigma_zdp','sx0'};
+        funs = {funs{:},'xBdot','qBdot','JB','JvB','JBandB','JSparseB','dydx','dzdx','dzdp','deltaxB','deltaqB','dsigma_ydp','dsigma_zdp','sx0','dJydx','dJydp','dJzdx','dJzdp'};
     end
     if(this.o2flag)
         funs = {funs{:},'s2root'};
@@ -102,7 +102,8 @@ function parseModel(this)
     
     this.funs = funs;
     
-    this.mfuns = {'J','dxdotdp'};
+    %this.mfuns = {'J','dxdotdp'};
+    this.mfuns = {};
     
     % compute symbolic expressions
     for ifun = 1:length(funs)
diff --git a/amiwrap.m b/amiwrap.m
@@ -96,6 +96,7 @@ function amiwrap( varargin )
     end
     
     % generate the matlab wrapper
+    disp('Generating M code ...')
     if(o2flag)
         model.generateM(modelo2);
     else
diff --git a/auxiliary/ami_mfun.m b/auxiliary/ami_mfun.m
@@ -182,6 +182,11 @@ function writeOutput(fid,expr)
         t = regexp(str,[var{1} '_([0-9]*)'],'tokens');
         ts = cellfun(@(x) [var{1} '_' x{1} ],t,'UniformOutput',false);
         td = cellfun(@(x) [var{1} '(' num2str(str2double(x)+1) ')'],t,'UniformOutput',false);
+        ts = unique(ts);
+        td = unique(td);
+        [~,idx_ts] = sort(cellfun(@(x) length(x),ts),'descend');
+        ts = ts(idx_ts);
+        td = ts(idx_ts);
         for it = 1:length(ts)
             str = strrep(str,ts{it},td{it});
         end
diff --git a/include/symbolic_functions.h b/include/symbolic_functions.h
@@ -7,10 +7,10 @@
 double amilog(double x);
 double heaviside(double x);
 
-double am_min(double a, double b);
-double Dam_min(int id,double a, double b);
-double am_max(double a, double b);
-double Dam_max(int id,double a, double b);
+double am_min(double a, double b, double c);
+double Dam_min(int id,double a, double b, double c);
+double am_max(double a, double b, double c);
+double Dam_max(int id,double a, double b, double c);
 
 /* sign */
 double sign(double x);
diff --git a/src/symbolic_functions.c b/src/symbolic_functions.c
@@ -102,24 +102,26 @@ double sign(double x) {
  *
  * @param a value1 @type double
  * @param b value2 @type double
+ * @param c bogus parameter do ensure correct parsing as a function @type double
  * @return if(a < b) then a else b @type double
  *
  */
-double am_min(double a, double b) {
+double am_min(double a, double b, double c) {
     return(fmin(a,b));
 }
 
 /**
  * parameter derivative of c implementation of matlab function min
  *
  * @param id argument index for differentiation
- * @param a bool1 @type double
- * @param b bool2 @type double
+ * @param a value1 @type double
+ * @param b value2 @type double
+ * @param c bogus parameter do ensure correct parsing as a function @type double
  * @return id == 1:  if(a < b) then 1 else 0 @type double
  * @return id == 2:  if(a < b) then 0 else 1 @type double
  *
  */
-double Dam_min(int id,double a, double b) {
+double Dam_min(int id,double a, double b, double c) {
     if (id == 1) {
         if (a < b) {
             return(1);
@@ -140,24 +142,26 @@ double Dam_min(int id,double a, double b) {
  *
  * @param a value1 @type double
  * @param b value2 @type double
+ * @param c bogus parameter do ensure correct parsing as a function @type double
  * @return if(a > b) then a else b @type double
  *
  */
-double am_max(double a, double b) {
+double am_max(double a, double b, double c) {
     return(fmax(a,b));
 }
 
 /**
  * parameter derivative of c implementation of matlab function max
  *
  * @param id argument index for differentiation
- * @param a bool1 @type double
- * @param b bool2 @type double
+ * @param a value1 @type double
+ * @param b value2 @type double
+ * @param c bogus parameter do ensure correct parsing as a function @type double
  * @return id == 1:  if(a > b) then 1 else 0 @type double
  * @return id == 2:  if(a > b) then 0 else 1 @type double
  *
  */
-double Dam_max(int id,double a, double b) {
+double Dam_max(int id,double a, double b, double c) {
     if (id == 1) {
         if (a > b) {
             return(1);
diff --git a/sundials-2.6.2/src/cvodes/cvodes.c b/sundials-2.6.2/src/cvodes/cvodes.c
@@ -2962,7 +2962,7 @@ int CVode(void *cvode_mem, realtype tout, N_Vector yout,
     next_q = q;
     
     /* Reset and check ewt, ewtQ, ewtS */   
-    if (nst > 0) {
+    if (!reinit) {
 
       ier = efun(zn[0], ewt, e_data);
       if(ier != 0) {
@@ -5113,7 +5113,7 @@ static int cvStep(CVodeMem cv_mem)
 
   /* If needed, adjust method parameters */
 
-  if ((nst > 0) && (hprime != h)) cvAdjustParams(cv_mem);
+  if ((!reinit) && (hprime != h)) cvAdjustParams(cv_mem);
 
   /* Looping point for attempts to take a step */
 
@@ -5670,7 +5670,7 @@ static void cvSet(CVodeMem cv_mem)
   rl1 = ONE / l[1];
   gamma = h * rl1;
   if (reinit) gammap = gamma;
-  gamrat = (nst > 0) ? gamma / gammap : ONE;  /* protect x / x != 1.0 */
+  gamrat = (!reinit) ? gamma / gammap : ONE;  /* protect x / x != 1.0 */
 }
 
 /*
diff --git a/symbolic/am_ceil.m b/symbolic/am_ceil.m
diff --git a/symbolic/am_if.m b/symbolic/am_if.m
@@ -2,5 +2,17 @@
 % syms x y z
 % f = symfun(sym('am_if(x,y,z)'),[x y z]);
 % fun = f(condition, truepart, falsepart);
-fun = falsepart + heaviside(condition)*(truepart-falsepart); 
+if(islogical(condition))
+    if(condition)
+        fun = truepart;
+    else
+        fun = falsepart;
+    end
+else
+    if(logical(condition~=0))
+        fun = falsepart + heaviside(condition)*(truepart-falsepart);
+    else
+        fun = falsepart;
+    end
+end
 end
diff --git a/symbolic/am_max.m b/symbolic/am_max.m
@@ -2,6 +2,6 @@
 % syms x y
 % f = symfun(sym('am_max(x,y)'),[x y]);
 % fun = f(a,b);
-fun = sym(['am_max(' char(vpa(a)) ',' char(vpa(b)) ')']);
+fun = sym(['am_max(' char(sym(a)) ',' char(sym(b)) ',0.0)']);
 % fun = b + heaviside(a-b)*(a-b);
 end
diff --git a/symbolic/am_min.m b/symbolic/am_min.m
@@ -2,6 +2,6 @@
 % syms x y
 % f = symfun(sym('am_min(x,y)'),[x y]);
 % fun = f(a,b);
-fun = sym(['am_min(' char(vpa(a)) ',' char(vpa(b)) ')']);
+fun = sym(['am_min(' char(sym(a)) ',' char(sym(b)) ',0.0)']);
 % fun = b + heaviside(b-a)*(a-b);
 end
