From 1773851d22f9de8d1f977f37b8f10c8899960ff2 Mon Sep 17 00:00:00 2001 From: nilm Date: Wed, 14 Oct 2009 22:55:33 +0000 Subject: [PATCH] Random stepping git-svn-id: https://bucket.mit.edu/svn/nilm/zoom@7989 ddd99763-3ecb-0310-9145-efcb8ce7c51f --- pc/Makefile | 6 +- pc/dctest.c | 61 +++++++++++++++-- pc/gpib.c | 2 +- pc/mt19937ar.c | 173 +++++++++++++++++++++++++++++++++++++++++++++++++ pc/mt19937ar.h | 9 +++ 5 files changed, 244 insertions(+), 7 deletions(-) create mode 100644 pc/mt19937ar.c create mode 100644 pc/mt19937ar.h diff --git a/pc/Makefile b/pc/Makefile index 7f63217..2fe95f1 100644 --- a/pc/Makefile +++ b/pc/Makefile @@ -10,16 +10,18 @@ zoom.o: serial-util.h read.o: serial-util.h +mt19937ar.o: mt19937ar.h + calibrate.o: serial-util.h gpib.h zoom.h -dctest.o: serial-util.h gpib.h zoom.h +dctest.o: serial-util.h gpib.h zoom.h mt19937ar.h read: read.o serial-util.o calibrate: calibrate.o serial-util.o gpib.o zoom.o dctest: LDFLAGS=-lpthread -dctest: dctest.o serial-util.o gpib.o zoom.o +dctest: dctest.o serial-util.o gpib.o zoom.o mt19937ar.o clean: rm -f *.o read calibrate dctest diff --git a/pc/dctest.c b/pc/dctest.c index def8782..b5ed899 100644 --- a/pc/dctest.c +++ b/pc/dctest.c @@ -16,6 +16,7 @@ #include "math.h" #include #include +#include "mt19937ar.h" #define info(x...) fprintf(stderr,x) @@ -29,6 +30,7 @@ int main(int argc, char *argv[]) char *zoomdev=strdup("/dev/serial/by-id/usb-FTDI_FT232R_USB_UART_A6007wc5-if00-port0"); char *gpibdev=strdup("/dev/serial/by-id/usb-Prologix_Prologix_GPIB-USB_Controller_PXQQY20G-if00-port0"); int rate=500000; + unsigned long seed = 1337; int getopt_index; int zoom, gpib; @@ -36,13 +38,14 @@ int main(int argc, char *argv[]) { "zoom-device", required_argument, NULL, 'Z' }, { "gpib-device", required_argument, NULL, 'G' }, { "rate", required_argument, NULL, 'r' }, + { "seed", required_argument, NULL, 's' }, { "help", no_argument, NULL, 'h' }, { 0, 0, 0, 0} }; int help=0; char c; - while ((c = getopt_long(argc, argv, "Z:G:r:h?", + while ((c = getopt_long(argc, argv, "Z:G:r:s:h?", long_opts, &getopt_index)) != -1) { switch(c) { @@ -57,7 +60,12 @@ int main(int argc, char *argv[]) case 'r': rate = atoi(optarg); if(rate == 0) - errx(1, "invalid rate: %s",optarg); + errx(1, "invalid rate: %s", optarg); + break; + case 's': + seed = atol(optarg); + if(seed == 0) + errx(1, "invalid seed: %s", optarg); break; case 'h': case '?': @@ -73,12 +81,16 @@ int main(int argc, char *argv[]) fprintf(stderr, " -Z, --zoom-device %-9s zoom NILM serial port\n", "/dev/xxx"); fprintf(stderr, " -G, --gpib-device %-9s GPIB serial port\n", "/dev/xxx"); fprintf(stderr, " -r, --rate %-16d baud rate\n", rate); + fprintf(stderr, " -s, --seed %-16ld random seed\n", seed); fprintf(stderr, " -h, --help this help\n"); return 1; } signal(SIGINT, handle_sig); + info("Initializing twister with seed %ld\n", seed); + init_genrand(seed); + /* open devices */ info("Opening Zoom NILM device %s\n", zoomdev); if ((zoom = serial_open(zoomdev, rate)) == -1) @@ -144,7 +156,7 @@ int process_adc_dac(const uint8_t *buf, struct threadinfo_t *ti) pthread_mutex_unlock(&ti->mutex); /* write it out */ - printf("%d %.15f %.15f %.8f %5d %5d %d\n", + printf("%d %.12f %.12f %.8f %5d %5d %d\n", stable, idesired, iactual, @@ -248,6 +260,15 @@ static int keithley_change(int gpib, double desired, struct threadinfo_t *ti) return 0; } +static double genrand(double min, double max) +{ + double x; + x = genrand_real1(); /* [0, 1] */ + x *= (max - min); /* [0, max-min] */ + x += min; /* [min, max] */ + return x; +} + static void dctest(int zoom, int gpib) { pthread_t thread; @@ -311,7 +332,39 @@ static void dctest(int zoom, int gpib) info("Running\n"); /* Change Keithley values */ while (!g_quit) { - sleep(1); + +#define DELAY 100000 /* after keithley settles, in us */ +#define STEPS 10 /* how many small steps to take */ + +#define K_MIN -0.5 /* keithley min, amps */ +#define K_MAX 0.5 /* keithley max, amps */ +#define STEPSIZE 0.03 /* max size of small step, in amps */ + + /* Choose any value within the Keithley range */ + double desired = genrand(K_MIN, K_MAX); + info("Big step: %.12f\n", desired); + if (keithley_change(gpib, desired, &ti) < 0) { + info("error setting keithley\n"); + break; + } + usleep(DELAY); + + /* Choose a few more values nearby */ + for (i = 0; i < STEPS && !g_quit; i++) { + desired += genrand(-STEPSIZE/2, STEPSIZE/2); + if (desired < -1.0) + desired = -1.0; + if (desired > 1.0) + desired = 1.0; + if (keithley_change(gpib, desired, &ti) < 0) { + info("error setting keithley\n"); + break; + } + usleep(DELAY); + } + + /* one extra delay just to separate things a bit */ + usleep(DELAY); } printf("main thread quitting\n"); diff --git a/pc/gpib.c b/pc/gpib.c index 78d79d0..7d7c116 100644 --- a/pc/gpib.c +++ b/pc/gpib.c @@ -50,7 +50,7 @@ int keithley_init(int fd) int keithley_current(int fd, double amps) { char s[128]; - sprintf(s, ":sour:curr %0.3f", amps); + sprintf(s, ":sour:curr %0.12f", amps); gput(s); return 0; } diff --git a/pc/mt19937ar.c b/pc/mt19937ar.c new file mode 100644 index 0000000..30138a1 --- /dev/null +++ b/pc/mt19937ar.c @@ -0,0 +1,173 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +#include "mt19937ar.h" + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void init_by_array(unsigned long init_key[], int key_length) +{ + int i, j, k; + init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,0x7fffffff]-interval */ +long genrand_int31(void) +{ + return (long)(genrand_int32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +double genrand_real1(void) +{ + return genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +double genrand_real2(void) +{ + return genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +double genrand_real3(void) +{ + return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on [0,1) with 53-bit resolution*/ +double genrand_res53(void) +{ + unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} +/* These real versions are due to Isaku Wada, 2002/01/09 added */ + diff --git a/pc/mt19937ar.h b/pc/mt19937ar.h new file mode 100644 index 0000000..e9abc32 --- /dev/null +++ b/pc/mt19937ar.h @@ -0,0 +1,9 @@ +void init_genrand(unsigned long s); +void init_by_array(unsigned long init_key[], int key_length); +unsigned long genrand_int32(void); +long genrand_int31(void); +double genrand_real1(void); +double genrand_real2(void); +double genrand_real3(void); +double genrand_res53(void); +