memo

This is a memo to document what I did for my master thesis on a daily
basis. It not only includes the stuff that I successfully implemented,
but also the problems that I encountered. Moreover, it will also contain
the discussion with my advisor and though process.

I hope this file could guide me to a clear goal of our project.


15.12.2020
Modified set_env.sh, CMakeLists.txt and wrote a new swig.sh to add a
linker file with liblapack, libblas libraries in the run time.

Figured out seg fault problem in call_Launchhelper(). We must initialize
a separate memeory block in runtime


16.12.2020

Reading day!

Read format in .i file about swig: %include, %inline, global var...,
static linking, extern


17.12.2020

Designed basic working case: first start with the implementation of
running dense matrix

wrote functions: SPDMATRIX_DENSE, launchhelper_denseSPD

problem encountered: see the issue in gitlab with the title 

        Cant read the binary matrix file from .sh


18.12.2020

Solved the issue from 17.12.2020 (see solution in gitlab)

modified tools.i file so that SPDMATRIX_DENSE function is callable in python


20.12.2020
1. rewrote launchhelper_denseSPD so that it could take as an parameter
a file path to the parameters file. It runs succesfully

2. Modified the constructor of gofmm::CommandLineHelper so that it can
take a vector as a parameter. This proves much cleaner and understandble
since we can avoid realloc (and general memory management) which typically happens when
we want to read file into char**

3. wrote file_to_argv: which reads a file into vector. Due to memory
issue related to local vars, this function is still under construction.

4. wrote `compile_log.txt`: document what the output of key cmds in the
compilitaion should look like. Users who intend to install gofmm + hmlp package
should consult it.


21.12.2020
1. Finished part 3 of the date 20.12.2020. Namely, wrote file_to_argv
class to wrap parameters from a txt file  into a class object.

  a. This class reduces hard coding, eg, when we need to run two different spd
     matrix, testsuit + spd dense, we don't need hard coding
  b. Use destructor to smartly clean up the memory space occupied by the vector

It runs successfully

2. Wrote test_gofmm.h file so that functions in test_gofmm.cpp could
reference each other without worrying about order


22.12.2020
1. Add all current functions in .cpp to its .h file so that
  a. The interface is cleaner
  b. No restriction of function call order in .cpp


05.01.2021
1. Reading day:)
2. Create infrastructure for storing matrix in a tree (see typedef in .h)
3. Finished the function compress
4. Finished the function evaluate


08.01.2021
1. Changed `cmds`

TASK: implement a function as folllows which takes `AA` as an input numpy
matrix declared in the python console

tools.launchhelper_denseSPD(AA, "parameters_dense.txt")


19.01.2021

1. Read where the matrix val was stored by looking at the `vector` and
`Data` class

2. Updated the `tools.i` file to implement the `numpy` usage

3. Wrote the function `load_denseSPD_from_console` to take a `numpy` array
as input and produces a tree-contained spd matrix as output.

Here are the stuff that I consulted for implementation

How to deal with numpy in c++
https://numpy.org/doc/stable/reference/swig.html

Case with 1D numpy array to work with:
https://stackoverflow.com/questions/52199057/c-array-to-numpy-with-swig-problem-with-typeerror

https://numpy.org/doc/stable/reference/swig.interface-file.html

// since the class inherits from vector (public vector<T, Allocator>),
// it has the method data(). Documentation is as follows:
// http://www.cplusplus.com/reference/vector/vector/data/
file.seekg( 0, std::ios::beg );
if ( !file.read( (char*)this->data(), size ) )
{
std::cerr << "ERROR: only " << file.gcount() << " bytes, expecting " << m * n * sizeof(T) << std::endl;
}
file.close();


20.01.2021

1. Further updated `tools.i` so that a `numpy` array could be imported
succesfully in `swig`

2. Wrote `test.py` to demonstrate how the task to implement `numpy` and
`gofmm` together works

3. For user-input lower dimensional `numpy` array, we must fine-tuning the
# of children of the tree.

In the end, the `test.py` runs succesfully for numpy array of size > 10
(both gofmm and accuracy report work) but for size < 10 only gofmm works.

[getianyi@Ubuntu build]$ python3 test.py
testsuit
partitioning ...
[ RT]     1 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
[ RT]     1 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
omp_get_max_threads() 4
Dependency clean up
[ RT]     1 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
Finish NearSamplesTask
Finish SymmetrizeNearInteractions
Skeletonization (HMLP Runtime) ...
[ RT]     4 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
MergeFarNodes ...
CacheFarNodes ...
========================================================
GOFMM compression phase
========================================================
NeighborSearch ------------------------  0.00s (  0.0%)
partitioning ----------------------  0.01s ( 88.2%)
Skeletonization -----------------------  0.00s ( 11.3%)
MergeFarNodes -------------------------  0.00s (  0.1%)
CacheFarNodes -------------------------  0.00s (  0.4%)
========================================================
Compress (1.00 not compressed) --------  0.02s (100.0%)
========================================================

Forward permute ...
N2S, S2S, S2N, L2L (HMLP Runtime) ...
[ RT]     7 [normal]     0 [nested] 9.216E+03 flops 9.243E+03 mops
d2h_t 0.000000s
aggregate_t 0.000000s
Backward permute ...
========================================================
GOFMM evaluation phase
========================================================
Allocate ------------------------------  0.00s (  3.8%)
Forward permute -----------------------  0.00s ( 14.9%)
N2S, S2S, S2N, L2L --------------------  0.00s ( 65.8%)
Backward permute ----------------------  0.00s ( 15.4%)
========================================================
Evaluate ------------------------------  0.00s (100.0%)
========================================================

========================================================
Accuracy report
========================================================
gid      0, ASKIT 0.0E+00, HODLR 0.0E+00, GOFMM 0.0E+00
gid      1, ASKIT 0.0E+00, HODLR 0.0E+00, GOFMM 0.0E+00
gid      2, ASKIT 0.0E+00, HODLR 0.0E+00, GOFMM 0.0E+00
========================================================
Elementwise ASKIT 0.0E+00, HODLR 0.0E+00, GOFMM 0.0E+00
========================================================
[ RT]     1 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
python3: ../gofmm/igofmm.hpp:652: void hmlp::gofmm::Factor<T>::Telescope(bool, hmlp::Data<T>&, hmlp::Data<T>&) [with T = float]: Assertion `Palr.col() == n' failed.
Aborted (core dumped)


Task: find a way to fine tuning the children #


22.01.2021

1. Find how to fine tuning the parameters `m` and `s` in accordance with
`n`


## problem size
n=5000
## maximum leaf node size
m=64
## maximum off-diagonal ranks
s=64
## number of neighbors
k=0
## number of right hand sides
nrhs=512
## user tolerance
stol=1E-5
## user computation budget [0,1]
budget=0.00
## distance type (geometry, kernel, angle)
distance="angle"
## spdmatrix type (testsuit, dense, ooc, kernel, userdefine)
matrixtype="testsuit"
## kernelmatrix type (gaussian, laplace)
kerneltype="gaussian"


25.01.2021

1. Combined all 3 cases: testsuit, dense spd and numpy into the file 
`test.py`

Task to do:

1.   auto u = Evaluate( tree, w );  // type: Data<T>

extract data from the container and do multiplication

question: multiplication 

should 
       Data<T> mul(Data<T> ..., Data<T> ...) {}



29.01.2021

# How to print a c matrix in numpy array form in python
https://stackoverflow.com/questions/36222455/swigcpython-passing-and-receiving-c-arrays

# An example of 2D spd matrix
# spdMatrix = np.array([[2, 6], [6, 20]], dtype=np.double)

1. Included mul_denseSPD in test_gofmm.cpp to calcuate the multiplication
of two dense SPD matrices

2. Wrote typemap for an 1D output matrix which will be resize to display
the instrinsic 2D layout (tools.i)
    Lesson learned:  (double* ARGOUT_ARRAY2, int DIM1, int DIM2) is unavailable
    Whenever we want to use a typemap in .i file, check first if such
    mapping is defined in numpy.i

3. Updated test.py to include the running example of multiplication of
two numpy spd matrices


05.02.2020

1. Inverse of dense SPD

https://www.cs.utexas.edu/ftp/techreports/tr06-20.pdf

https://stackoverflow.com/questions/40703042/more-efficient-way-to-invert-a-matrix-knowing-it-is-symmetric-and-positive-semi


06.02.2020
1. Read about how blas implement inverse of a matrix, solve a linear system,

http://www.netlib.org/lapack/double/


07.02.2021

1. Use the following method to compute the inverse of K


    void Factorize( Data<T> &Kaa ) 
    {
      assert( isLeaf() );
      assert( Kaa.row() == n ); assert( Kaa.col() == n );

      /** Initialize with Kaa. */
      Z = Kaa;

      /** Record the partial pivoting order. */
      ipiv.resize( n, 0 );

      /** Compute 1-norm of Z. */
      T nrm1 = 0.0;
      for ( auto &z : Z ) nrm1 += z;

      /** Pivoted LU factorization. */
      xgetrf( n, n, Z.data(), n, ipiv.data() );

      /** Compute 1-norm condition number. */
      T rcond1 = 0.0;
      Data<T> work( Z.row(), 4 );
      vector<int> iwork( Z.row() );
      xgecon( "1", Z.row(), Z.data(), Z.row(), nrm1, 
          &rcond1, work.data(), iwork.data() );
      if ( 1.0 / rcond1 > 1E+6 )
        printf( "Warning! large 1-norm condition number %3.1E, nrm1( Z ) %3.1E\n", 
            1.0 / rcond1, nrm1 );
    }; /** end Factorize() */


2. Problem: how to convert Tree class to Data class


08.02.2021 - 10.02.2021
1. Read intensively about the class functions, fortran files and blas
prototypes to figure out which functions should be used to compute the
inverse (It's time consuming because the existing code don't have adequate
commentary. Besides, the coding style makes it hard to use egrep to grap
key fields. Eg, public fields spread randomly inside a class. Random
commas here and there -> making it extremly time consuming to use
egrep espcially in such large project)

2. Attempted to write a new Factorize class that takes in the root node
and do LUV or LV computation. But failed because something goes wrong
with the input.


11.02.2021
1. Wrote `invert_denseSPD` to convert a SPD matrix.

2. Updated `test.py` to cover the inverse user case


12.02.2021
1. Debugging with Severin

2. Task to do:
   a. figure out why inverse is instable
   b. Create github


15.02.2021
1. Tried printing out inverse matrix value from all nodes. But the same
problem after the 4th row or somewhere there, entries start to contain
very small vals -> det(inverse) close to 0

2. Task for tomorrow: will not use factorization but pseudo-inverse and
see how that goes


16.02.2021
1. Tried using xgetrf and inv(k) = w * inv(u). The result is a lot closer
than just using the factorize. But still unstable


24.02.2021 -26.02.2021
1. Wrote the error file but lost it during merge :(


02.03.2021
1. Rewrote multiplication function using gofmm

03.03.2021
1. Work on error file again but this time incorporate egrep into python
2. Wrote 3 functions to compute the rses

04.03.2021
1. Incorporate the functions wrote on the date 03.03.2021 into a class,
rse_calculator (benefits of a class :))

2. Notable performance: take a look at the last entry which corresponds
to the leaf node size 1024

[0.8612689165197026, 0.842358990265323, 0.8481395839236431, 0.8183800151701816, 0.751166024988613, 0.7121029139939241, 0.6134476531553748, 7.63604173110931e-05]

e-05: this is a notable error drop compared to all previous entries corresponding
to leaf node size from 8 to 512.

3. Further results

The relative standard error at the leaf node size 8, 32, 128, 512 is
[0.878310270424629, 0.845326314465862, 0.7778199424492407, 0.3471068224167622]

The relative standard error at the tolerance 1E-3, 1E-5, 1E-7 is
[0.34710689297322533, 0.33541029792654214, 0.33688449797138287]


05.03.2021
1. Calculate the rses for the configuration of two sets, stols x leafNodeSize
            8                32              128             512             1024
1E-3 array([[4.41262073e-01, 4.39609106e-01, 4.60542097e-01, 2.68897813e-01, 3.72566878e-04],
1E-5       [4.63088328e-01, 4.55703526e-01, 4.68613841e-01, 2.68898622e-01, 6.48993133e-05],
1E-7       [4.40821804e-01, 4.42426887e-01, 4.59727435e-01, 2.65221226e-01, 7.00785465e-05]])

note: leafNodeSize = 1024 -> notable drop in rse

2. stuff to read
a. smth to do with factorization
https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution
        
   factorize: calculate LUVP... (lower triangular, upper triangular, V,
              permutation,...) for K, the SPD
   solve: use LUV to solve the linear eq system Kw = u where K, u known

b. New thought to calculate the inverse of the SPD:
   Let the potential be n x n identity matrix and use the solver to get
   w. This w is effectively inv(K)


08.03.2021
1. Wrote 2 child classes of rse_calculator:
   a. spdMultiplication_rseCalulator
   b. spdInverse_rseCalulator
to calculate the rses of matrix mul and inv

row: stols = [1E-3, 1E-5, 1E-7]  
col: rank = [8, 32, 128, 256, 512, 1024] (also is leaf node size)

multiplication error

[[8.28216761e-01 8.17690475e-01 7.30507649e-01 4.25563318e-01
  5.33191277e-05 5.06015293e-05]
 [8.21753934e-01 8.31821189e-01 7.50877147e-01 4.25563018e-01
  5.22699376e-05 5.21239007e-05]
 [8.14985462e-01 8.25909278e-01 7.44181603e-01 4.25562888e-01
  5.05020967e-05 5.09992211e-05]]

inverse error

[[5.00332286e+00 5.00652777e+00 4.36099554e+00 2.42939717e+00
  2.27446961e-04 2.27446961e-04]
 [5.00156093e+00 5.00886389e+00 4.35608443e+00 2.42939559e+00
  2.27446961e-04 2.27446961e-04]
 [5.00321689e+00 5.01030233e+00 4.35736420e+00 2.42939697e+00
  2.27446961e-04 2.27446961e-04]]


09.03.2021
1. Tasks: a. Use the files from Python/ from Severin to construct a dense SPD
          b. Run rses on them

2. Finished the 2 tasks as required in the point 1


12.03.2021
1. https://packaging.python.org/tutorials/packaging-projects/
2. https://python-packaging-tutorial.readthedocs.io/en/latest/uploading_pypi.html
tools.so + tools.py + ...
pip install


15.03.2021
1. The following name has the last successful commit before we start working
on replace the configuration file with inputs in local script. If somehow
the replacement doesn't work, git pull the repo with the following commit
name

last Copy before ditching configFile

2. Working on the gofmmTree. Finished its default and user-defined
constructors


16.03.2021
1. Finished both mul + inv

2. modified rse_calculator and its derived classes, spdMultiplication_rseCalulator
and spdInverse_rseCalulator



18.03.2021

1. Shit was going on, was stupid enough to upgrade my Ubunut to 20.04 with a broken
upgrade. Now goes back to the stone age.

2. Strugging with the pt 1 for the most part of the day but somehow managed to
almost finish the pip install gofmm1

3. Read the follwoing pipy upload stuff:

https://medium.com/@joel.barmettler/how-to-upload-your-python-package-to-pypi-65edc5fe9c56


19.03.2021

1. figure out relative import issue now gofmm1 is good to go


26.04.2021



Size = 4096, random generation

-----------------------------------------------------------
RSE of matrix multiplication
-----------------------------------------------------------
row: stol, col: rank
[[ 0.68789079  0.69009854  0.65449571]
 [ 0.69010893  0.68689731  0.6522179 ]
 [ 0.68935427  0.69502891  0.65675001]]


-----------------------------------------------------------
RSE of matrix inverse
-----------------------------------------------------------
row: stol, col: rank
[[ 5.25490516  5.28997791  5.20992468]
 [ 5.25620159  5.2886143   5.20659196]
 [ 5.2557793   5.2896646   5.20770443]]



Size = 512, random generation

-----------------------------------------------------------
RSE of matrix multiplication
-----------------------------------------------------------
row: stol, col: rank
[[ 0.8657503   0.92984818  0.50276503]
 [ 0.84694693  0.85230478  0.50276433]
 [ 0.85074576  0.90364598  0.50276519]]


-----------------------------------------------------------
RSE of matrix inverse
-----------------------------------------------------------
row: stol, col: rank
[[ 4.91141171  4.98151767  2.27484366]
 [ 4.91254531  4.95907969  2.27484344]
 [ 4.909484    4.95820075  2.27484331]]
 

size = 512, kde

-----------------------------------------------------------
RSE of matrix multiplication
-----------------------------------------------------------
row: stol, col: rank
[[  9.78082888e-02   9.98502079e-02   9.14858713e-02]
 [  9.25327829e-02   2.57159204e-03   2.94960048e-03]
 [  1.04849655e-01   5.73086316e-05   3.79355086e-05]]


-----------------------------------------------------------
RSE of matrix inverse
-----------------------------------------------------------
row: stol, col: rank
[[  1.52237376e-03   1.38722317e-03   1.40763783e-03]
 [  1.48353771e-03   4.95824899e-05   4.82194998e-05]
 [  1.44550701e-03   3.59018948e-05   5.84606678e-05]]


size = 512, kde

========================================================
GOFMM compression phase
========================================================
NeighborSearch ------------------------  0.00s (  0.0%)
partitioning ----------------------  0.01s (  4.7%)
Skeletonization -----------------------  0.28s ( 89.4%)
MergeFarNodes -------------------------  0.00s (  0.2%)
CacheFarNodes -------------------------  0.02s (  5.7%)
========================================================
Compress (0.12 not compressed) --------  0.32s (100.0%)
========================================================

[ RT]    15 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
[ RT]    15 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
do ulv[ RT]    53 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
[ RT]     8 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
[ RT]     0 [normal]     0 [nested] 0.000E+00 flops 0.000E+00 mops
swig/python detected a memory leak of type 'DATA_s *', no destructor found.
swig/python detected a memory leak of type 'SPDMATRIX_DENSE *', no destructor found.
swig/python detected a memory leak of type 'SPDMATRIX_DENSE *', no destructor found.


-----------------------------------------------------------
RSE of matrix multiplication
-----------------------------------------------------------
row: stol, col: rank
[[  9.78082888e-02   9.98502079e-02   9.14858713e-02]
 [  9.25327829e-02   2.57159204e-03   2.94960048e-03]
 [  1.04849655e-01   5.73086316e-05   3.79355086e-05]]


-----------------------------------------------------------
RSE of matrix inverse
-----------------------------------------------------------
row: stol, col: rank
[[  1.52237376e-03   1.38722317e-03   1.40763783e-03]
 [  1.48353771e-03   4.95824899e-05   4.82194998e-05]
 [  1.44550701e-03   3.59018948e-05   5.84606678e-05]]


size: 2048, kde

========================================================
GOFMM compression phase
========================================================
NeighborSearch ------------------------  0.00s (  0.0%)
partitioning ----------------------  0.08s (  3.0%)
Skeletonization -----------------------  2.63s ( 92.9%)
MergeFarNodes -------------------------  0.01s (  0.2%)
CacheFarNodes -------------------------  0.11s (  3.9%)
========================================================
Compress (0.03 not compressed) --------  2.83s (100.0%)
========================================================



-----------------------------------------------------------
RSE of matrix multiplication
-----------------------------------------------------------
row: stol, col: rank
[[  1.02383007e-01   6.47416878e-02   6.21142321e-02]
 [  1.10221497e-01   5.71518148e-04   4.91981278e-04]
 [  1.12021029e-01   8.71955749e-05   7.54597947e-05]]


-----------------------------------------------------------
RSE of matrix inverse
-----------------------------------------------------------
row: stol, col: rank
[[  2.70073274e-03   1.33993375e-03   1.52413063e-03]
 [  3.19061439e-03   5.05631021e-05   4.36224789e-05]
 [  2.73955893e-03   3.70981401e-05   6.99600012e-05]]  


real    11m25.654s
user    11m23.796s
sys     0m1.488s


size: 4096, random

-----------------------------------------------------------                                          
RSE of matrix multiplication                                                                         
-----------------------------------------------------------                                          
row: stol, col: rank                                                                                 
[[ 0.26455808  0.28034735]                                                                           
 [ 0.26529322  0.2829262 ]                                                                           
 [ 0.26451882  0.28591586]]                                                                          
                                                                                                     
                                                                                                     
-----------------------------------------------------------                                          
RSE of matrix inverse                                                                                
-----------------------------------------------------------                                          
row: stol, col: rank                                                                                 
[[ 5.27702629  5.19596347]                                                                           
 [ 5.27678235  5.20095586]                                                                           
 [ 5.27707929  5.19531117]]

real    14m2.558s
user    13m58.782s
sys     0m3.283s

11.05.2021
1. Read Gaussian KDE
https://freelunchtheorem.github.io/Conditional_Density_Estimation/docs/html/density_estimator/nkde.html


2. What parameters to inspect: rank + leaf node size

   add D.  Complexity Analysis to the theory
https://mediatum.ub.tum.de/doc/1522878/1522878.pdf

3. bib docu
https://web.mit.edu/rsi/www/pdfs/bibtex-format.pdf