diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/potentials/Tables')
21 files changed, 1514 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries new file mode 100644 index 00000000..f550f801 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Entries @@ -0,0 +1,19 @@ +/divide_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/divide_by_table.m/1.1.1.1/Thu Aug 5 15:25:54 2004// +/extend_domain_table.m/1.1.1.1/Wed Aug 4 15:53:26 2004// +/marg_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/marg_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/marg_table.m/1.1.1.1/Wed Aug 4 15:51:48 2004// +/marg_tableC.c/1.1.1.1/Wed Oct 2 15:39:02 2002// +/marg_tableM.m/1.1.1.1/Tue Oct 1 17:39:08 2002// +/mult_by_sparse_table.c/1.1.1.1/Wed May 29 15:59:58 2002// +/mult_by_table.c/1.1.1.1/Tue Oct 1 21:23:22 2002// +/mult_by_table.m/1.1.1.1/Wed Aug 4 15:52:58 2004// +/mult_by_table2.m/1.1.1.1/Wed Oct 2 15:30:32 2002// +/mult_by_tableC.c/1.1.1.1/Tue Oct 1 21:33:50 2002// +/mult_by_tableM.m/1.1.1.1/Wed Oct 2 15:28:48 2002// +/mult_by_table_global.m/1.1.1.1/Sun Sep 29 10:21:30 2002// +/rep_mult.c/1.1.1.1/Wed May 29 15:59:58 2002// +/repmat_and_mult.c/1.1.1.1/Tue Oct 1 21:20:00 2002// +D diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository new file mode 100644 index 00000000..e52428a3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/potentials/Tables diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c new file mode 100644 index 00000000..d7645ab1 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_sparse_table.c @@ -0,0 +1,126 @@ +/* divide_by_sparse_table.c ../potential/tables*/ + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex; + int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + + plhs[0] = mxDuplicateArray(prhs[0]); + pbDomain = mxGetPr(prhs[1]); + bdim = mxGetNumberOfElements(prhs[1]); + psDomain = mxGetPr(prhs[4]); + sdim = mxGetNumberOfElements(prhs[4]); + + pbSize = mxGetPr(prhs[2]); + psSize = mxGetPr(prhs[5]); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + bpr = mxGetPr(plhs[0]); + bir = mxGetIr(plhs[0]); + bjc = mxGetJc(plhs[0]); + NZB = bjc[1]; + + spr = mxGetPr(prhs[3]); + sir = mxGetIr(prhs[3]); + sjc = mxGetJc(prhs[3]); + NZS = sjc[1]; + + if(sdim == 0){ + value = *spr; + if(value == 0)value = 1; + for(i=0; i<NZB; i++){ + bpr[i] /= value; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + bpr[i] /= spr[position]; + } + } + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c new file mode 100644 index 00000000..6942eddd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.c @@ -0,0 +1,120 @@ +/* divide_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs, value; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + value = *sp; + if(value == 0) value = 1; + for(i=0; i<NB; i++){ + zp[i] /= value; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + value = sp[i]; + if(value == 0) value = 1; + zp[i] /= value; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + value = *sp; + if(value == 0) value = 1; + *zp++ /= value; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m new file mode 100644 index 00000000..6f10718b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/divide_by_table.m @@ -0,0 +1,12 @@ +function bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% DIVIDE_BY_TABLE +% bigT = divide_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +% Replace 0s by 1s before dividing. This is valid, Ts(i)=0 iff Tbig(i)=0. +Ts = Ts + (Ts==0); +%Tbig.T(:) = Tbig.T(:) ./ Ts(:); +bigT(:) = bigT(:) ./ Ts(:); + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m new file mode 100644 index 00000000..4ed0b2cd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/extend_domain_table.m @@ -0,0 +1,24 @@ +function B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz) +% EXTEND_DOMAIN_TABLE Expand an array so it has the desired size. +% B = extend_domain_table(A, smalldom, smallsz, bigdom, bigsz) +% +% A is the array with domain smalldom and sizes smallsz. +% bigdom is the desired domain, with sizes bigsz. +% +% Example: +% smalldom = [1 3], smallsz = [2 4], bigdom = [1 2 3 4], bigsz = [2 1 4 5], +% so B(i,j,k,l) = A(i,k) for i in 1:2, j in 1:1, k in 1:4, l in 1:5 + +if isequal(size(A), [1 1]) % a scalar + B = A; % * myones(bigsz); + return; +end + +map = find_equiv_posns(smalldom, bigdom); +sz = ones(1, length(bigdom)); +sz(map) = smallsz; +B = myreshape(A, sz); % add dimensions for the stuff not in A +sz = bigsz; +sz(map) = 1; % don't replicate along A's dimensions +B = myrepmat(B, sz(:)'); + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c new file mode 100644 index 00000000..f6f7aa26 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_sparse_table.c @@ -0,0 +1,181 @@ +/* marg_sparse_table.c ../potential/tables*/ + +/******************************************/ +/* 5 input & 1 output */ +/* Big sparse table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small sparse table */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr; + const char *field_names[] = {"domain", "T", "sizes"}; + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + bdim = mxGetNumberOfElements(prhs[1]); + sdim = mxGetNumberOfElements(prhs[3]); + pbSize = mxGetPr(prhs[2]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[3]); + bpr = mxGetPr(prhs[0]); + bir = mxGetIr(prhs[0]); + bjc = mxGetJc(prhs[0]); + NZB = bjc[1]; + + if(sdim == 0){ + plhs[0] = mxCreateSparse(1, 1, 1, mxREAL); + spr = mxGetPr(plhs[0]); + bir = mxGetIr(plhs[0]); + bjc = mxGetJc(plhs[0]); + *spr = 0; + *bir = 0; + bjc[0] = 0; + bjc[1] = 1; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + return; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + ssize = malloc(sdim * sizeof(int)); + + NS = 1; + for(i=0; i<count; i++){ + ssize[i] = (int)pbSize[mask[i]]; + NS *= ssize[i]; + } + + for(i=0; i<NZB; i++)sTable[i] = 0; + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * ssize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); + free(ssize); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c new file mode 100644 index 00000000..16329dee --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.c @@ -0,0 +1,175 @@ +/* marg_table.c ../potential/tables */ + + +/******************************************/ +/* 5 input & 1 output */ +/* Big table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small table */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize; + double *pb, *ps, *bp, *sp, *pbd; + + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[3]); + pb = mxGetPr(prhs[1]); + ps = mxGetPr(prhs[3]); + + NB = mxGetNumberOfElements(prhs[0]); + bp = mxGetPr(prhs[0]); + + pbd = mxGetPr(prhs[2]); + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + if(siz_s == 0){ + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + sp = mxGetPr(plhs[0]); + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + ssize = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(ps[i] == pb[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbd[i]; + sy[i] = 1; + } + for(i=0; i<siz_s; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + ssize[i] = sx[temp]; + } + + NS = 1; + for(i=0; i<ndim; i++){ + NS *= sy[i]; + } + + plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL); + sp = mxGetPr(plhs[0]); + + if(NS == 1){ + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++) *sp++ = *bp++; + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + if(maximize){ + for(j=0; j<NB; j++){ + *sp = (*sp < *bp)? *bp : *sp; + bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + else{ + for(j=0; j<NB; j++){ + *sp += *bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); + free(ssize); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m new file mode 100644 index 00000000..0f1e5be5 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_table.m @@ -0,0 +1,27 @@ +function smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) +% MARG_TABLE Marginalize a table +% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) + +if nargin < 5, maximize = 0; end + + +smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array +sum_over = mysetdiff(bigdom, onto); +ndx = find_equiv_posns(sum_over, bigdom); +if maximize + for i=1:length(ndx) + smallT = max(smallT, [], ndx(i)); + end +else + for i=1:length(ndx) + smallT = sum(smallT, ndx(i)); + end +end + + +ns = zeros(1, max(bigdom)); +%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1 +ns(bigdom) = bigsz; + +smallT = squeeze(smallT); % remove all dimensions of size 1 +smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1 diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c new file mode 100644 index 00000000..16329dee --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableC.c @@ -0,0 +1,175 @@ +/* marg_table.c ../potential/tables */ + + +/******************************************/ +/* 5 input & 1 output */ +/* Big table */ +/* Big domain */ +/* Big sizes */ +/* onto */ +/* maximize, if missed, maximize=0 */ +/* */ +/* small table */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize; + double *pb, *ps, *bp, *sp, *pbd; + + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[3]); + pb = mxGetPr(prhs[1]); + ps = mxGetPr(prhs[3]); + + NB = mxGetNumberOfElements(prhs[0]); + bp = mxGetPr(prhs[0]); + + pbd = mxGetPr(prhs[2]); + + if(nrhs < 5) maximize = 0; + else maximize = (int)mxGetScalar(prhs[4]); + + if(siz_s == 0){ + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + sp = mxGetPr(plhs[0]); + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + ssize = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(ps[i] == pb[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbd[i]; + sy[i] = 1; + } + for(i=0; i<siz_s; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + ssize[i] = sx[temp]; + } + + NS = 1; + for(i=0; i<ndim; i++){ + NS *= sy[i]; + } + + plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL); + sp = mxGetPr(plhs[0]); + + if(NS == 1){ + if(maximize){ + for(i=0; i<NB; i++){ + *sp = (*sp < bp[i])? bp[i] : *sp; + } + } + else{ + for(i=0; i<NB; i++){ + *sp += bp[i]; + } + } + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++) *sp++ = *bp++; + free(mask); + free(sx); + free(sy); + free(ssize); + return; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + if(maximize){ + for(j=0; j<NB; j++){ + *sp = (*sp < *bp)? *bp : *sp; + bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + else{ + for(j=0; j<NB; j++){ + *sp += *bp++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + } + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); + free(ssize); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m new file mode 100644 index 00000000..22487dff --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/marg_tableM.m @@ -0,0 +1,30 @@ +function smallT = marg_tableM(bigT, bigdom, bigsz, onto, maximize) +% MARG_TABLE Marginalize a table +% smallT = marg_table(bigT, bigdom, bigsz, onto, maximize) + + % marg_tableM is the same as marg_table.m, but we are sure + % it is not marg_table.c + +if nargin < 5, maximize = 0; end + + +smallT = myreshape(bigT, bigsz); % make sure it is a multi-dim array +sum_over = mysetdiff(bigdom, onto); +ndx = find_equiv_posns(sum_over, bigdom); +if maximize + for i=1:length(ndx) + smallT = max(smallT, [], ndx(i)); + end +else + for i=1:length(ndx) + smallT = sum(smallT, ndx(i)); + end +end + + +ns = zeros(1, max(bigdom)); +%ns(bigdom) = mysize(bigT); % ignores trailing dimensions of size 1 +ns(bigdom) = bigsz; + +smallT = squeeze(smallT); % remove all dimensions of size 1 +smallT = myreshape(smallT, ns(onto)); % put back relevant dims of size 1 diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c new file mode 100644 index 00000000..c1331508 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_sparse_table.c @@ -0,0 +1,155 @@ +/* mult_by_sparse_table.c ../potential/tables*/ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; + + pbDomain = mxGetPr(prhs[1]); + bdim = mxGetNumberOfElements(prhs[1]); + psDomain = mxGetPr(prhs[4]); + sdim = mxGetNumberOfElements(prhs[4]); + + pbSize = mxGetPr(prhs[2]); + psSize = mxGetPr(prhs[5]); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + bpr = mxGetPr(prhs[0]); + bir = mxGetIr(prhs[0]); + bjc = mxGetJc(prhs[0]); + NZB = bjc[1]; + + spr = mxGetPr(prhs[3]); + sir = mxGetIr(prhs[3]); + sjc = mxGetJc(prhs[3]); + NZS = sjc[1]; + + plhs[0] = mxDuplicateArray(prhs[0]); + rpr = mxGetPr(plhs[0]); + rir = mxGetIr(plhs[0]); + rjc = mxGetJc(plhs[0]); + rjc[0] = 0; + rjc[1] = NZB; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + rpr[i] *= *spr; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + rpr[nzCounts] = bpr[i] * spr[position]; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + reset_nzmax(plhs[0], NZB, nzCounts); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c new file mode 100644 index 00000000..3f67983d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.c @@ -0,0 +1,114 @@ +/* mult_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + zp[i] *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + zp[i] *= sp[i]; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m new file mode 100644 index 00000000..a92e340f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table.m @@ -0,0 +1,7 @@ +function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m new file mode 100644 index 00000000..92fb4ac7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table2.m @@ -0,0 +1,30 @@ +function bigT = mult_by_table2(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +%Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +%bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape + +% extend_domain_table has a lot of overhead for small tables, +% since it calls myreshape and myrepmat, which check for 1 dimensional case. +% Here, we check up front. + +if length(bigdom)==1 % vector + bigT = bigT .* smallT; % smallT can be scalar or vector +else + if (length(bigsz) == length(smallsz)) & all(bigsz == smallsz) + bigT = bigT .* smallT; + else + map = find_equiv_posns(smalldom, bigdom); + sz = ones(1, length(bigdom)); + sz(map) = smallsz; + smallT = reshape(smallT, sz); % add dimensions of size 1 for missing domain + % we can use reshape instead of myreshape, because we know length(sz)>1 + sz = bigsz; + sz(map) = 1; % don't replicate along small domain, which is shared + % we can use repmat instead of myrepmat, because we know length(sz)>1 + smallT = repmat(smallT, sz(:)'); + bigT(:) = bigT(:) .* smallT(:); + end +end diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c new file mode 100644 index 00000000..3f67983d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableC.c @@ -0,0 +1,114 @@ +/* mult_by_table.c ../potential/tables */ + + +/******************************************/ +/* 6 input & 1 output */ +/* Big table [0] */ +/* Big domain [1] */ +/* big sizes [2] */ +/* Small table [3] */ +/* small domain [4] */ +/* small sizes [5] */ +/* */ +/* New big table[0] */ +/******************************************/ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, temp; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2; + double *pbDomain, *psDomain, *sp, *zp, *bs; + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + + siz_b = mxGetNumberOfElements(prhs[1]); + siz_s = mxGetNumberOfElements(prhs[4]); + pbDomain = mxGetPr(prhs[1]); + psDomain = mxGetPr(prhs[4]); + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[3]); + sp = mxGetPr(prhs[3]); + + bs = mxGetPr(prhs[2]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + zp[i] *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + zp[i] *= sp[i]; + } + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)bs[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + temp = mask[i]; + sy[temp] = sx[temp]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m new file mode 100644 index 00000000..a92e340f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_tableM.m @@ -0,0 +1,7 @@ +function bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% MULT_BY_TABLE +% bigT = mult_by_table(bigT, bigdom, bigsz, smallT, smalldom, smallsz) +% + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +bigT(:) = bigT(:) .* Ts(:); % must have bigT(:) on LHS to preserve shape diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m new file mode 100644 index 00000000..cc6d7124 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/mult_by_table_global.m @@ -0,0 +1,7 @@ +function mult_by_table_global(bigT, bigdom, bigsz, smallT, smalldom, smallsz) + +% all arguments are read only +global NEWBIGT_GLOBAL + +Ts = extend_domain_table(smallT, smalldom, smallsz, bigdom, bigsz); +NEWBIGT_GLOBAL = bigT(:) .* Ts(:); diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c new file mode 100644 index 00000000..c7b687a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/rep_mult.c @@ -0,0 +1,92 @@ +/* rep_mult.c repmat first two operands to the size provided by */ +/* the third operand, then perform point multiply */ +/* 3 input, 1 output */ +/* C = rep_mult(A, B, sizes) */ + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + double *xp, *yp, *zp, *pSizes; + int xnd, ynd, numElements = 1; + const int *xdim, *ydim; + int i, j, ndim; + int *s, *sx, *sy, *cpsx, *cpsy; + int *subs, *s1, *cpsx2, *cpsy2; + + if (nrhs != 3) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + xnd = mxGetNumberOfDimensions(prhs[0]); + ynd = mxGetNumberOfDimensions(prhs[1]); + xdim = mxGetDimensions(prhs[0]); + ydim = mxGetDimensions(prhs[1]); + ndim = mxGetNumberOfElements(prhs[2]); + + pSizes = mxGetPr(prhs[2]); + + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + s = (int *)malloc(sizeof(int)*ndim); + s1 = (int *)malloc(sizeof(int)*ndim); + *(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1; + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsx2 = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + subs[i] = 0; + sx[i] = (i < xnd) ? xdim[i] : 1; + sy[i] = (i < ynd) ? ydim[i] : 1; + s[i] = (int)pSizes[i]; + s1[i] = s[i] - 1; + numElements *= s[i]; + } + + for(i=0; i<ndim-1; i++){ + cpsx[i+1] = cpsx[i]*sx[i]--; + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsx2[i] = cpsx[i]*sx[i]; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]); + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL); + zp = mxGetPr(plhs[0]); + xp = mxGetPr(prhs[0]); + yp = mxGetPr(prhs[1]); + + for(j=0; j<numElements; j++){ + *zp++ = *xp * *yp; + for(i=0; i<ndim; i++){ + if(subs[i] == s1[i]){ + subs[i] = 0; + if(sx[i]) + xp -= cpsx2[i]; + if(sy[i]) + yp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sx[i]) + xp += cpsx[i]; + if(sy[i]) + yp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(s1); + free(cpsx); + free(cpsy); + free(subs); + free(cpsx2); + free(cpsy2); +} diff --git a/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c new file mode 100644 index 00000000..a7bdba35 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/potentials/Tables/repmat_and_mult.c @@ -0,0 +1,97 @@ +/**************************************************** +A = mult_by_array(big, small) +implicitely copies small |big|/|small| times +and then does element-wise multiplication. + +i.e., +C = repmat(small(:), 1, length(big(:))/length(small(:))) +A = reshape(big(:) .* C(:), size(big)) + +However, this C version avoids the expense of the repmat. + +Written by wei.hu@intel.com, 28 Jan 2002. +/****************************************************/ + + +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + double *sp, *zp; + int i, j, NB, NS, xnd, ynd, ndim; + const int *xdim, *ydim; + int *s, *sx, *sy, *cpsy, *subs, *cpsy2; + + if (nrhs != 2) + mexErrMsgTxt("Incorrect number of inputs."); + + if (nlhs > 1) + mexErrMsgTxt("Too many output arguments."); + + plhs[0] = mxDuplicateArray(prhs[0]); + zp = mxGetPr(plhs[0]); + sp = mxGetPr(prhs[1]); + + xnd = mxGetNumberOfDimensions(prhs[0]); + ynd = mxGetNumberOfDimensions(prhs[1]); + xdim = mxGetDimensions(prhs[0]); + ydim = mxGetDimensions(prhs[1]); + ndim = xnd; + + NB = mxGetNumberOfElements(prhs[0]); + NS = mxGetNumberOfElements(prhs[1]); + + if(NS == 1){ + for(i=0; i<NB; i++){ + *zp++ *= *sp; + } + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + *zp++ *= *sp++; + } + return; + } + + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + subs[i] = 0; + sx[i] = xdim[i]; + sy[i] = (i < ynd) ? ydim[i] : 1; + s[i] = sx[i] - 1; + } + + for (i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + *zp++ *= *sp; + for(i=0; i<ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) sp -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) sp += cpsy[i]; + break; + } + } + } + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); +} |