Merge pull request #52 from McMasterU/improve_interface

Improve interface

HashedExpression.cabal

@@ -4,7 +4,7 @@ cabal-version: 1.12
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: 4917059d4c504182bd983e5ed47c3a6f19b916a0a8a03765addee61e72e7f4a7
+-- hash: 9d0766922c44fc23016fa53343138e21016edc5091f5bb2d2141d69eccd6972c
 
 name:           HashedExpression
 version:        0.0.9
@@ -32,7 +32,9 @@ library
       HashedExpression.Codegen.CSimple
       HashedExpression.Differentiation.Reverse
       HashedExpression.Differentiation.Reverse.State
+      HashedExpression.Dot
       HashedExpression.Embed
+      HashedExpression.Interface
       HashedExpression.Internal
       HashedExpression.Internal.Base
       HashedExpression.Internal.Builder

app/Examples/Brain.hs

@@ -1,34 +1,46 @@
 module Examples.Brain where
 
 import HashedExpression
+import HashedExpression.Modeling.Typed
 import System.FilePath ((</>))
 import Prelude hiding ((**), (^))
-import HashedExpression.Modeling.Typed
 
-brain_reconstructFromMRI :: OptimizationProblem
-brain_reconstructFromMRI =
+brainReconstructFromMRI :: OptimizationProblem
+brainReconstructFromMRI =
   let -- variables
       x = variable2D @128 @128 "x"
+      --- bound
+      xLowerBound = bound2D @128 @128 "x_lb"
+      xUpperBound = bound2D @128 @128 "x_ub"
       -- parameters
       im = param2D @128 @128 "im"
       re = param2D @128 @128 "re"
       mask = param2D @128 @128 "mask"
       -- regularization
+      regularization =
+        norm2square (rotate (0, 1) x - x)
+          + norm2square (rotate (1, 0) x - x)
       lambda = 3000
-      regularization = lambda * (norm2square (rotate (0, 1) x - x) + norm2square (rotate (1, 0) x - x))
    in OptimizationProblem
-        { objective = norm2square ((mask +: 0) * (ft (x +: 0) - (re +: im))) + regularization,
+        { objective =
+            norm2square ((mask +: 0) * (ft (x +: 0) - (re +: im)))
+              + lambda * regularization,
           constraints =
-            [ x .<= VFile (HDF5 "bound.h5" "ub"),
-              x .>= VFile (HDF5 "bound.h5" "lb")
+            [ x .<= xUpperBound,
+              x .>= xLowerBound
             ],
           values =
             [ im :-> VFile (HDF5 "kspace.h5" "im"),
               re :-> VFile (HDF5 "kspace.h5" "re"),
-              mask :-> VFile (HDF5 "mask.h5" "mask")
-            ],
-          workingDir = "examples" </> "Brain"
+              mask :-> VFile (HDF5 "mask.h5" "mask"),
+              xLowerBound :-> VFile (HDF5 "bound.h5" "lb"),
+              xUpperBound :-> VFile (HDF5 "bound.h5" "ub")
+            ]
         }
 
 brain :: IO ()
-brain = proceed brain_reconstructFromMRI CSimpleConfig {output = OutputHDF5, maxIteration = Nothing}
+brain =
+  proceed
+    brainReconstructFromMRI
+    CSimpleConfig {output = OutputHDF5, maxIteration = Nothing}
+    ("examples" </> "Brain")

app/Examples/LinearRegression.hs

@@ -18,9 +18,8 @@ ex1_linearRegression =
           values =
             [ x :-> VFile (TXT "x.txt"),
               y :-> VFile (TXT "y.txt")
-            ],
-          workingDir = "examples" </> "LinearRegression"
+            ]
         }
 
 ex1 :: IO ()
-ex1 = proceed ex1_linearRegression CSimpleConfig {output = OutputText, maxIteration = Nothing}
+ex1 = proceed ex1_linearRegression CSimpleConfig {output = OutputText, maxIteration = Nothing} ("examples" </> "LinearRegression")

app/Examples/LogisticRegression.hs

@@ -1,9 +1,9 @@
 module Examples.LogisticRegression where
 
 import HashedExpression
+import HashedExpression.Modeling.Typed
 import System.FilePath ((</>))
 import Prelude hiding ((**), (^))
-import HashedExpression.Modeling.Typed
 
 sigmoid :: (ToShape d) => TypedExpr d R -> TypedExpr d R
 sigmoid x = 1.0 / (1.0 + exp (- x))
@@ -26,9 +26,12 @@ ex2_logisticRegression =
           values =
             [ x :-> VFile (TXT "x_expanded.txt"),
               y :-> VFile (TXT "y.txt")
-            ],
-          workingDir = "examples" </> "LogisticRegression"
+            ]
         }
 
 ex2 :: IO ()
-ex2 = proceed ex2_logisticRegression CSimpleConfig {output = OutputText, maxIteration = Nothing}
+ex2 =
+  proceed
+    ex2_logisticRegression
+    CSimpleConfig {output = OutputText, maxIteration = Nothing}
+    ("examples" </> "LogisticRegression")

app/Examples/NeuralNetwork.hs

@@ -5,9 +5,9 @@ module Examples.NeuralNetwork where
 import Data.Function ((&))
 import GHC.TypeLits (KnownNat, type (+), type (-))
 import HashedExpression
+import HashedExpression.Modeling.Typed
 import System.FilePath ((</>))
 import Prelude hiding ((**), (^))
-import HashedExpression.Modeling.Typed
 
 sigmoid :: (ToShape d) => TypedExpr d R -> TypedExpr d R
 sigmoid x = 1.0 / (1.0 + exp (- x))
@@ -43,9 +43,12 @@ ex4_neuralNetwork =
           values =
             [ x :-> VFile (HDF5 "data.h5" "x"),
               y :-> VFile (HDF5 "data.h5" "y")
-            ],
-          workingDir = "examples" </> "NeuralNetwork"
+            ]
         }
 
 ex4 :: IO ()
-ex4 = proceed ex4_neuralNetwork CSimpleConfig {output = OutputHDF5, maxIteration = Just 400}
+ex4 =
+  proceed
+    ex4_neuralNetwork
+    CSimpleConfig {output = OutputHDF5, maxIteration = Just 400}
+    ("examples" </> "NeuralNetwork")

app/Main.hs

@@ -9,4 +9,4 @@ import Examples.LogisticRegression
 import Examples.NeuralNetwork
 
 main :: IO ()
-main = ex4
+main = brain 

embed/csimple.c

@@ -8,42 +8,46 @@
 %{fftUtils}
 
 // number of (higher dimensional) variables 
-#define NUM_VARIABLES %{numHigherOrderVariables}
-// number of scalar variables (because each higher dimensional var is a grid of scalar variables)
-#define NUM_ACTUAL_VARIABLES %{numActualVariables}
+#define NUM_HIGH_DIMENSIONAL_VARIABLES %{numHighDimensionalVariables}
+// number of scalar variables (because each high dimensional var is a grid of scalar variables)
+#define NUM_VARIABLES %{numVariables}
 #define MEMORY_NUM_DOUBLES %{totalDoubles}
 #define MEMORY_NUM_COMPLEX_DOUBLES %{totalComplexes}
 // all the actual double variables are allocated one after another, starts from here
 #define VARS_START_OFFSET %{varStartOffset}
 #define MAX_NUM_ITERATIONS %{maxNumIterations}
-#define NUM_SCALAR_CONSTRAINT %{numScalarConstraints}
+#define NUM_GENERAL_CONSTRAINT %{numGeneralConstraints}
 
-const char* var_name[NUM_VARIABLES] = { %{varNames} };
-const int var_size[NUM_VARIABLES] = { %{varSizes} };
+const char* var_name[NUM_HIGH_DIMENSIONAL_VARIABLES] = { %{varNames} };
+const int var_size[NUM_HIGH_DIMENSIONAL_VARIABLES] = { %{varSizes} };
 
-const int var_offset[NUM_VARIABLES] = { %{varOffsets} };
-const int partial_derivative_offset[NUM_VARIABLES] = { %{partialDerivativeOffsets} };
+const int var_offset[NUM_HIGH_DIMENSIONAL_VARIABLES] = { %{varOffsets} };
+const int partial_derivative_offset[NUM_HIGH_DIMENSIONAL_VARIABLES] = { %{partialDerivativeOffsets} };
 const int objective_offset = %{objectiveOffset};
 
 double ptr[MEMORY_NUM_DOUBLES];
 complex double ptr_c[MEMORY_NUM_COMPLEX_DOUBLES];
 
-double lower_bound[NUM_ACTUAL_VARIABLES];
-double upper_bound[NUM_ACTUAL_VARIABLES];
+double lower_bound[NUM_VARIABLES];
+double upper_bound[NUM_VARIABLES];
 
-double sc_lower_bound[NUM_SCALAR_CONSTRAINT];
-double sc_upper_bound[NUM_SCALAR_CONSTRAINT];
-const int sc_offset[NUM_SCALAR_CONSTRAINT] = { %{scalarConstraintOffsets} };
+double sc_lower_bound[NUM_GENERAL_CONSTRAINT];
+double sc_upper_bound[NUM_GENERAL_CONSTRAINT];
+const int sc_offset[NUM_GENERAL_CONSTRAINT] = { %{generalConstraintOffsets} };
 
-const int sc_partial_derivative_offset[NUM_SCALAR_CONSTRAINT][NUM_VARIABLES] = {  %{scalarConstraintPartialDerivativeOffsets} };
+const int sc_partial_derivative_offset[NUM_GENERAL_CONSTRAINT][NUM_HIGH_DIMENSIONAL_VARIABLES] = {  %{generalConstraintPartialDerivativeOffsets} };
 
 void read_bounds() {
-  for (int i = 0; i < NUM_ACTUAL_VARIABLES; i++) {
+  for (int i = 0; i < NUM_VARIABLES; i++) {
     lower_bound[i] = -INFINITY;
     upper_bound[i] = INFINITY;
   }
+  for (int i = 0; i < NUM_GENERAL_CONSTRAINT; i++) {
+    sc_lower_bound[i] = -INFINITY;
+    sc_upper_bound[i] = INFINITY;
+  }
   %{readBounds}
-  %{readBoundScalarConstraints}
+  %{readBoundGeneralConstraints}
 }
 
 void read_values() {
@@ -70,4 +74,4 @@ void evaluate_scalar_constraints() {
 }
 void evaluate_scalar_constraints_jacobian() {
   %{evaluateScalarConstraintsJacobian}
-}
+}

examples/brain/plot.py

@@ -19,6 +19,10 @@ def plot_image(data):
 
 plt.ion()
 
+print('Real part acquired by MRI...')
+bound = read_hdf5("bound.h5", "lb")
+plot_image(bound)
+input('Program paused. Press ENTER to continue')
 
 print('Real part acquired by MRI...')
 re = read_hdf5("kspace.h5", "re")
@@ -30,6 +34,11 @@ def plot_image(data):
 plot_image(np.log(np.abs(im) + 1e-10))
 input('Program paused. Press ENTER to continue')
 
+print('Mask...')
+mask = read_hdf5("mask.h5", "mask")
+plot_image(mask)
+input('Program paused. Press ENTER to continue')
+
 print('Naively reconstruct by taking inverse Fourier Transform...')
 naive = np.abs(np.fft.ifft2(re + 1j * im))
 plot_image(naive)

solvers/gradient_descent/gradient_descent.c

@@ -8,12 +8,12 @@
 
 #define oo 10000000
 
-extern const char* var_name[NUM_VARIABLES];
-extern const int var_num_dim[NUM_VARIABLES];
-extern const int var_shape[NUM_VARIABLES][3];
-extern const int var_size[NUM_VARIABLES];
-extern const int var_offset[NUM_VARIABLES];
-extern const int partial_derivative_offset[NUM_VARIABLES];
+extern const char* var_name[NUM_HIGH_DIMENSIONAL_VARIABLES];
+extern const int var_num_dim[NUM_HIGH_DIMENSIONAL_VARIABLES];
+extern const int var_shape[NUM_HIGH_DIMENSIONAL_VARIABLES][3];
+extern const int var_size[NUM_HIGH_DIMENSIONAL_VARIABLES];
+extern const int var_offset[NUM_HIGH_DIMENSIONAL_VARIABLES];
+extern const int partial_derivative_offset[NUM_HIGH_DIMENSIONAL_VARIABLES];
 extern const int objective_offset;
 extern double ptr[MEM_SIZE];
 
@@ -34,7 +34,7 @@ double random_in(double min, double max) {
 
 bool any_partial_derivative_NaN() {
   int i, j;
-  for (i = 0; i < NUM_VARIABLES; i++) {
+  for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
     for (j = 0; j < var_size[i]; j++) {
       if (isnan(ptr[partial_derivative_offset[i] + j])) {
         return true;
@@ -49,7 +49,7 @@ int main() {
   srand(time(NULL));
   int i, j, cnt, k;
   int total_variable_size = 0;
-  for (i = 0; i < NUM_VARIABLES; i++) {
+  for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
     total_variable_size = total_variable_size + var_size[i];
   }
 
@@ -61,7 +61,7 @@ int main() {
 
   while (true) {
     // initialize optimizing variables
-    for (i = 0; i < NUM_VARIABLES; i++) {
+    for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
       for (j = 0; j < var_size[i]; j++) {
         ptr[var_offset[i] + j] = random_in(0, 1);
       }
@@ -77,7 +77,7 @@ int main() {
 
       // save the state of all variables
       cnt = 0;
-      for (i = 0; i < NUM_VARIABLES; i++) {
+      for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
         for (j = 0; j < var_size[i]; j++) {
           var_temp[cnt] = ptr[var_offset[i] + j];
           grad_temp[cnt] = ptr[partial_derivative_offset[i] + j];
@@ -88,7 +88,7 @@ int main() {
       // fx
       double fx = ptr[objective_offset];
       double minus_dot_grad = 0;
-      for (i = 0; i < NUM_VARIABLES; i++) {
+      for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
         for (j = 0; j < var_size[i]; j++) {
           minus_dot_grad = minus_dot_grad - (ptr[partial_derivative_offset[i] + j] * ptr[partial_derivative_offset[i] + j]);
         }
@@ -102,7 +102,7 @@ int main() {
       // after while, x still the same
       while (true) {
         // write updated variables
-        for (i = 0; i < NUM_VARIABLES; i++) {
+        for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
           for (j = 0; j < var_size[i]; j++) {
             ptr[var_offset[i] + j] -= t * ptr[partial_derivative_offset[i] + j];
           }
@@ -111,7 +111,7 @@ int main() {
         double fnew = ptr[objective_offset];
 
         cnt = 0;
-        for (i = 0; i < NUM_VARIABLES; i++) {
+        for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
           for (j = 0; j < var_size[i]; j++) {
             ptr[var_offset[i] + j] = var_temp[cnt];
             cnt++;
@@ -124,7 +124,7 @@ int main() {
           // acc = - grad(old) <.> grad(new)
           double acc = 0;
           cnt = 0;
-          for (i = 0; i < NUM_VARIABLES; i++) {
+          for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
             for (j = 0; j < var_size[i]; j++) {
               acc = acc - grad_temp[cnt] * ptr[partial_derivative_offset[i] + j];
               cnt++;
@@ -154,23 +154,23 @@ int main() {
       while (t > 0 && any_partial_derivative_NaN()) {
         t = t / 2;
         cnt = 0;
-        for (i = 0; i < NUM_VARIABLES; i++) {
+        for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
           for (j = 0; j < var_size[i]; j++) {
             var_temp[cnt] = ptr[var_offset[i] + j];
             cnt++;
           }
         }
 
         cnt = 0;
-        for (i = 0; i < NUM_VARIABLES; i++) {
+        for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
           for (j = 0; j < var_size[i]; j++) {
             ptr[var_offset[i] + j] -= t * grad_temp[cnt];
             cnt++;
           }
         }
         evaluate_partial_derivatives_and_objective();
         cnt = 0;
-        for (i = 0; i < NUM_VARIABLES; i++) {
+        for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
           for (j = 0; j < var_size[i]; j++) {
             ptr[var_offset[i] + j] = var_temp[cnt];
             cnt++;
@@ -187,7 +187,7 @@ int main() {
       double max_step = 0;
       // OK now we have a good t, let's update
       cnt = 0;
-      for (i = 0; i < NUM_VARIABLES; i++) {
+      for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
         for (j = 0; j < var_size[i]; j++) {
           double step = -t * grad_temp[cnt];
           ptr[var_offset[i] + j] += step;
@@ -209,7 +209,7 @@ int main() {
 
     if (!isnan(ptr[objective_offset])) {
       printf("f_min = %f\n", ptr[objective_offset]);
-      for (i = 0; i < NUM_VARIABLES; i++) {
+      for (i = 0; i < NUM_HIGH_DIMENSIONAL_VARIABLES; i++) {
         char* var_file_name = malloc(strlen(var_name[i]) + 4);
         strcpy(var_file_name, var_name[i]);
         strcat(var_file_name, ".txt");
@@ -248,4 +248,4 @@ int main() {
 
   free(var_temp);
   free(grad_temp);
-}
+}