Commit 87851431 authored by Carl Poelking's avatar Carl Poelking
Browse files

Embedded soap.linalg.permanent and soap.soapy.momo modules.

parent 26d685d3
......@@ -3,6 +3,7 @@ file(GLOB LOCAL_SOURCES *.cpp)
# SUBDIRECTORIES
add_subdirectory(base)
add_subdirectory(linalg)
add_subdirectory(linalg/permanent)
add_subdirectory(tools)
add_subdirectory(soapy)
......
from _linalg import *
from permanent import *
file(GLOB LOCAL_SOURCES *.cpp)
# SUMMARIZE SOURCES
message("Sources: soap/linalg/permanent")
foreach(item ${LOCAL_SOURCES})
message(STATUS " o ${item}")
endforeach()
# ADD & LINK LIBRARY
add_library(_permanent ${LOCAL_SOURCES})
target_compile_options(_permanent PRIVATE "-O3" PRIVATE "-Ofast" PRIVATE "-march=native")
target_link_libraries(_permanent ${PYTHON_LIBRARIES})
set_target_properties(_permanent PROPERTIES PREFIX "" SUFFIX ".so" LIBRARY_OUTPUT_DIRECTORY .)
# INSTALL
install(TARGETS _permanent LIBRARY DESTINATION ${LOCAL_INSTALL_DIR}/linalg/permanent)
install(FILES __init__.py DESTINATION ${LOCAL_INSTALL_DIR}/linalg/permanent)
__all__= [ "permanent_mc", "permanent_ryser", "regmatch" ]
from _permanent import *
// (c) Sandip De
// 16th Oct 2015
// Lausanne
// C++ implementation of python module to compute permanent of a matrix by random montecarlo
/* Functions to compute the permanent, given a numpy array */
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>
#include <vector>
#include <stdlib.h>
#include <algorithm>
#include <random>
#include <cmath>
#include <iostream>
// Array access macros.
#define SM(x0, x1) (*(npy_double*) (( (char*) PyArray_DATA(matrix) + \
(x0) * PyArray_STRIDES(matrix)[0] + \
(x1) * PyArray_STRIDES(matrix)[1])))
#define SM_shape(x0) (int) PyArray_DIM(matrix, x0)
template <typename T>
class Matrix
{
std::vector<T> inner_;
unsigned int dimx_, dimy_;
public:
unsigned int size() const { return dimx_; }
Matrix (unsigned int dimx, unsigned int dimy)
: dimx_ (dimx), dimy_ (dimy)
{
inner_.resize (dimx_*dimy_);
}
inline T operator()(unsigned int x, unsigned int y) const
{
if (x >= dimx_ || y>= dimy_)
throw 0; // ouch
return inner_[dimx_*y + x];
}
inline T& operator()(unsigned int x, unsigned int y)
{
if (x >= dimx_ || y>= dimy_)
throw 0; // ouch
return inner_[dimx_*y + x];
}
};
// Forward function declaration
static PyObject *permanent_mc(PyObject *self, PyObject *args);
// Forward function declaration
static PyObject *permanent_ryser(PyObject *self, PyObject *args);
// Forward function declaration
static PyObject *regmatch(PyObject *self, PyObject *args);
// Method list
static PyMethodDef methods[] = {
{ "permanent_mc", permanent_mc, METH_VARARGS, "Computes the permanent of a numpy matrix by random montecarlo method upto given accuracy"},
{ "permanent_ryser", permanent_ryser, METH_VARARGS, "Computes the permanent of a numpy matrix by Ryser algorithm"},
{ "regmatch", regmatch, METH_VARARGS, "Computes the permanent of a numpy matrix by Ryser algorithm"},
{ NULL, NULL, 0, NULL } // Sentinel
};
// Module initialization
PyMODINIT_FUNC init_permanent(void) {
(void) Py_InitModule("_permanent", methods);
import_array();
}
double fact(int n)
{
double fn=1.0; for (int i=2; i<=n; ++i) fn*=double(i);
return fn;
}
static npy_double _mcperm(PyArrayObject *matrix, PyFloatObject *eps, PyIntObject *ntry, PyIntObject *seed)
{
// int n=mtx.size();
int n = (int) PyArray_DIM(matrix, 0);
std::vector<int> idx(n);
// double eps=1e-3;
double eps1=PyFloat_AS_DOUBLE(eps);
int ntry1=PyInt_AS_LONG(ntry);
int seed1=PyInt_AS_LONG(seed);
for (int i=0; i<n; ++i) idx[i]=i;
double pi, prm=0, prm2=0, fn=fact(n), ti=0;
int i=0, istride=0, pstride=n*100;
if (seed1>0) std::srand(seed1);
//std::cerr<<eps<<" "<<ntry1<<" "<<seed1<<"\n";
while (true)
{
// combines shuffles and cyclic permutations (which are way cheaper!)
if (i%n==0) std::random_shuffle(idx.begin(),idx.end());
else { for (int i=0; i<n; ++i) idx[i]=(idx[i]+1)%n; }
//if (i%10000==0) { for (int j=0; j<n; ++j) std::cerr<<idx[j]<<" "; std::cerr<<"\n"; }
// computes the product of elements for the selected permutation
pi = SM(0, idx[0]);;
for (int j=1; j<n; ++j)
pi *= SM(j, idx[j]);
// accumulates mean and mean square
prm += pi;
prm2 += pi*pi;
++i;
if (ntry1>0 && i >=ntry1) { ti=i; break; }
if (ntry1==0 && i==pstride) // check if we are converged
{
++istride; i=0; ti=double(istride)*double(pstride);
double err=sqrt((prm2-prm*prm/ti)/ti/(ti-1) ) / (prm/ti);
//std::cerr <<istride<< " "<<fn*prm/ti<< " "<<err<< "\n";
if (err< eps1) break;
}
}
//std::cerr <<i<< " "<<fn*prm/ti<< "\n";
return prm/ti*fn;
}
// sinkhorn regularized best match
// NB this assumes that the input matrix is a kernel matrix with entries \in [0,1],
// NB this also works on rectangular matrices
static npy_double _shmatch(PyArrayObject* matrix, PyFloatObject *gamma, PyFloatObject *eps)
{
int nx = (int) PyArray_DIM(matrix, 0);
int ny = (int) PyArray_DIM(matrix, 1);
std::vector<double> u(nx), ou(nx), v(ny);
double ax = 1.0/nx, ay=1.0/ny;
Matrix<double> Kg(nx,ny);
for (int i=0; i<nx; ++i) u[i]=1.0;
for (int i=0; i<ny; ++i) v[i]=1.0;
double lambda=1.0/PyFloat_AS_DOUBLE(gamma), terr=PyFloat_AS_DOUBLE(eps)*PyFloat_AS_DOUBLE(eps), derr;
for (int i=0; i<nx; ++i) for (int j=0; j<ny; ++j) Kg(i,j)=std::exp(-(1-SM(i,j))*lambda);
do
{
// u<-1.0/Kg.v
for (int i=0; i<nx; ++i) { ou[i]=u[i]; u[i]=0.0; }
for (int i=0; i<nx; ++i) for (int j=0; j<ny; ++j) u[i]+=Kg(i,j)*v[j];
// at this point we can compute how far off unity we are
derr = 0.0;
for (int i=0; i<nx; ++i) derr+=(ax-ou[i]*u[i])*(ax-ou[i]*u[i]);
for (int i=0; i<nx; ++i) u[i]=ax/u[i];
// v<-1.0/Kg.u
for (int i=0; i<ny; ++i) v[i]=0.0;
for (int i=0; i<ny; ++i) for (int j=0; j<nx; ++j) v[i]+=Kg(j,i)*u[j];
for (int i=0; i<ny; ++i) v[i]=ay/v[i];
//std::cerr<<derr<<"\n";
} while (derr>terr);
double rval=0, rrow;
for (int i=0; i<nx; ++i)
{
rrow=0;
for (int j=0; j<ny; ++j) rrow+=Kg(i,j)*SM(i,j)*v[j];
rval+=u[i]*rrow;
}
//std::cerr<<"regmatch "<< rval/n <<"\n";
return rval;
}
// Count the number of set bits in a binary string
inline int countbits(unsigned int n)
{
int q=n;
q = (q & 0x5555555555555555) + ((q & 0xAAAAAAAAAAAAAAAA) >> 1);
q = (q & 0x3333333333333333) + ((q & 0xCCCCCCCCCCCCCCCC) >> 2);
q = (q & 0x0F0F0F0F0F0F0F0F) + ((q & 0xF0F0F0F0F0F0F0F0) >> 4);
q = (q & 0x00FF00FF00FF00FF) + ((q & 0xFF00FF00FF00FF00) >> 8);
q = (q & 0x0000FFFF0000FFFF) + ((q & 0xFFFF0000FFFF0000) >> 16);
q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isn't strictly qecessary.
return q;
}
inline int bitparity (unsigned int n) { return 1 - (countbits(n) & 1)*2; }
// Ryser's algorithm
// Adapted from a complex-argument version by Pete Shadbolt
static npy_double _ryperm(PyArrayObject *matrix) {
unsigned int n = (unsigned int) PyArray_DIM(matrix, 0);
npy_double sumz, prod;
npy_double perm = 0;
unsigned long two2n = 1 << n;
unsigned long i, y, z;
for (i=0; i<two2n; ++i) {
prod = 1.0;
for (y=0; y<n; ++y) {
sumz = 0;
for (z=0; z<n; ++z) {
if ((i & (1 << z)) != 0) { sumz += SM(z, y); }
}
prod*=sumz;
}
perm += prod * bitparity(i);
}
if (n%2 == 1) { perm*=-1; }
return perm;
}
// Computes the permanent using a monte carlo scheme
static PyObject *permanent_mc(PyObject *self, PyObject *args) {
// Parse the input
PyArrayObject *matrix;
PyFloatObject *eps;
PyIntObject *ntry;
PyIntObject *seed;
if (!PyArg_ParseTuple(args, "O!O!O!O!", &PyArray_Type, &matrix, &PyFloat_Type, &eps, &PyInt_Type, &ntry, &PyInt_Type, &seed)) {return NULL;}
// Compute the permanent
npy_double p = _mcperm(matrix,eps,ntry,seed);
return PyFloat_FromDouble(p);
}
// Exact permanent based on Ryser algorithm
static PyObject *permanent_ryser(PyObject *self, PyObject *args) {
// Parse the input
PyArrayObject *matrix;
if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &matrix)) {return NULL;}
// Compute the permanent
npy_double p = _ryperm(matrix);
return PyFloat_FromDouble(p);
}
// Computes regularised best-match usin Sinkhorn algorithm
static PyObject *regmatch(PyObject *self, PyObject *args) {
// Parse the input
PyArrayObject *matrix;
PyFloatObject *eps;
PyFloatObject *gamma;
if (!PyArg_ParseTuple(args, "O!O!O!", &PyArray_Type, &matrix, &PyFloat_Type, &gamma, &PyFloat_Type, &eps)) {return NULL;}
// Compute the permanent
npy_double p = _shmatch(matrix,gamma,eps);
return PyFloat_FromDouble(p);
}
install(FILES __init__.py elements.py util.py npthreads.py kernel.py simspace.py pca.py dimred.py lagraph.py lamatch.py math.py DESTINATION ${LOCAL_INSTALL_DIR}/soapy)
install(FILES __init__.py elements.py util.py npthreads.py kernel.py simspace.py pca.py dimred.py lagraph.py lamatch.py math.py momo.py DESTINATION ${LOCAL_INSTALL_DIR}/soapy)
import numpy as np
import numpy.linalg as la
import soap
import permanent as perm
from .. import linalg as perm
import kernel as kern
import lagraph
......
import os
import sys
import commands
import argparse
import time
import numpy as np
from lxml import etree
boolean_dict = \
{'true' : True, '1' : True, 'yes' : True,
'false' : False, '0' : False, 'no' : False, 'none' : False }
# =============================================================================
# XML WRAPPERS
# =============================================================================
class ExtendableNamespace(argparse.Namespace):
def AddNamespace(self, **kwargs):
for name in kwargs:
att = getattr(self, name, None)
if att is None:
setattr(self, name, kwargs[name])
else:
setattr(self, name, kwargs[name].As(type(att)))
return
def Add(self, name, value):
att = getattr(self, name, None)
if att is None:
setattr(self, name, value)
else:
att.Add(name, value)
return value
def GenerateTreeDict(tree, element, path='', paths_rel_to=None):
if type(element) == etree._Comment: return [], {}
# Update path
if path == '':
if element.tag != paths_rel_to:
path += element.tag
else:
path += '/' + element.tag
# Containers for lower levels
tag_node = {}
nodes = []
# Construct Node
xmlnode = XmlNode(element, path) # tree.getpath(element))
nodes.append(xmlnode)
if len(element) == 0:
tag_node[path] = xmlnode
#else:
# print "len 0", xmlnode.path
# tag_node[path] = xmlnode
# Iterate over children
for child in element:
child_elements, childtag_element = GenerateTreeDict(tree, child, path)
nodes = nodes + child_elements
for key in childtag_element.keys():
if tag_node.has_key(key):
if type(tag_node[key]) != list:
tag_node[key] = [ tag_node[key], childtag_element[key] ]
else:
tag_node[key].append(childtag_element[key])
else:
tag_node[key] = childtag_element[key]
return nodes, tag_node
def NamespaceFromDict(tree_dict):
nspace = ExtendableNamespace()
for key in tree_dict.keys():
#print "NSPACE Path = %s" % key
sections = key.split('/')
values = [ None for s in sections ]
values[-1] = tree_dict[key]
add_to_nspace = nspace
for s,v in zip(sections, values):
if v == None:
#print "NSPACE Check for existing"
if getattr(add_to_nspace, s, None):
#print "NSPACE '%s' already exists" % s
add_to_nspace = getattr(add_to_nspace, s, None)
else:
#print "NSPACE '%s' created" % s
sub_nspace = ExtendableNamespace()
add_to_nspace = add_to_nspace.Add(s, sub_nspace)
else:
#print "NSPACE '%s' value added" % s
add_to_nspace.Add(s, v)
return nspace
class XmlTree(list):
def __init__(self, xmlfile, paths_rel_to=None):
self.xmlfile = xmlfile
self.xtree = etree.parse(xmlfile)
self.xroot = self.xtree.getroot()
self.nodes, self.tag_node = GenerateTreeDict(self.xtree, self.xroot, '', paths_rel_to)
self.xspace = NamespaceFromDict(self.tag_node)
def SelectByTag(self, tag):
selection = [ e for e in self.nodes if e.tag == tag ]
return selection
def __getitem__(self, key):
return self.tag_node[key]
def keys(self):
return self.tag_node.keys()
class XmlNode(object):
def __init__(self, element, path):
self.path = path
self.node = element
self.tag = element.tag
self.value = element.text
self.attributes = element.attrib
def As(self, typ):
if typ == np.array:
sps = self.value.split()
return typ([ float(sp) for sp in sps ])
elif typ == bool:
return boolean_dict.get(self.value.lower())
else:
return typ(self.value)
def AsArray(self, typ, sep=' ', rep='\t\n'):
for r in rep:
self.value = self.value.replace(r, sep)
sp = self.value.split(sep)
return [ typ(s) for s in sp if str(s) != '' ]
def SetNodeValue(self, new_value):
self.value = new_value
if self.node != None:
self.node.firstChild.nodeValue = new_value
return
def __getitem__(self, key):
return self.node.get(key)
# =============================================================================
# COMMAND LINE & XML INPUT INTERFACE
# =============================================================================
class CLIO_HelpFormatter(argparse.HelpFormatter):
def _format_usage(self, usage, action, group, prefix):
return "%s : Command Line Interface\n" % sys.argv[0]
class OptionsInterface(object):
def __init__(self):
# COMMAND-LINE ARGUMENTS
self.is_connected_to_cmd_ln = False
self.cmd_ln_args = None
self.cmd_ln_opts = None
self.cmd_ln_nicknames = [ '-h' ]
self.boolean_translator = boolean_dict
self.subtype = str
# XML OPTIONS FILE
self.is_connected_to_xml = False
self.xmlfile = None
self.tree = None
self.xdict = None
self.xspace = None
# JOINED OPTIONS
self.opts = ExtendableNamespace()
def Connect(self, xmlfile=None):
self.ConnectToCmdLn()
self.ConnectToOptionsFile(xmlfile)
def Parse(self, xkey='options'):
if self.is_connected_to_cmd_ln:
self.ParseCmdLn()
if self.is_connected_to_xml:
self.ParseOptionsFileXml(xkey)
if self.is_connected_to_cmd_ln and not self.is_connected_to_xml:
return self.cmd_ln_opts
elif self.is_connected_to_xml and not self.is_connected_to_cmd_ln:
return self.xspace
else:
return self.cmd_ln_opts, self.xspace
def ParseOptionsFile(self, xmlfile, xkey):
self.xmlfile = xmlfile
self.is_connected_to_xml = True
self.ParseOptionsFileXml(xkey)
return self.xspace
# COMMAND-LINE CONVENIENCE ================================================
def __call__(self):
return self.cmd_ln_opts
def ConnectToCmdLn(self, prog=sys.argv[0], descr=None):
self.cmd_ln_args = argparse.ArgumentParser(prog=sys.argv[0],
formatter_class=lambda prog: CLIO_HelpFormatter(prog,max_help_position=70))
self.is_connected_to_cmd_ln = True
return
def ParseCmdLn(self):
self.cmd_ln_opts = self.cmd_ln_args.parse_args()
def InterpretAsBoolean(self, expr):
try:
return self.boolean_translator.get(expr.lower())
except KeyError:
raise ValueError('CLIO does not know how to convert %s into a boolean.' % expr)
def InterpretAsNumpyArray(self, expr):
print "Interpret", expr
array = [ float(e) for e in expr ]
array = np.array(array)
return array
def InterpretAsList(self, expr):
array = [ self.subtype(e) for e in expr ]
return array
def AddArg(self, name, typ=str, nickname=None,
default=None, destination=None, help=None):
# Sort out <name> (e.g. --time) vs <destination> (e.g., time)
if '--' != name[0:2]:
dest = name
name = '--' + name
else:
dest = name[2:]
# Sort out <default> vs <required>
if default == None: required = True
else: required = False
# Construct default <help> it not given
if help == None: help = str(typ)
# Construct <nickname> if not given
if nickname == None:
nickname = '-'
for char in dest:
nickname += char
if not nickname in self.cmd_ln_nicknames:
break
if nickname in self.cmd_ln_nicknames:
raise ValueError('CLIO could not construct nickname from %s option'\
% name)
self.cmd_ln_nicknames.append(nickname)
# Process type
if typ in [int, float, str]:
nargs = None
elif typ == bool:
typ = self.InterpretAsBoolean
nargs = None
elif typ == np.array:
raise NotImplementedError
typ = float # self.InterpretAsNumpyArray
nargs = 3
elif typ == list:
typ = str
nargs = '*'
elif len(typ) == 2 and typ[0] == list:
typ = typ[1]
nargs = '*'
else:
raise NotImplementedError("CLIO does not know how to generate type '%s'"\
% typ)
self.cmd_ln_args.add_argument(nickname, name,
dest=dest,
action='store',
nargs=nargs,
required=required,
type=typ,
metavar=dest[0:1].upper(),
default=default,
help=help)
return
# COMMAND-LINE CONVENIENCE ================================================
def ConnectToOptionsFile(self, xmlfile):
if xmlfile == None or xmlfile == '': return
self.xmlfile = xmlfile
self.is_connected_to_xml = True
return
def ParseOptionsFileXml(self, xkey='options'):
if self.xmlfile == None: return
self.tree = XmlTree(self.xmlfile, paths_rel_to=xkey)
self.xdict = self.tree.tag_node
self.xspace = self.tree.xspace
return
def __getitem__(self, key):
try:
return self.xspace.__dict__[key]
except KeyError:
return self.cmd_ln_opts.__dict__[key]
except KeyError:
raise AttributeError('No such option registered: \'%s\'' % key)
return None
# =============================================================================
# OS SHELL INTERFACE
# =============================================================================
class ShellInterface(object):
def __init__(self):
# =====================================================================
# PRINTER ATTRIBUTES
self.color_dict = { \
'pp' : '\033[95m',
'mb' : '\033[34m',
'lb' : '\033[1;34m',
'my' : '\033[1;33m',
'mg' : '\033[92m',
'mr' : '\033[91m',
'ww' : '\033[0;1m',
'ok' : '\033[92m',
'xx' : '\033[91m',
'warning' : '\033[93m',
'error' : '\033[95m',