Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/marg_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_table.c ../potential/tables */
wolffd@0 2
wolffd@0 3
wolffd@0 4 /******************************************/
wolffd@0 5 /* 5 input & 1 output */
wolffd@0 6 /* Big table */
wolffd@0 7 /* Big domain */
wolffd@0 8 /* Big sizes */
wolffd@0 9 /* onto */
wolffd@0 10 /* maximize, if missed, maximize=0 */
wolffd@0 11 /* */
wolffd@0 12 /* small table */
wolffd@0 13 /******************************************/
wolffd@0 14
wolffd@0 15 #include "mex.h"
wolffd@0 16
wolffd@0 17 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 18 int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
wolffd@0 19 int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
wolffd@0 20 double *pb, *ps, *bp, *sp, *pbd;
wolffd@0 21
wolffd@0 22
wolffd@0 23 siz_b = mxGetNumberOfElements(prhs[1]);
wolffd@0 24 siz_s = mxGetNumberOfElements(prhs[3]);
wolffd@0 25 pb = mxGetPr(prhs[1]);
wolffd@0 26 ps = mxGetPr(prhs[3]);
wolffd@0 27
wolffd@0 28 NB = mxGetNumberOfElements(prhs[0]);
wolffd@0 29 bp = mxGetPr(prhs[0]);
wolffd@0 30
wolffd@0 31 pbd = mxGetPr(prhs[2]);
wolffd@0 32
wolffd@0 33 if(nrhs < 5) maximize = 0;
wolffd@0 34 else maximize = (int)mxGetScalar(prhs[4]);
wolffd@0 35
wolffd@0 36 if(siz_s == 0){
wolffd@0 37 plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
wolffd@0 38 sp = mxGetPr(plhs[0]);
wolffd@0 39 if(maximize){
wolffd@0 40 for(i=0; i<NB; i++){
wolffd@0 41 *sp = (*sp < bp[i])? bp[i] : *sp;
wolffd@0 42 }
wolffd@0 43 }
wolffd@0 44 else{
wolffd@0 45 for(i=0; i<NB; i++){
wolffd@0 46 *sp += bp[i];
wolffd@0 47 }
wolffd@0 48 }
wolffd@0 49 return;
wolffd@0 50 }
wolffd@0 51
wolffd@0 52 mask = malloc(siz_s * sizeof(int));
wolffd@0 53 ssize = malloc(siz_s * sizeof(int));
wolffd@0 54 count = 0;
wolffd@0 55 for(i=0; i<siz_s; i++){
wolffd@0 56 for(j=0; j<siz_b; j++){
wolffd@0 57 if(ps[i] == pb[j]){
wolffd@0 58 mask[count] = j;
wolffd@0 59 count++;
wolffd@0 60 break;
wolffd@0 61 }
wolffd@0 62 }
wolffd@0 63 }
wolffd@0 64
wolffd@0 65 ndim = siz_b;
wolffd@0 66 sx = (int *)malloc(sizeof(int)*ndim);
wolffd@0 67 sy = (int *)malloc(sizeof(int)*ndim);
wolffd@0 68 for(i=0; i<ndim; i++){
wolffd@0 69 sx[i] = (int)pbd[i];
wolffd@0 70 sy[i] = 1;
wolffd@0 71 }
wolffd@0 72 for(i=0; i<siz_s; i++){
wolffd@0 73 temp = mask[i];
wolffd@0 74 sy[temp] = sx[temp];
wolffd@0 75 ssize[i] = sx[temp];
wolffd@0 76 }
wolffd@0 77
wolffd@0 78 NS = 1;
wolffd@0 79 for(i=0; i<ndim; i++){
wolffd@0 80 NS *= sy[i];
wolffd@0 81 }
wolffd@0 82
wolffd@0 83 plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
wolffd@0 84 sp = mxGetPr(plhs[0]);
wolffd@0 85
wolffd@0 86 if(NS == 1){
wolffd@0 87 if(maximize){
wolffd@0 88 for(i=0; i<NB; i++){
wolffd@0 89 *sp = (*sp < bp[i])? bp[i] : *sp;
wolffd@0 90 }
wolffd@0 91 }
wolffd@0 92 else{
wolffd@0 93 for(i=0; i<NB; i++){
wolffd@0 94 *sp += bp[i];
wolffd@0 95 }
wolffd@0 96 }
wolffd@0 97 free(mask);
wolffd@0 98 free(sx);
wolffd@0 99 free(sy);
wolffd@0 100 free(ssize);
wolffd@0 101 return;
wolffd@0 102 }
wolffd@0 103
wolffd@0 104 if(NS == NB){
wolffd@0 105 for(i=0; i<NB; i++) *sp++ = *bp++;
wolffd@0 106 free(mask);
wolffd@0 107 free(sx);
wolffd@0 108 free(sy);
wolffd@0 109 free(ssize);
wolffd@0 110 return;
wolffd@0 111 }
wolffd@0 112
wolffd@0 113 s = (int *)malloc(sizeof(int)*ndim);
wolffd@0 114 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
wolffd@0 115 subs = (int *)malloc(sizeof(int)*ndim);
wolffd@0 116 cpsy2 = (int *)malloc(sizeof(int)*ndim);
wolffd@0 117 for(i = 0; i < ndim; i++){
wolffd@0 118 subs[i] = 0;
wolffd@0 119 s[i] = sx[i] - 1;
wolffd@0 120 }
wolffd@0 121
wolffd@0 122 for(i = 0; i < ndim-1; i++){
wolffd@0 123 cpsy[i+1] = cpsy[i]*sy[i]--;
wolffd@0 124 cpsy2[i] = cpsy[i]*sy[i];
wolffd@0 125 }
wolffd@0 126 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
wolffd@0 127
wolffd@0 128 if(maximize){
wolffd@0 129 for(j=0; j<NB; j++){
wolffd@0 130 *sp = (*sp < *bp)? *bp : *sp;
wolffd@0 131 bp++;
wolffd@0 132 for(i = 0; i < ndim; i++){
wolffd@0 133 if(subs[i] == s[i]){
wolffd@0 134 subs[i] = 0;
wolffd@0 135 if(sy[i])
wolffd@0 136 sp -= cpsy2[i];
wolffd@0 137 }
wolffd@0 138 else{
wolffd@0 139 subs[i]++;
wolffd@0 140 if(sy[i])
wolffd@0 141 sp += cpsy[i];
wolffd@0 142 break;
wolffd@0 143 }
wolffd@0 144 }
wolffd@0 145 }
wolffd@0 146 }
wolffd@0 147 else{
wolffd@0 148 for(j=0; j<NB; j++){
wolffd@0 149 *sp += *bp++;
wolffd@0 150 for(i = 0; i < ndim; i++){
wolffd@0 151 if(subs[i] == s[i]){
wolffd@0 152 subs[i] = 0;
wolffd@0 153 if(sy[i])
wolffd@0 154 sp -= cpsy2[i];
wolffd@0 155 }
wolffd@0 156 else{
wolffd@0 157 subs[i]++;
wolffd@0 158 if(sy[i])
wolffd@0 159 sp += cpsy[i];
wolffd@0 160 break;
wolffd@0 161 }
wolffd@0 162 }
wolffd@0 163 }
wolffd@0 164 }
wolffd@0 165
wolffd@0 166 free(sx);
wolffd@0 167 free(sy);
wolffd@0 168 free(s);
wolffd@0 169 free(cpsy);
wolffd@0 170 free(subs);
wolffd@0 171 free(cpsy2);
wolffd@0 172 free(mask);
wolffd@0 173 free(ssize);
wolffd@0 174 }
wolffd@0 175