Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/marg_sparse_table.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* marg_sparse_table.c ../potential/tables*/
wolffd@0 2
wolffd@0 3 /******************************************/
wolffd@0 4 /* 5 input & 1 output */
wolffd@0 5 /* Big sparse table */
wolffd@0 6 /* Big domain */
wolffd@0 7 /* Big sizes */
wolffd@0 8 /* onto */
wolffd@0 9 /* maximize, if missed, maximize=0 */
wolffd@0 10 /* */
wolffd@0 11 /* small sparse table */
wolffd@0 12 /******************************************/
wolffd@0 13
wolffd@0 14 #include <math.h>
wolffd@0 15 #include <stdlib.h>
wolffd@0 16 #include "mex.h"
wolffd@0 17
wolffd@0 18 int compare(const void* src1, const void* src2){
wolffd@0 19 int i1 = *(int*)src1 ;
wolffd@0 20 int i2 = *(int*)src2 ;
wolffd@0 21 return i1-i2 ;
wolffd@0 22 }
wolffd@0 23
wolffd@0 24 void ind_subv(int index, const int *cumprod, int n, int *bsubv){
wolffd@0 25 int i;
wolffd@0 26
wolffd@0 27 for (i = n-1; i >= 0; i--) {
wolffd@0 28 bsubv[i] = ((int)floor(index / cumprod[i]));
wolffd@0 29 index = index % cumprod[i];
wolffd@0 30 }
wolffd@0 31 }
wolffd@0 32
wolffd@0 33 int subv_ind(const int n, const int *cumprod, const int *subv){
wolffd@0 34 int i, index=0;
wolffd@0 35
wolffd@0 36 for(i=0; i<n; i++){
wolffd@0 37 index += subv[i] * cumprod[i];
wolffd@0 38 }
wolffd@0 39 return index;
wolffd@0 40 }
wolffd@0 41
wolffd@0 42 mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){
wolffd@0 43 mxArray *spTable;
wolffd@0 44 int i, temp, *irs, *jcs, count=0;
wolffd@0 45 double *sr;
wolffd@0 46
wolffd@0 47 spTable = mxCreateSparse(N, 1, nzCounts, mxREAL);
wolffd@0 48 sr = mxGetPr(spTable);
wolffd@0 49 irs = mxGetIr(spTable);
wolffd@0 50 jcs = mxGetJc(spTable);
wolffd@0 51
wolffd@0 52 jcs[0] = 0;
wolffd@0 53 jcs[1] = nzCounts;
wolffd@0 54
wolffd@0 55 for(i=0; i<nzCounts; i++){
wolffd@0 56 irs[i] = sequence[count];
wolffd@0 57 count++;
wolffd@0 58 temp = sequence[count];
wolffd@0 59 sr[i] = Table[temp];
wolffd@0 60 count++;
wolffd@0 61 }
wolffd@0 62 return spTable;
wolffd@0 63 }
wolffd@0 64
wolffd@0 65 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 66 int i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0;
wolffd@0 67 int *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0 68 double *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr;
wolffd@0 69 const char *field_names[] = {"domain", "T", "sizes"};
wolffd@0 70
wolffd@0 71 if(nrhs < 5) maximize = 0;
wolffd@0 72 else maximize = (int)mxGetScalar(prhs[4]);
wolffd@0 73
wolffd@0 74 bdim = mxGetNumberOfElements(prhs[1]);
wolffd@0 75 sdim = mxGetNumberOfElements(prhs[3]);
wolffd@0 76 pbSize = mxGetPr(prhs[2]);
wolffd@0 77 pbDomain = mxGetPr(prhs[1]);
wolffd@0 78 psDomain = mxGetPr(prhs[3]);
wolffd@0 79 bpr = mxGetPr(prhs[0]);
wolffd@0 80 bir = mxGetIr(prhs[0]);
wolffd@0 81 bjc = mxGetJc(prhs[0]);
wolffd@0 82 NZB = bjc[1];
wolffd@0 83
wolffd@0 84 if(sdim == 0){
wolffd@0 85 plhs[0] = mxCreateSparse(1, 1, 1, mxREAL);
wolffd@0 86 spr = mxGetPr(plhs[0]);
wolffd@0 87 bir = mxGetIr(plhs[0]);
wolffd@0 88 bjc = mxGetJc(plhs[0]);
wolffd@0 89 *spr = 0;
wolffd@0 90 *bir = 0;
wolffd@0 91 bjc[0] = 0;
wolffd@0 92 bjc[1] = 1;
wolffd@0 93 if(maximize){
wolffd@0 94 for(i=0; i<NZB; i++){
wolffd@0 95 *spr = (*spr < bpr[i])? bpr[i] : *spr;
wolffd@0 96 }
wolffd@0 97 }
wolffd@0 98 else{
wolffd@0 99 for(i=0; i<NZB; i++){
wolffd@0 100 *spr += bpr[i];
wolffd@0 101 }
wolffd@0 102 }
wolffd@0 103 return;
wolffd@0 104 }
wolffd@0 105
wolffd@0 106 mask = malloc(sdim * sizeof(int));
wolffd@0 107 count = 0;
wolffd@0 108 for(i=0; i<sdim; i++){
wolffd@0 109 for(j=0; j<bdim; j++){
wolffd@0 110 if(psDomain[i] == pbDomain[j]){
wolffd@0 111 mask[count] = j;
wolffd@0 112 count++;
wolffd@0 113 break;
wolffd@0 114 }
wolffd@0 115 }
wolffd@0 116 }
wolffd@0 117
wolffd@0 118 sTable = malloc(NZB * sizeof(double));
wolffd@0 119 sequence = malloc(NZB * 2 * sizeof(double));
wolffd@0 120 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 121 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 122 bsubv = malloc(bdim * sizeof(int));
wolffd@0 123 ssubv = malloc(sdim * sizeof(int));
wolffd@0 124 ssize = malloc(sdim * sizeof(int));
wolffd@0 125
wolffd@0 126 NS = 1;
wolffd@0 127 for(i=0; i<count; i++){
wolffd@0 128 ssize[i] = (int)pbSize[mask[i]];
wolffd@0 129 NS *= ssize[i];
wolffd@0 130 }
wolffd@0 131
wolffd@0 132 for(i=0; i<NZB; i++)sTable[i] = 0;
wolffd@0 133
wolffd@0 134 bCumprod[0] = 1;
wolffd@0 135 for(i=0; i<bdim-1; i++){
wolffd@0 136 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 137 }
wolffd@0 138 sCumprod[0] = 1;
wolffd@0 139 for(i=0; i<sdim-1; i++){
wolffd@0 140 sCumprod[i+1] = sCumprod[i] * ssize[i];
wolffd@0 141 }
wolffd@0 142
wolffd@0 143 count = 0;
wolffd@0 144 for(i=0; i<NZB; i++){
wolffd@0 145 bindex = bir[i];
wolffd@0 146 ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0 147 for(j=0; j<sdim; j++){
wolffd@0 148 ssubv[j] = bsubv[mask[j]];
wolffd@0 149 }
wolffd@0 150 sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0 151 result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare);
wolffd@0 152 if(result){
wolffd@0 153 position = (result - sequence) / 2;
wolffd@0 154 if(maximize)
wolffd@0 155 sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position];
wolffd@0 156 else sTable[position] += bpr[i];
wolffd@0 157 }
wolffd@0 158 else {
wolffd@0 159 if(maximize)
wolffd@0 160 sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts];
wolffd@0 161 else sTable[nzCounts] += bpr[i];
wolffd@0 162 sequence[count] = sindex;
wolffd@0 163 count++;
wolffd@0 164 sequence[count] = nzCounts;
wolffd@0 165 nzCounts++;
wolffd@0 166 count++;
wolffd@0 167 }
wolffd@0 168 }
wolffd@0 169
wolffd@0 170 qsort(sequence, nzCounts, sizeof(int) * 2, compare);
wolffd@0 171 plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS);
wolffd@0 172
wolffd@0 173 free(sTable);
wolffd@0 174 free(sequence);
wolffd@0 175 free(mask);
wolffd@0 176 free(bCumprod);
wolffd@0 177 free(sCumprod);
wolffd@0 178 free(bsubv);
wolffd@0 179 free(ssubv);
wolffd@0 180 free(ssize);
wolffd@0 181 }