Everything builds after re-doing to use ibmisc/linear

Author: Elizabeth Fischer

CMakeLists.txt

@@ -102,6 +102,10 @@ find_package(Ibmisc REQUIRED)
 list(APPEND EXTERNAL_LIBS ${IBMISC_LIBRARY})
 include_directories(${IBMISC_INCLUDE_DIR})
 
+find_package(ZLIB REQUIRED)
+list(APPEND EXTERNAL_LIBS ${ZLIB_LIBRARIES})
+include_directories(${ZLIB_INCLUDE_DIRS})
+
 
 # --- Proj.4
 find_package(PROJ4 REQUIRED)

modele/combine_global_ec.cpp

@@ -3,7 +3,8 @@
 #include <ibmisc/netcdf.hpp>
 #include <ibmisc/runlength.hpp>
 #include <icebin/error.hpp>
-#include <icebin/RLWeightedSparse.hpp>
+#include <ibmisc/linear/linear.hpp>
+#include <ibmisc/linear/compressed.hpp>
 
 using namespace std;
 using namespace ibmisc;
@@ -62,11 +63,11 @@ void combine_chunks(
 
 
     // Allocate
-    RLWeightedSparse ret(sparse_extents);
+    linear::Weighted_Compressed ret;
 
-    {auto wM(ret.wM.accum());
+    {auto wM(ret.weights[0].accum());
     auto M(ret.M.accum());
-    auto Mw(ret.Mw.accum());
+    auto Mw(ret.weights[1].accum());
 
         // Aggregate them together
         int iM=0;
@@ -78,16 +79,16 @@ void combine_chunks(
             auto dimA(nc_read_blitz<int,1>(ncio.nc, "dim"+sgrids[1]));
 
             auto wM_d(nc_read_blitz<double,1>(ncio.nc, BvA+".wM"));
-            for (int i=0; i<wM_d.extent(0); ++i) wM->add({dimB(i)}, wM_d(i));
+            for (int i=0; i<wM_d.extent(0); ++i) wM.add({dimB(i)}, wM_d(i));
 
             auto indices_d(nc_read_blitz<int,2>(ncio.nc, BvA+".M.indices"));
             auto values_d(nc_read_blitz<double,1>(ncio.nc, BvA+".M.values"));
             for (int i=0; i<values_d.extent(0); ++i) {
-                M->add({dimB(indices_d(i,0)), dimA(indices_d(i,1))}, values_d(i));
+                M.add({dimB(indices_d(i,0)), dimA(indices_d(i,1))}, values_d(i));
             }
 
             auto Mw_d(nc_read_blitz<double,1>(ncio.nc, BvA+".Mw"));
-            for (int i=0; i<Mw_d.extent(0); ++i) Mw->add({dimA(i)}, Mw_d(i));
+            for (int i=0; i<Mw_d.extent(0); ++i) Mw.add({dimA(i)}, Mw_d(i));
         }
 
         // Check

pylib/_icebin.pyx

@@ -53,8 +53,7 @@ cdef class RegridMatrices:
     def __dealloc__(self):
         del self.cself
 
-    def matrix(self, str spec_name, bool scale=True, bool correctA=True,
-        sigma=(0,0,0), conserve=True):
+    def matrix(self, str spec_name):
         """Compute a regrid matrix.
         spec_name:
             Type of regrid matrix to obtain.  Choice are:
@@ -68,12 +67,11 @@ cdef class RegridMatrices:
             for the output grid).
         returns: WeightedSparse
         """
-        cdef cicebin.CythonWeightedSparse *crm
-        crm = cicebin.RegridMatrices_matrix(
-            self.cself, spec_name.encode(), scale, correctA,
-            sigma[0], sigma[1], sigma[2], conserve)
-        ret = WeightedSparse()
-        ret.cself = crm
+        cdef cibmisc.linear_Weighted *lw
+        lw = cicebin.RegridMatrices_matrix(self.cself, spec_name.encode())
+        cdef ibmisc.linear_Weighted ret
+        ret = ibmisc.linear_Weighted()
+        ret.cself = lw
         return ret
 
 cdef class GCMRegridder:
@@ -167,33 +165,19 @@ cdef class GCMRegridder:
             exgrid_fname.encode(), exgrid_vname.encode(),
             interp_style.encode())
 
-    def regrid_matrices(self, str sheet_name, elevmaskI):
+    def regrid_matrices(self, str sheet_name, elevmaskI,
+        bool scale=True, bool correctA=True,
+        sigma=(0,0,0), conserve=True):
+
         elevmaskI = elevmaskI.reshape(-1)
         cdef cicebin.RegridMatrices *crm = \
             cicebin.new_regrid_matrices(self.cself.get(), sheet_name.encode(),
-            <PyObject *>elevmaskI)   # PyObject=Borrowed reference, object = owned reference
+            <PyObject *>elevmaskI,   # PyObject=Borrowed reference, object = owned reference
+            scale, correctA, sigma[0], sigma[1], sigma[2], conserve)
         rm = RegridMatrices()
         rm.cself = crm
         return rm
 
-def coo_multiply(M, xx, double fill=np.nan, ignore_nan=False):
-    """M:
-        SciPy sparse matrix"""
-
-    warnings.warn(
-        "coo_multiply is deprecated; use WeightedSparse.apply() instead.",
-        DeprecationWarning)
-
-    xx = xx.reshape(-1)
-    yy = np.zeros(M._shape[0])
-    yy[:] = fill
-
-    cicebin.coo_matvec(<PyObject *>yy, <PyObject *>xx, ignore_nan,
-        M._shape[0], M._shape[1],
-        <PyObject *>M.row, <PyObject *>M.col, <PyObject *>M.data)
-
-    return yy
-
 # ============================================================
 
 cdef class HntrSpec:

pylib/cicebin.pxd

@@ -74,9 +74,6 @@ cdef extern from "icebin/GCMRegridder.hpp" namespace "icebin":
         GCMRegridder_Standard() except +
 
 cdef extern from "icebin_cython.hpp" namespace "icebin::cython":
-    cdef cppclass CythonWeightedSparse:
-        object shape() except +
-
     cdef void read_fgrid(
         cibmisc.unique_ptr[Grid] &fgridA,
         string &gridA_fname,
@@ -108,22 +105,14 @@ cdef extern from "icebin_cython.hpp" namespace "icebin::cython":
         string &exgrid_fname, string &exgrid_vname,
         string &sinterp_style) except +
 
-    cdef CythonWeightedSparse *RegridMatrices_matrix(
-        RegridMatrices *self, string spec_name, bool scale, bool correctA,
-        double sigma_x, double sigma_y, double sigma_z, bool conserve) except +
-
-    cdef object CythonWeightedSparse_dense_extent(CythonWeightedSparse *self) except +
-    cdef object CythonWeightedSparse_sparse_extent(CythonWeightedSparse *self) except +
-    cdef object CythonWeightedSparse_to_tuple(CythonWeightedSparse *self) except +
-
-    cdef object CythonWeightedSparse_apply(CythonWeightedSparse *BvA, PyObject *A, double fill, bool force_conservationlmake) except +
-
-    cdef void coo_matvec(PyObject *yy_py, PyObject *xx_py, bool ignore_nan,
-        int M_nrow, int M_ncol, PyObject *M_row_py, PyObject *M_col_py, PyObject *M_data_py) except +
+    cdef cibmisc.linear_Weighted *RegridMatrices_matrix(
+        RegridMatrices *self, string spec_name) except +
 
     cdef Hntr_regrid(Hntr *hntr, object WTA_py, object A_py, bool mean_polar) except +
 
-    cdef RegridMatrices *new_regrid_matrices(GCMRegridder *gcm, string &sheet_name, PyObject *elevmaskI_py) except +
+    cdef RegridMatrices *new_regrid_matrices(GCMRegridder *gcm, string &sheet_name, PyObject *elevmaskI_py,
+        bool scale, bool correctA,
+        double sigma_x, double sigma_y, double sigma_z, bool conserve) except +
 
     cdef void update_topo(GCMRegridder *_gcmA, string &topoO_fname,
         PyObject *elev_sigmas_py, bool initial_timestep, string &segments_py, string &primary_segment_py,

pylib/icebin_cython.cpp

@@ -185,129 +185,11 @@ void coo_matvec(PyObject *yy_py, PyObject *xx_py, bool ignore_nan,
 }
 
 
-extern CythonWeightedSparse *RegridMatrices_matrix(RegridMatrices *cself, std::string const &spec_name, bool scale, bool correctA, double sigma_x, double sigma_y, double sigma_z, bool conserve)
+extern linear::Weighted *RegridMatrices_matrix(RegridMatrices *cself, std::string const &spec_name)
 {
-    std::unique_ptr<CythonWeightedSparse> CRM(new CythonWeightedSparse());
-
-    CRM->RM = cself->matrix_d(spec_name,
-        {&CRM->dims[0], &CRM->dims[1]},
-        RegridParams(scale, correctA, {sigma_x, sigma_y, sigma_z}));
-
-    CythonWeightedSparse *ret = CRM.release();
-    return ret;
+    return cself->matrix(spec_name).release();
 }
 
-extern PyObject *CythonWeightedSparse_apply(
-    CythonWeightedSparse *BvA,
-    PyObject *A_s_py,            // A_b{nj_s} One row per variable
-    double fill, bool force_conservation)
-{
-    // |j_s| = size of sparse input vector space (A_s)
-    // |j_d] = size of dense input vector space (A_d)
-    // |n| = number of variables being processed together
-
-    // Allocate dense A matrix
-    auto &bdim(BvA->dims[0]);
-    auto &adim(BvA->dims[1]);
-    auto A_s(np_to_blitz<double,2>(A_s_py, "A", {-1,-1}));    // Sparse indexed dense vector
-    int n_n = A_s.extent(0);
-
-    // Densify the A matrix
-    blitz::Array<double,2> A_d(n_n, adim.dense_extent());
-    for (int j_d=0; j_d < adim.dense_extent(); ++j_d) {
-        int j_s = adim.to_sparse(j_d);
-        for (int n=0; n < n_n; ++n) {
-            A_d(n,j_d) = A_s(n,j_s);
-        }
-    }
-
-    // Apply...
-    auto B_d_eigen(BvA->RM->apply_e(A_d, fill, force_conservation));    // Column major indexing
-
-    // Allocate output vector and get a Blitz view
-    // We will copy from the Eigen data structure to the Python
-    PyObject *B_s_py = ibmisc::cython::new_pyarray<double,2>(
-        std::array<int,2>{n_n, bdim.sparse_extent()});
-    auto B_s(np_to_blitz<double,2>(B_s_py, "B_s_py", {-1,-1}));
-
-    // Sparsify the output B
-    for (int n=0; n < n_n; ++n)
-    for (int j_s=0; j_s < bdim.sparse_extent(); ++j_s) {
-        B_s(n,j_s) = fill;
-    }
-
-    for (int j_d=0; j_d < bdim.dense_extent(); ++j_d) {
-        int j_s = bdim.to_sparse(j_d);
-        for (int n=0; n < n_n; ++n) {
-            B_s(n,j_s) = B_d_eigen(j_d,n);
-        }
-    }
-
-    return B_s_py;
-}
-
-PyObject *CythonWeightedSparse_dense_extent(CythonWeightedSparse const *cself)
-{
-    return Py_BuildValue("ll",
-        (long)cself->dims[0].dense_extent(),
-        (long)cself->dims[1].dense_extent());
-}
-
-PyObject *CythonWeightedSparse_sparse_extent(CythonWeightedSparse const *cself)
-{
-    return Py_BuildValue("ll",
-        (long)cself->dims[0].sparse_extent(),
-        (long)cself->dims[1].sparse_extent());
-}
-
-
-PyObject *CythonWeightedSparse_to_tuple(CythonWeightedSparse *cself)
-{
-    // ----- Convert a dense vector w/ dense indices to a dense vector with sparse indices
-    if (cself->dims[0].sparse_extent() <= 0) (*icebin_error)(-1,
-        "dims[0].sparse_extent() must be > 0");
-
-    // ---------- wM
-    // Allocate the output Python vector
-    PyObject *wM_py = ibmisc::cython::new_pyarray<double,1>(
-        std::array<int,1>{cself->dims[0].sparse_extent()});
-    // Get a Blitz view on it
-    auto wM_pyb(np_to_blitz<double,1>(wM_py, "wM_py", {-1}));
-    // Copy, while translating the dimension
-    spcopy(
-        accum::to_sparse(make_array(&cself->dims[0]),
-        accum::blitz_existing(wM_pyb)),
-        cself->RM->wM);
-
-    // -------------- M
-    // Convert a sparse matrix w/ dense indices to a sparse matrix with sparse indices
-    TupleListT<2> RM_sp;
-printf("DD1 about to copy: sparse_extent=(%ld, %ld), dense_extent=(%d, %d)\n",
-cself->dims[0].sparse_extent(),
-cself->dims[1].sparse_extent(),
-cself->dims[0].dense_extent(),
-cself->dims[1].dense_extent());
-
-    spcopy(
-        accum::to_sparse(make_array(&cself->dims[0], &cself->dims[1]),
-        accum::ref(RM_sp)),
-        *cself->RM->M);
-    PyObject *M_py = ibmisc::cython::spsparse_to_tuple(RM_sp);
-
-    // ---------- Mw
-    // Allocate the output Python vector
-    PyObject *Mw_py = ibmisc::cython::new_pyarray<double,1>(
-        std::array<int,1>{cself->dims[1].sparse_extent()});
-    // Get a Blitz view on it
-    auto Mw_pyb(np_to_blitz<double,1>(Mw_py, "Mw_py", {-1}));
-    // Copy, while translating the dimension
-    spcopy(
-        accum::to_sparse(make_array(&cself->dims[1]),
-        accum::blitz_existing(Mw_pyb)),
-        cself->RM->Mw);
-
-    return Py_BuildValue("OOO", wM_py, M_py, Mw_py);
-}
 
 // ------------------------------------------------------------
 PyObject *Hntr_regrid(Hntr const *hntr, PyObject *WTA_py, PyObject *A_py, bool mean_polar)
@@ -323,13 +205,27 @@ PyObject *Hntr_regrid(Hntr const *hntr, PyObject *WTA_py, PyObject *A_py, bool m
     return B_py;
 }
 
-RegridMatrices *new_regrid_matrices(GCMRegridder const *gcm, std::string const &sheet_name, PyObject *elevmaskI_py)
+RegridMatrices *new_regrid_matrices(
+    GCMRegridder const *gcm,
+    std::string const &sheet_name,
+    PyObject *elevmaskI_py,
+    // --------- Params
+    bool scale,
+    bool correctA,
+    double sigma_x,
+    double sigma_y,
+    double sigma_z,
+    bool conserve)
 {
+
+
     auto sheet_index = gcm->ice_regridders().index.at(sheet_name);
     IceRegridder *ice_regridder = &*gcm->ice_regridders()[sheet_index];
     auto elevmaskI(np_to_blitz<double,1>(elevmaskI_py, "elevmaskI", {ice_regridder->nI()}));
 
-    return new RegridMatrices(gcm->regrid_matrices(sheet_index, elevmaskI));
+    return gcm->regrid_matrices(
+        sheet_index, elevmaskI,
+        RegridParams(scale, correctA, {sigma_x, sigma_y, sigma_z})).release();
 }
 
 std::string to_string(PyObject *str, std::string const &vname)

pylib/icebin_cython.hpp

@@ -65,39 +65,24 @@ extern void GCMRegridder_add_sheet(GCMRegridder *cself,
     std::string const &sinterp_style);
 
 
-/** Wraps WeightedSparse to keep around the dense/sparse dimension
-    translators */
-struct CythonWeightedSparse {
-    std::array<SparseSetT,2> dims;
-    std::unique_ptr<WeightedSparse> RM;
+extern ibmisc::linear::Weighted *RegridMatrices_matrix(RegridMatrices *cself,
+    std::string const &spec_name);
 
-    /** @return 2-D shape of the stored matrix */
-    PyObject *shape();
-};
-
-
-extern CythonWeightedSparse *RegridMatrices_matrix(RegridMatrices *cself,
-    std::string const &spec_name, bool scale, bool correctA,
-    double sigma_x, double sigma_y, double sigma_z, bool conserve);
-
-extern PyObject *CythonWeightedSparse_apply(
-    CythonWeightedSparse *BvA,
-    PyObject *A_s_py,
-    double fill, bool force_conservation);            // A_b{nj_s} One row per variable
-
-
-PyObject *CythonWeightedSparse_dense_extent(CythonWeightedSparse const *cself);
-PyObject *CythonWeightedSparse_sparse_extent(CythonWeightedSparse const *cself);
-
-PyObject *CythonWeightedSparse_to_tuple(CythonWeightedSparse *cself);
-
-void coo_matvec(PyObject *yy_py, PyObject *xx_py, bool ignore_nan,
-    size_t M_nrow, size_t M_ncol, PyObject *M_row_py, PyObject *M_col_py, PyObject *M_data_py);
 
 PyObject *Hntr_regrid(modele::Hntr const *hntr, PyObject *WTA_py, PyObject *A_py, bool mean_polar);
 
 
-RegridMatrices *new_regrid_matrices(GCMRegridder const *gcm, std::string const &sheet_name, PyObject *elevmaskI_py);
+RegridMatrices *new_regrid_matrices(
+    GCMRegridder const *gcm,
+    std::string const &sheet_name,
+    PyObject *elevmaskI_py,
+    // --------- Params
+    bool scale,
+    bool correctA,
+    double sigma_x,
+    double sigma_y,
+    double sigma_z,
+    bool conserve);
 
 /** Allows Python users access to GCMCoupler_Modele::update_topo().
 Starting from output of Gary's program (on the Ocean grid), this subroutine

slib/icebin/IceCoupler.cpp

@@ -475,10 +475,10 @@ bool run_ice)
         dimA1.ncio(ncio, "dimA");
         dimE1.ncio(ncio, "dimE");
 
-        AvE1->ncio_nodim(ncio, "AvE");
-        E1vI_nc->ncio_nodim(ncio, "EvI_nc");
-        A1vI->ncio_nodim(ncio, "AvI");
-        IvE1->ncio_nodim(ncio, "IvE");
+        AvE1->ncio(ncio, "AvE", {"dimA", "dimE"});
+        E1vI_nc->ncio(ncio, "EvI_nc", {"dimE", "dimI"});
+        A1vI->ncio(ncio, "AvI", {"dimA", "dimI"});
+        IvE1->ncio(ncio, "IvE", {"dimI", "dimE"});
     }
 
     // Save stuff for next time around