Additional low-rank decomposition methods (by Matt, on top of a few uncommitted updates on my side)

.gitignore

@@ -8,11 +8,13 @@
 *.npz
 *.obj
 *.pdf
+*.ply
 *.png
 *.so
 *.tif
 *.txt
 *.zip
+*.mat
 *~
 .DS_Store
 .DS_Store
@@ -27,4 +29,6 @@ dist
 src/flux/cgal/aabb.cpp
 src/flux/thermal.c
 src/flux/view_factor/view_factor.c
+src/flux/cgal/CGAL/
+src/flux/cgal/boost/
 stats

README.md

@@ -4,38 +4,52 @@
 
 ## Installation ##
 
-### Set up a virtual environment and clone the repository ###
+If using Windows, first install Microsoft Visual Studio C++ Build Tools following 
+the [directions at this link](https://github.com/bycloudai/InstallVSBuildToolsWindows). This includes
+an installation of MSVC, Windows SDK, and C++ CMake tools for Windows followed by adding MSBuild Tools
+to your system path.
 
-Make a new directory and clone this repository to it. Then, inside the
-directory that you've just created, run `python -m venv .`. This will
-create a "virtual environment", which is a useful tool that Python
-provides for managing dependencies for projects. The new directory
-"contains" the virtual environment.
+### Setting up a virtual environment ###
 
-### Activate the virtual environment ###
-
-To activate the virtual environment, run `source bin/activate` from
-the directory containing it. Make sure to do this before doing
-anything else below.
+Using [Anaconda](https://www.anaconda.com/),
+a virtual environment management tool for Python, create and activate a new conda environment and download
+pip via the anaconda channel by running:
+``` shell
+conda create --name radiosity
+conda activate radiosity
+conda install -c anaconda pip
+conda install -c anaconda cython
+conda install -c conda-forge rioxarray
+```
+Alternatively, run `python -m venv .` to create a virtual environment in standard Python without using Anaconda. Activating the virtual 
+environment (via `conda activate radiosity` or `source bin/activate`) means that all python and pip operations will be
+managed using the packages installed into the environment. Therefore, the command
+must be run before proceeding to the next steps (and before running any code).
 
 ### Getting the dependencies ###
 
 First, install
 [python-embree](https://github.com/sampotter/python-embree) in this
-virtual environment:
+conda environment. Begin by installing Embree from the [Embree website](https://www.embree.org/). Typically, the Embree binaries,
+headers, and libraries will all be installed to C:\Program Files\Intel\Embree3 by default. Then, clone, build, and install python-embree
+to the radiosity conda environment:
 ``` shell
 git clone https://github.com/sampotter/python-embree
 cd python-embree
+python setup.py build_ext -I "/c/Program\ Files/Intel/Embree3/include" -L "/c/Program\ Files/Intel/Embree3/lib"
 python setup.py install
-cd -
+cd ..
 ```
 
-Next, install the rest of the dependencies:
+Next, clone this repository and install the rest of the dependencies, including Boost and CGAL.
 ``` shell
 git clone https://github.com/sampotter/python-flux
 cd python-flux
 pip install -r requirements.txt
 ```
+Install the latest versions of [Boost](https://www.boost.org/) 
+and [CGAL](https://www.cgal.org/download.html) from the corresponding websites. Finally,
+link the CGAL and boost directories in `python-flux/setup.py` (via the `include_dirs` argument of the `flux.cgal.aabb` extension).
 
 ### Installing this package ###
 
@@ -44,24 +58,23 @@ Finally, from the cloned repository, run:
 python setup.py install
 ```
 
-## Running the unit tests ##
-
-To run python-flux's unit tests, just run:
-``` shell
-./run_tests.sh
-```
-from the root of this repository. All this script does is execute `python -m unittest discover ./tests`; i.e., it discovers and runs all Python unit tests found in the directory `./tests`.
-
 ## Running the examples ##
 
 The `examples` directory contains simple Python scripts illustrating
-how to use this package. Each example directory contains a README with
-more information about that particulra example.
+how to use this package. The most up-to-date implementations can be
+found in the `examples/shackleton_vary_outer` and `examples/blurred_south_pole`
+sub-directories.
 
 For instance, try running:
 ``` shell
-cd examples/lunar_south_pole
-python haworth.py
+conda install -c conda-forge arrow
+conda install -c conda-forge tqdm
+
+cd examples/shackleton_vary_outer
+python make_mesh.py 5.0 40
+python make_compressed_form_factor_matrix.py --compression_type "svd" --max_area 5.0 --outer_radius 40 --tol 1e-1 --k0 40 --add_residuals
+python make_block_plots.py --compression_type "svd" --max_area 5.0 --outer_radius 40 --tol 1e-1 --k0 40 --add_residuals
+python collect_time_dep_data_nomemory.py --compression_type "svd" --max_area 5.0 --outer_radius 40 --tol 1e-1 --k0 40 --add_residuals
 ```
 after following the installation instructions
 above.

examples/blurred_south_pole/make_block_plots.py

@@ -0,0 +1,175 @@
+import itertools as it
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+import numpy as np
+
+import flux.compressed_form_factors as cff
+from flux.compressed_form_factors import CompressedFormFactorMatrix
+
+import scipy
+import os
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--compression_type', type=str, default="svd", choices=["svd", "rand_svd", "aca", "rand_id"])
+parser.add_argument('--max_area', type=float, default=3.0)
+parser.add_argument('--outer_radius', type=int, default=80)
+parser.add_argument('--sigma', type=float, default=5.0)
+
+parser.add_argument('--tol', type=float, default=1e-1)
+parser.add_argument('--min_depth', type=int, default=1)
+parser.add_argument('--max_depth', type=int, default=0)
+parser.add_argument('--compress_sparse', action='store_true')
+
+parser.add_argument('--add_residuals', action='store_true')
+parser.add_argument('--k0', type=int, default=40)
+parser.add_argument('--p', type=int, default=5)
+parser.add_argument('--q', type=int, default=1)
+
+parser.add_argument('--cliques', action='store_true')
+parser.add_argument('--n_cliques', type=int, default=25)
+parser.add_argument('--obb', action='store_true')
+
+parser.add_argument('--plot_labels', action='store_true')
+
+parser.set_defaults(feature=False)
+
+args = parser.parse_args()
+
+
+def plot_blocks(block, fig, **kwargs):
+
+    if 'figsize' in kwargs:
+        fig.set_size_inches(*kwargs['figsize'])
+    else:
+        fig.set_size_inches(12, 12)
+
+    ax = fig.add_subplot()
+    ax.axis('off')
+
+    ax.set_xlim(-0.001, 1.001)
+    ax.set_ylim(-0.001, 1.001)
+
+    rect = patches.Rectangle((0, 0), 1, 1, linewidth=1, edgecolor='k',
+                             facecolor='none')
+    ax.add_patch(rect)
+
+    def get_ind_offsets(inds):
+        return np.concatenate([[0], np.cumsum([len(I) for I in inds])])
+
+    def add_rects(block, c0=(0, 0), w0=1, h0=1):
+        row_offsets = get_ind_offsets(block._row_block_inds)
+        col_offsets = get_ind_offsets(block._col_block_inds)
+        for (i, row), (j, col) in it.product(
+            enumerate(row_offsets[:-1]),
+            enumerate(col_offsets[:-1])
+        ):
+            w = w0*(row_offsets[i + 1] - row)/row_offsets[-1]
+            h = h0*(col_offsets[j + 1] - col)/col_offsets[-1]
+            i0, j0 = c0
+            c = (i0 + w0*row/row_offsets[-1], j0 + h0*col/col_offsets[-1])
+
+            child = block._blocks[i, j]
+            if child.is_leaf:
+                if isinstance(child, cff.FormFactorSvdBlock) or isinstance(child, cff.FormFactorSparseSvdBlock) or \
+                isinstance(child, cff.FormFactorAcaBlock) or isinstance(child, cff.FormFactorSparseAcaBlock) or \
+                isinstance(child, cff.FormFactorIdBlock) or isinstance(child, cff.FormFactorSparseIdBlock):
+                    facecolor = 'cyan' if child.compressed else 'orange'
+                    rect = patches.Rectangle(
+                        c, w, h, edgecolor='none', facecolor=facecolor)
+                    ax.add_patch(rect)
+                elif isinstance(child, cff.FormFactorZeroBlock):
+                    rect = patches.Rectangle(
+                        c, w, h, edgecolor='none', facecolor='black')
+                    ax.add_patch(rect)
+                elif isinstance(child, cff.FormFactorSparseBlock):
+                    rect = patches.Rectangle(
+                        c, w, h, edgecolor='none', facecolor='white')
+                    ax.add_patch(rect)
+                elif isinstance(child, cff.FormFactorDenseBlock):
+                    rect = patches.Rectangle(
+                        c, w, h, edgecolor='none', facecolor='magenta')
+                    ax.add_patch(rect)
+                elif isinstance(child, cff.FormFactorNullBlock):
+                    continue
+                else:
+                    raise Exception('TODO: add %s to _plot_block' % type(child))
+
+                if args.plot_labels:
+                    ax.text(c[0] + 0.5*w, c[1] + 0.5*h, "{:.0e}".format(np.prod(child.shape)), transform=ax.transAxes, fontsize=8, verticalalignment='center', horizontalalignment='center')
+            else:
+                add_rects(child, c, w, h)
+
+            rect = patches.Rectangle(
+                c, w, h, linewidth=1, edgecolor='k', facecolor='none')
+            ax.add_patch(rect)
+
+    add_rects(block)
+
+    ax.invert_xaxis()
+
+    return fig, ax
+
+
+
+
+compression_type = args.compression_type
+max_area_str = str(args.max_area)
+outer_radius_str = str(args.outer_radius)
+sigma_str = str(args.sigma)
+
+max_depth = args.max_depth if args.max_depth != 0 else None
+
+if compression_type == "svd" or compression_type == "aca":
+    if args.add_residuals:
+        FF_dir = "{}_resid_{}_{}_{}_{:.0e}_{}k0".format(compression_type, max_area_str, outer_radius_str, sigma_str, args.tol,
+        args.k0)
+    else:
+        FF_dir = "{}_{}_{}_{}_{:.0e}_{}k0".format(compression_type, max_area_str, outer_radius_str, sigma_str, args.tol,
+            args.k0)
+
+elif compression_type == "rand_svd" or compression_type == "rand_id":
+    if args.add_residuals:
+        FF_dir = "{}_resid_{}_{}_{}_{:.0e}_{}p_{}q_{}k0".format(compression_type, max_area_str, outer_radius_str, sigma_str, args.tol,
+            args.p, args.q, args.k0)
+    else:
+        FF_dir = "{}_{}_{}_{}_{:.0e}_{}p_{}q_{}k0".format(compression_type, max_area_str, outer_radius_str, sigma_str, args.tol,
+            args.p, args.q, args.k0)
+
+
+if not args.cliques and args.min_depth != 1:
+    FF_dir += "_{}mindepth".format(args.min_depth)
+
+if max_depth is not None:
+    FF_dir += "_{}maxdepth".format(max_depth)
+
+if args.compress_sparse:
+    FF_dir += "_cs"
+
+result_dir = "results"
+if args.cliques:
+    FF_dir += "_{}nc".format(args.n_cliques)
+    result_dir += "_cliques"
+
+if args.cliques and args.obb:
+    FF_dir = FF_dir + "_obb"
+
+if args.add_residuals:
+    FF_path =  result_dir+"/"+ FF_dir + "/FF_{}_{}_{:.0e}_{}_resid.bin".format(max_area_str, outer_radius_str, args.tol, compression_type)
+else:
+    FF_path =  result_dir+"/"+ FF_dir + "/FF_{}_{}_{:.0e}_{}.bin".format(max_area_str, outer_radius_str, args.tol, compression_type)
+
+if not os.path.exists(FF_path):
+    print("PATH DOES NOT EXIST " + FF_path)
+    assert False
+FF = CompressedFormFactorMatrix.from_file(FF_path)
+
+
+fig = plt.figure(figsize=(18, 6))  # , figsize=(18, 6))
+print(f'- {FF_path}')
+fig, ax = plot_blocks(FF._root, fig)
+if args.plot_labels:
+    fig.savefig(result_dir+"/"+FF_dir+"/block_plot_labeled.png")
+else:
+    fig.savefig(result_dir+"/"+FF_dir+"/block_plot.png")
+plt.close(fig)