From dd09264b03ee9a7534c8532d2fa119da81a710f4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 2 Apr 2019 16:39:32 +0200
Subject: [PATCH] k-order DLL: return unfolded matrices in 5th output argument

Thus we can remove unfolding code from k_order_pert.m.
---
 matlab/k_order_pert.m                         | 81 ++-----------------
 .../k_order_perturbation.cc                   | 11 ++-
 2 files changed, 15 insertions(+), 77 deletions(-)

diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index 596dc100f..892cec68b 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -59,15 +59,13 @@ switch(order)
     dr.g_2 = g_2;
   case 3
     if options.pruning
-        [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr, ...
-                                                          M,options);
+        [err, g_0, g_1, g_2, g_3, derivs] = k_order_perturbation(dr,M,options);
         if err
             info(1)=9;
             return
         end
     else
-        [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr, ...
-                                                         M,options);
+        [err, g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
         if err
             info(1)=9;
             return
@@ -86,14 +84,14 @@ end
 if options.pruning && order == 3
     dr.ghx = derivs.gy;
     dr.ghu = derivs.gu;
-    dr.ghxx = unfold2(derivs.gyy,nspred);
+    dr.ghxx = derivs.gyy;
     dr.ghxu = derivs.gyu;
-    dr.ghuu = unfold2(derivs.guu,exo_nbr);
+    dr.ghuu = derivs.guu;
     dr.ghs2 = derivs.gss;
-    dr.ghxxx = unfold3(derivs.gyyy,nspred);
-    dr.ghxxu = unfold21(derivs.gyyu,nspred,exo_nbr);
-    dr.ghxuu = unfold12(derivs.gyuu,nspred,exo_nbr);
-    dr.ghuuu = unfold3(derivs.guuu,exo_nbr);
+    dr.ghxxx = derivs.gyyy;
+    dr.ghxxu = derivs.gyyu;
+    dr.ghxuu = derivs.gyuu;
+    dr.ghuuu = derivs.guuu;
     dr.ghxss = derivs.gyss;
     dr.ghuss = derivs.guss;
 else
@@ -145,66 +143,3 @@ else
         dr.ghuu = ghuu;
     end
 end
-
-function y = unfold2(x,n)
-y=zeros(size(x,1),n*n);
-m = 1;
-for i=1:n
-    for j=i:n
-        y(:,(i-1)*n+j)=x(:,m);
-        if j ~= i
-            y(:,(j-1)*n+i)=x(:,m);
-        end
-        m = m+1;
-    end
-end
-
-function y = unfold3(x,n)
-y = zeros(size(x,1),n*n*n);
-m = 1;
-for i=1:n
-    for j=i:n
-        for k=j:n
-            xx = x(:,m);
-            y(:,(i-1)*n*n+(j-1)*n+k) = xx;
-            y(:,(i-1)*n*n+(k-1)*n+j) = xx;
-            y(:,(j-1)*n*n+(k-1)*n+i) = xx;
-            y(:,(j-1)*n*n+(i-1)*n+k) = xx;
-            y(:,(k-1)*n*n+(i-1)*n+j) = xx;
-            y(:,(k-1)*n*n+(j-1)*n+i) = xx;
-            m = m + 1;
-        end
-    end
-end
-
-function y = unfold21(x,n1,n2)
-y = zeros(size(x,1),n1*n1*n2);
-m = 1;
-for i=1:n1
-    for j=i:n1
-        for k=1:n2
-            xx = x(:,m);
-            y(:,(i-1)*n1*n2+(j-1)*n2+k) = xx;
-            if j ~= i
-                y(:,(j-1)*n1*n2+(i-1)*n2+k) = xx;
-            end
-            m = m + 1;
-        end
-    end
-end
-
-function y = unfold12(x,n1,n2)
-y = zeros(size(x,1),n1*n2*n2);
-m = 1;
-for i=1:n1
-    for j=1:n2
-        for k=j:n2
-            xx = x(:,m);
-            y(:,(i-1)*n2*n2+(j-1)*n2+k) = xx;
-            if k ~= j
-                y(:,(i-1)*n2*n2+(k-1)*n2+j) = xx;
-            end
-            m = m + 1;
-        end
-    end
-end
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 19fc3bbfe..27da80e9a 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -59,11 +59,12 @@ DynareMxArrayToString(const mxArray *mxFldp, int len, int width, std::vector<std
 void
 copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &derivs, const std::string &fieldname)
 {
-  const TwoDMatrix &x = derivs.get(sym);
-  int n = x.numRows();
-  int m = x.numCols();
+  const FGSTensor &x = derivs.get(sym);
+  auto x_unfolded = x.unfold();
+  int n = x_unfolded->numRows();
+  int m = x_unfolded->numCols();
   mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
-  std::copy_n(x.getData().base(), n*m, mxGetPr(tmp));
+  std::copy_n(x_unfolded->getData().base(), n*m, mxGetPr(tmp));
   mxSetField(destin, 0, fieldname.c_str(), tmp);
 }
 
@@ -259,6 +260,8 @@ extern "C" {
               }
             if (kOrder == 3 && nlhs > 5)
               {
+                /* Return as 5th argument a struct containing *unfolded* matrices
+                   for 3rd order decision rule */
                 const FGSContainer &derivs = app.get_rule_ders();
                 const std::string fieldnames[] = {"gy", "gu", "gyy", "gyu", "guu", "gss",
                                                   "gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss"};