Merge pull request #17 from fruzsinaagocs/fix_reverse_integration

Fixing code s.t. unit tests for reverse integration are passed again
fruzsinaagocs · Feb 15, 2022 · 22dcfaf · 22dcfaf
2 parents d008f63 + 0b94c23
commit 22dcfaf
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 26 deletions.
diff --git a/include/solver.hpp b/include/solver.hpp
@@ -217,7 +217,7 @@ a_tol, double h_0, const char* full_output){
             sign = 1;
     }
     else{
-        throw "Direction of integration in conflict with direction of initial step, terminating. Please check your values for ti, tf, and h. ";
+        throw std::logic_error("Direction of integration in conflict with direction of initial step, terminating. Please check your values for ti, tf, and h. ");
         return;
     }
 
@@ -226,20 +226,21 @@ a_tol, double h_0, const char* full_output){
     dotimes.resize(dosize);
     dosol.resize(dosize);
     dodsol.resize(dosize);
-    int docount = 0;
     auto doit = do_times.begin();
-    if(de_sys_->Winterp.sign_ == 1 or sign == 1){
+
+    if((de_sys_->is_interpolated == 1 and de_sys_->Winterp.sign_ == 1) or (de_sys_->is_interpolated == 0 and sign == 1)){
         for(auto it=dotimes.begin(); it!=dotimes.end(); it++){
             *it = *doit;
-            docount++; doit++;
+            doit++;
         }
     }
     else{
-         for(auto it=dotimes.rbegin(); it!=dotimes.rend(); ++it){
+        for(auto it=dotimes.rbegin(); it!=dotimes.rend(); ++it){
             *it = *doit;
-            docount++; ++doit;
+            ++doit;
         }
     }
+
     dotit = dotimes.begin();
     switch(order){
         case 1: wkbsolver1 = WKBSolver1(*de_sys_, order);
@@ -470,21 +471,6 @@ void Solution::solve(){
         };
     };
 
-    // If integrating backwards, reverse dense output (because it will have been
-    // reversed at the start)
-    if(de_sys_->is_interpolated == 1){
-        if(de_sys_->Winterp.sign_ == 0){
-            dosol.reverse();
-            dodsol.reverse();
-        }
-    }
-    else{
-        if(sign == 0){
-            dosol.reverse();
-            dodsol.reverse();
-        }
-    }
-
     // Write output to file if prompted
     if(not (*fo==0)){
         std::string output(fo);
@@ -520,5 +506,21 @@ void Solution::solve(){
 
         f.close();
     }
-
+
+    // If integrating backwards, reverse dense output (because it will have been
+    // reversed at the start)
+    if(de_sys_->is_interpolated == 1){
+        if(de_sys_->Winterp.sign_ == 0){
+            dosol.reverse();
+            dodsol.reverse();
+        }
+    }
+    else{
+        if(sign == 0){
+            dosol.reverse();
+            dodsol.reverse();
+        }
+    }
+
+
 };
diff --git a/setup.py b/setup.py
@@ -20,7 +20,7 @@ def readme(short=False):
 
 setup(
     name="pyoscode",
-    version="1.0.4",
+    version="1.0.5",
     description=readme(short=True),
     long_description=readme(),
     url="https://github.com/fruzsinaagocs/oscode",

diff --git a/tests/test_airy.py b/tests/test_airy.py
@@ -113,7 +113,7 @@ def test_airy_backward_uneven():
     dense = np.asarray(sol['x_eval'])
     dense_d = np.asarray(sol['dx_eval'])   
     types = np.asarray(sol['types'])
-    ana_t = np.linspace(ti,tf,5000)
+    ana_t = t_eval
     ana_x = np.asarray([airy(-T)[0]+1j*airy(-T)[2] for T in t])
     dense_ana_x = np.asarray([airy(-T)[0]+1j*airy(-T)[2] for T in ana_t])
     dense_ana_dx = np.asarray([-airy(-T)[1]-1j*airy(-T)[3] for T in ana_t])
@@ -123,5 +123,3 @@ def test_airy_backward_uneven():
     assert (np.any(dense_error) == False and np.any(dense_error_d) == False and np.any(error) == False )
 
 
-
-