Changeset 2286 for research


Ignore:
Timestamp:
Apr 16, 2008, 12:11:32 AM (12 years ago)
Author:
Sam Hocevar
Message:
  • More scripts for part 3.
  • Implemented main.c as a seccomp bytecode for CPUShare.
Location:
research/2008-displacement
Files:
3 added
7 edited
2 copied

Legend:

Unmodified
Added
Removed
  • research/2008-displacement/.gitignore

    r2284 r2286  
    44*.log
    55*.pdf
     6*.o
    67xy2d
    78vote
     9bytecode
     10bytecode.*.bin
     11bytecode.lds
     12bytecode.lds.s
  • research/2008-displacement/Makefile

    r2284 r2286  
    1 all: main xy2d vote
     1ifeq ($(shell uname -m), x86_64)
     2ifeq ($(shell [ -d /lib64 ] && echo 64), 64)
     3CFLAGS = -m32 -march=i686 -mcpu=pentium4 -fomit-frame-pointer
     4CPPFLAGS = -m32
     5LDFLAGS = -melf_i386
     6endif
     7else # ! x86_64
     8
     9ifeq ($(shell uname -m), i686)
     10CFLAGS = -march=i686 -mcpu=i686 -fomit-frame-pointer
     11CPPFLAGS = -m32
     12LDFLAGS = -melf_i386
     13else # ! i686
     14
     15ifeq ($(shell uname -m), ppc64)
     16# mcpu means march and mtune means mcpu, oh well
     17CFLAGS = -m32 -mcpu=750 -mtune=970
     18CPPFLAGS = -m32
     19LDFLAGS = -melf32ppclinux
     20else # ! ppc64
     21
     22ifeq ($(shell uname -m), ppc)
     23# mcpu means march and mtune means mcpu, oh well
     24CFLAGS = -mcpu=750 -mtune=970
     25CPPFLAGS = -m32
     26LDFLAGS = -melf32ppclinux
     27endif # ppc
     28
     29endif # ppc64
     30endif # i686
     31endif # x86_64
     32
     33CFLAGS += -fno-common -O2 -Wall -Iglibc -I.
     34CPPFLAGS += -Iglibc -I.
     35LDFLAGS += -O2
     36
     37CC = $(CROSS_COMPILE)gcc
     38CPP = $(CROSS_COMPILE)cpp
     39LD = $(CROSS_COMPILE)ld
     40
     41OBJCOPY = $(CROSS_COMPILE)objcopy
     42
     43LIBGCC = $(shell $(CC) $(CFLAGS) -print-libgcc-file-name)
     44LIBGCCEH = $(shell $(CC) $(CFLAGS) -print-file-name=libgcc_eh.a)
     45LIBM = $(shell $(CC) $(CFLAGS) -print-file-name=libm.a)
     46LIBC = $(shell $(CC) $(CFLAGS) -print-file-name=libc.a)
     47
     48all: bytecode.text.bin bytecode.data.bin main xy2d vote
     49
     50bytecode.o: main.c
     51        $(CC) -DBYTECODE -c $(CFLAGS) $< -o $@
     52
     53bytecode.lds.s: bytecode.lds.S
     54        $(CPP) $(CPPFLAGS) $< -o $@
     55
     56bytecode.lds: bytecode.lds.s
     57        grep -A100000000 SECTION $< > $@
     58
     59bytecode: bytecode.o bytecode.lds
     60        $(LD) $(LDFLAGS) -static -T bytecode.lds $< --start-group $(LIBGCC) $(LIBC) $(LIBM) $(LIBGCCEH) --end-group -o $@
     61
     62bytecode.text.bin: bytecode
     63        $(OBJCOPY) -O binary $< -j .text $@
     64
     65bytecode.data.bin: bytecode
     66        $(OBJCOPY) -O binary $< -j .data $@
     67
     68clean:
     69        rm -f bytecode bytecode.lds bytecode.lds.s
     70        rm -f xy2d vote main
     71        rm -f *.pyc *.bin *.o
    272
    373xy2d: xy2d.c
     
    1080        $(CC) -Wall -O3 -funroll-loops -ffast-math -W -Wall $^ -o $@ $$(pkg-config --cflags --libs sdl) -lSDL_image -lm
    1181
    12 clean:
    13         rm -f xy2d vote main
    14 
     82.PHONY: clean
  • research/2008-displacement/README

    r2280 r2286  
    1 # Trouver des images au pif
    2 % find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9].???' | rev | sort -k2 -t. | rev | xargs -n 1 ./main | tee fs-4chan.txt
    3 % find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9].???' | rev | sort -k2 -t. | rev | xargs -n 1 ./main-jajuni | tee jajuni-4chan.txt
    4 % cat /tmp/4chanlist.txt | xargs -n 1 ./main | tee -a fs-4chan.txt
     1# List all my 4chan images
     2find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9].???' | rev | sort -k2 -t. | rev > 4chan-list.txt
     3
     4# Put all my 4chan images in 100 separate /tmp directories
     5for x in $(seq -w 00 09); do echo $x; mkdir -p /tmp/4chan/$x; cp $(find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]'$x'.???') /tmp/4chan/$x; done
     6
     7# Condorcet voting for phase 2 results
     8#  - raster + E
     9#  - raster + E_min
     10#  - serpentine + E
     11#  - serpentine + E_min
     12for x in part2/*txt ; do grep '^.1' $x | awk '{ print $3,$5 }' >| $x.clean; done ; ./vote part2/*clean | sort -rnk3 | head -30 ; rm -f part2/*clean
     13for x in part2/*txt ; do grep '^.1' $x | awk '{ print $3,$9 }' >| $x.clean; done ; ./vote part2/*clean | sort -rnk3 | head -30 ; rm -f part2/*clean
     14for x in part2/*txt ; do grep '^.2' $x | awk '{ print $3,$5 }' >| $x.clean; done ; ./vote part2/*clean | sort -rnk3 | head -30 ; rm -f part2/*clean
     15for x in part2/*txt ; do grep '^.2' $x | awk '{ print $3,$9 }' >| $x.clean; done ; ./vote part2/*clean | sort -rnk3 | head -30 ; rm -f part2/*clean
     16
     17# Get phase 3 and phase 4 stuff
     18ssh canard.zoy.org "cd test-20080329; tar cz *raster.txt *serp.txt" | tar xz
     19for x in *-raster.txt; do y="$x"; y="${y%%-raster.txt}"; y="${y%%.tiff}"; y="${y##usc-sipi}"; \mv "$x" part3/"$y".txt; done
     20for x in *-serp.txt; do y="$x"; y="${y%%-serp.txt}"; y="${y%%.tiff}"; y="${y##usc-sipi}"; \mv "$x" part4/"$y".txt; done
     21
     22# Condorcet voting for part 3 and 4
     23for x in part3/*txt ; do cat $x | awk '{ print $2,$4 }' >| $x.clean; done ; ./vote part3/*clean | sort -rnk3 | head -20 ; rm -f part3/*clean
     24for x in part3/*txt ; do cat $x | awk '{ print $2,$8 }' >| $x.clean; done ; ./vote part3/*clean | sort -rnk3 | head -20 ; rm -f part3/*clean
     25for x in part4/*txt ; do cat $x | awk '{ print $2,$4 }' >| $x.clean; done ; ./vote part4/*clean | sort -rnk3 | head -20 ; rm -f part4/*clean
     26for x in part4/*txt ; do cat $x | awk '{ print $2,$8 }' >| $x.clean; done ; ./vote part4/*clean | sort -rnk3 | head -20 ; rm -f part4/*clean
     27
     28# Mean voting for part 3 and 4
     29cat part3/aerials2.1.01.txt | while read x k y ; do echo "$k $(grep $k part3/* | awk '{ a+=$4; dx+=$10; dy+=$12; n+=1 } END { print a/n, dx/n, dy/n }')"; done | sort -nk2 | head -20
     30cat part3/aerials2.1.01.txt | while read x k y ; do echo "$k $(grep $k part3/* | awk '{ a+=$8; dx+=$10; dy+=$12; n+=1 } END { print a/n, dx/n, dy/n }')"; done | sort -nk2 | head -20
     31cat part4/aerials2.1.01.txt | while read x k y ; do echo "$k $(grep $k part4/* | awk '{ a+=$4; dx+=$10; dy+=$12; n+=1 } END { print a/n, dx/n, dy/n }')"; done | sort -nk2 | head -20
     32cat part4/aerials2.1.01.txt | while read x k y ; do echo "$k $(grep $k part4/* | awk '{ a+=$8; dx+=$10; dy+=$12; n+=1 } END { print a/n, dx/n, dy/n }')"; done | sort -nk2 | head -20
     33
     34# Clever stuff (or not)
     35cat part3/aerials2.1.01.txt | grep K | while read x k y ; do grep $k part3/* | awk '{ dx+=$10; dy+=$12; n+=1 } END { print dx/n, dy/n }' | read dx dy; echo "$k $(grep $k part3/* | awk 'BEGIN { dx='$dx'; dy='$dy' } { x=dx-$10; y=dy-$12; t+=2.^-5*(x*x+y*y); a+=t*$4; n+=t } END { print a/n, n }')"; done | sort -nk2 | head -20
     36cat part3/aerials2.1.01.txt | grep K | while read x k y ; do grep $k part3/* | awk '{ dx+=$10; dy+=$12; n+=1 } END { print dx/n, dy/n }' | read dx dy; echo "$k $(grep $k part3/* | awk 'BEGIN { dx='$dx'; dy='$dy' } { x=dx-$10; y=dy-$12; t+=2.^-5*(x*x+y*y); a+=t*$8; n+=t } END { print a/n, n }')"; done | sort -nk2 | head -20
     37cat part4/aerials2.1.01.txt | grep K | while read x k y ; do grep $k part4/* | awk '{ dx+=$10; dy+=$12; n+=1 } END { print dx/n, dy/n }' | read dx dy; echo "$k $(grep $k part4/* | awk 'BEGIN { dx='$dx'; dy='$dy' } { x=dx-$10; y=dy-$12; t+=2.^-5*(x*x+y*y); a+=t*$4; n+=t } END { print a/n, n }')"; done | sort -nk2 | head -20
     38cat part4/aerials2.1.01.txt | grep K | while read x k y ; do grep $k part4/* | awk '{ dx+=$10; dy+=$12; n+=1 } END { print dx/n, dy/n }' | read dx dy; echo "$k $(grep $k part4/* | awk 'BEGIN { dx='$dx'; dy='$dy' } { x=dx-$10; y=dy-$12; t+=2.^-5*(x*x+y*y); a+=t*$8; n+=t } END { print a/n, n }')"; done | sort -nk2 | head -20
     39
     40###
     41###
     42###
     43###
     44###  Stuff below here is deprecated or unsorted
     45###
     46###
     47###
     48###
     49
     50#% find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9].???' | rev | sort -k2 -t. | rev | xargs -n 1 ./main | tee fs-4chan.txt
     51#% find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9].???' | rev | sort -k2 -t. | rev | xargs -n 1 ./main-jajuni | tee jajuni-4chan.txt
     52#% cat /tmp/4chanlist.txt | xargs -n 1 ./main | tee -a fs-4chan.txt
    553
    654# Lena
     
    190238for x in out-*.txt; do sort -k7 $x | head -20 ; done | cut -f1 -d: | sort | uniq -c | sort -n
    191239
    192 ### Pour faire des répertoires
    193 for x in $(seq -w 00 09); do echo $x; mkdir -p /tmp/4chan/$x; cp $(find ~/4chan/unsorted-4chan/http* -name '1[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]'$x'.???') /tmp/4chan/$x; done
    194 
    195 ### Conversion du vieux format de résultats
     240###
     241### Conversion du vieux format de résultats (deprecated)
    196242cat -n $x | sed 's/^  *[0-9]*\([0-9]\)[[:space:]]*/[\1] /; s/....###/###/; s/\[[27]/[1/; s/\[[38]/[2/; s/\[[49]/[3/; s/\[[50]/[4/'
    197243
     244#
  • research/2008-displacement/main.c

    r2281 r2286  
    44#include <stdlib.h>
    55#include <string.h>
    6 #include <stdlib.h>
     6#include <stdarg.h>
    77
    88#include <math.h>
    99
    10 #include <SDL_image.h>
     10#ifdef BYTECODE
     11#   include <seccomp-loader.h>
     12#else
     13#   include <SDL_image.h>
     14#endif
    1115
    1216#define MAXWIDTH 512
     
    1822static int WIDTH, HEIGHT;
    1923
    20 static inline float get(float const *src, int x, int y)
     24int main(int, char *[]);
     25
     26#ifdef BYTECODE
     27#   define MAXIMAGES 128
     28static int slots[MAXIMAGES];
     29static double slotbuf[MAXIMAGES * MAXWIDTH * MAXHEIGHT];
     30
     31volatile int stop;
     32
     33void sighandler(int signal)
     34{
     35    if(sys_write(1, "done", 4) != 4)
     36        sys_exit(3);
     37    stop = 1;
     38}
     39
     40void bytecode(unsigned char * mem, int heap_size, int stack_size)
     41{
     42    char mode[] = "-0";
     43    char *argv[] = { "program", mode, "arg" };
     44    char c;
     45
     46    if(sys_read(0, mode + 1, 1) != 1)
     47        sys_exit(-5);
     48
     49    main(sizeof(argv)/sizeof(*argv), argv);
     50
     51    c = 0;
     52    if(sys_write(1, &c, 1) != 1)
     53        sys_exit(-6);
     54}
     55
     56static int myatoi(const char *str)
     57{
     58    int i = 0;
     59    while(*str)
     60    {
     61        i *= 10;
     62        i += (unsigned char)*str++ - (unsigned char)'0';
     63    }
     64    return i;
     65}
     66#define atoi myatoi
     67#endif
     68
     69#define WRITE_INT(i, base) \
     70    do \
     71    { \
     72        char buf[128], *b = buf + 127; \
     73        if(i <= 0) \
     74            sys_write(1, (i = -i) ? "-" : "0", 1); /* XXX: hack here */ \
     75        while(i) \
     76        { \
     77            *b-- = hex2char[i % base]; \
     78            i /= base; \
     79        } \
     80        sys_write(1, b + 1, (int)(buf + 127 - b)); \
     81    } while(0)
     82
     83static void msg(const char *f, ...)
     84{
     85    va_list args;
     86    va_start(args, f);
     87#ifdef BYTECODE
     88    static char const *hex2char = "0123456789abcdef";
     89
     90    for( ; *f; f++)
     91    {
     92        if(*f != '%')
     93        {
     94            sys_write(1, f, 1);
     95            continue;
     96        }
     97
     98        f++;
     99        if(!*f)
     100            break;
     101
     102        if(*f == 'c')
     103        {
     104            char i = (char)(unsigned char)va_arg(args, int);
     105            if(i >= 0x20 && i < 0x7f)
     106                sys_write(1, &i, 1);
     107            else if(i == '\n')
     108                sys_write(1, "\\n", 2);
     109            else if(i == '\t')
     110                sys_write(1, "\\t", 2);
     111            else if(i == '\r')
     112                sys_write(1, "\\r", 2);
     113            else
     114            {
     115                sys_write(1, "\\x", 2);
     116                sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
     117                sys_write(1, hex2char + (i & 0x0f), 1);
     118            }
     119        }
     120        else if(*f == 'i' || *f == 'd')
     121        {
     122            int i = va_arg(args, int);
     123            WRITE_INT(i, 10);
     124        }
     125        else if(*f == 'x')
     126        {
     127            int i = va_arg(args, int);
     128            WRITE_INT(i, 16);
     129        }
     130        else if(f[0] == 'l' && (f[1] == 'i' || f[1] == 'd'))
     131        {
     132            long int i = va_arg(args, long int);
     133            WRITE_INT(i, 10);
     134            f++;
     135        }
     136        else if(f[0] == 'l' && f[1] == 'l' && (f[2] == 'i' || f[1] == 'd'))
     137        {
     138            long long int i = va_arg(args, long long int);
     139            WRITE_INT(i, 10);
     140            f += 2;
     141        }
     142        else if(f[0] == 'g')
     143        {
     144            double g = va_arg(args, double), h = 9.9f;
     145            int i;
     146            if(g < 0.f)
     147            {
     148                sys_write(1, "-", 1);
     149                g = -g;
     150            }
     151            for(i = 0; i < 7; i++)
     152            {
     153                sys_write(1, hex2char + (int)g, 1);
     154                if(i == 0)
     155                    sys_write(1, ".", 1);
     156                g = (g - (int)g) * 10;
     157                h = h / 10.f;
     158                if(g < h)
     159                    break;
     160            }
     161        }
     162        else if(f[0] == 'p')
     163        {
     164            uintptr_t i = va_arg(args, uintptr_t);
     165            if(!i)
     166                sys_write(1, "NULL", 5);
     167            else
     168            {
     169                sys_write(1, "0x", 2);
     170                WRITE_INT(i, 16);
     171            }
     172        }
     173        else if(f[0] == 's')
     174        {
     175            char *s = va_arg(args, char *);
     176            if(!s)
     177                sys_write(1, "(nil)", 5);
     178            else
     179            {
     180                int l = 0;
     181                while(s[l])
     182                    l++;
     183                sys_write(1, s, l);
     184            }
     185        }
     186        else if(f[0] == '0' && f[1] == '2' && f[2] == 'x')
     187        {
     188            int i = va_arg(args, int);
     189            sys_write(1, hex2char + ((i & 0xf0) >> 4), 1);
     190            sys_write(1, hex2char + (i & 0x0f), 1);
     191            f += 2;
     192        }
     193        else
     194        {
     195            sys_write(1, f - 1, 2);
     196        }
     197    }
     198#else
     199    vprintf(f, args);
     200    fflush(stdout);
     201#endif
     202    va_end(args);
     203}
     204
     205static inline double get(double const *src, int x, int y)
    21206{
    22207    return src[y * WIDTH + x];
    23208}
    24209
    25 static inline void put(float *src, int x, int y, float p)
     210static inline void put(double *src, int x, int y, double p)
    26211{
    27212    src[y * WIDTH + x] = p;
    28213}
    29214
    30 static float *new(void)
    31 {
    32     return malloc(WIDTH * HEIGHT * sizeof(float));
    33 }
    34 
    35 static float *copy(float const *src)
    36 {
    37     float *dest = malloc(WIDTH * HEIGHT * sizeof(float));
    38     memcpy(dest, src, WIDTH * HEIGHT * sizeof(float));
     215static double *new(void)
     216{
     217#ifdef BYTECODE
     218    int i;
     219    for(i = 0; i < MAXIMAGES; i++)
     220        if(slots[i] == 0)
     221            break;
     222    if(i == MAXIMAGES)
     223        return NULL;
     224    slots[i] = 1;
     225    return slotbuf + i * MAXWIDTH * MAXHEIGHT;
     226#else
     227    return malloc(WIDTH * HEIGHT * sizeof(double));
     228#endif
     229}
     230
     231static void del(double *img)
     232{
     233#ifdef BYTECODE
     234    int i;
     235    for(i = 0; i < MAXIMAGES; i++)
     236        if(slotbuf + i * MAXWIDTH * MAXHEIGHT == img)
     237            break;
     238    if(i == MAXIMAGES)
     239        return;
     240    slots[i] = 0;
     241#else
     242    free(img);
     243#endif
     244}
     245
     246static double *copy(double const *src)
     247{
     248    double *dest = new();
     249    memcpy(dest, src, WIDTH * HEIGHT * sizeof(double));
    39250    return dest;
    40251}
     
    43254#define NN ((N * 2 + 1))
    44255
    45 static void makegauss(float mat[NN][NN], float sigma, float dx, float dy)
    46 {
    47     float t = 0;
     256static void makegauss(double mat[NN][NN], double sigma, double dx, double dy)
     257{
     258    double t = 0;
    48259    int i, j;
    49260
     
    53264        for(i = 0; i < NN; i++)
    54265        {
    55             float a = (float)(i - N) + dx;
    56             float b = (float)(j - N) + dy;
     266            double a = (double)(i - N) + dx;
     267            double b = (double)(j - N) + dy;
    57268            mat[i][j] = pow(M_E, - (a * a + b * b) / sigma);
    58269            t += mat[i][j];
     
    64275}
    65276
    66 static float *gauss(float const *src, float mat[NN][NN])
    67 {
    68     float *dest = new();
     277static double *gauss(double const *src, double mat[NN][NN])
     278{
     279    double *dest = new();
    69280    int x, y, i, j;
    70281
     
    72283        for(x = N; x < WIDTH - N; x++)
    73284        {
    74             float p = 0;
     285            double p = 0;
    75286
    76287            for(j = 0; j < NN; j++)
     
    84295}
    85296
    86 static float *fullgauss(float const *src, float mat[NN][NN])
    87 {
    88     float *dest = new();
     297static double *fullgauss(double const *src, double mat[NN][NN])
     298{
     299    double *dest = new();
    89300    int x, y, i, j;
    90301
     
    92303        for(x = 0; x < WIDTH; x++)
    93304        {
    94             float p = 0;
     305            double p = 0;
    95306
    96307            for(j = 0; j < NN; j++)
     
    107318
    108319#if 0
    109 static float fulldist(float const *p1, const float *p2)
    110 {
    111     float error = 0.;
     320static double fulldist(double const *p1, const double *p2)
     321{
     322    double error = 0.;
    112323    int x, y;
    113324
     
    115326        for(x = 0; x < WIDTH; x++)
    116327        {
    117             float t = get(p1, x, y) - get(p2, x, y);
     328            double t = get(p1, x, y) - get(p2, x, y);
    118329            error += t * t;
    119330        }
     
    123334#endif
    124335
    125 static float dist(float const *p1, float const *p2, float max)
    126 {
    127     float error = 0.;
     336static double dist(double const *p1, double const *p2, double max)
     337{
     338    double error = 0.;
    128339    int x, y;
    129340
     
    134345        for(x = N; x < WIDTH - N; x++)
    135346        {
    136             float t = get(p1, x, y) - get(p2, x, y);
     347            double t = get(p1, x, y) - get(p2, x, y);
    137348            error += t * t;
    138349        }
     
    148359}
    149360
    150 static float *load(char const *name)
    151 {
     361static double *load(char const *name)
     362{
     363    double *floats;
     364    int x, y;
     365
     366#ifdef BYTECODE
     367    char c;
     368
     369    if(sys_read(0, &c, 1) != 1)
     370        sys_exit(-5);
     371    WIDTH = ((int)(unsigned char)c) << 8;
     372    if(sys_read(0, &c, 1) != 1)
     373        sys_exit(-5);
     374    WIDTH |= (int)(unsigned char)c;
     375
     376    if(sys_read(0, &c, 1) != 1)
     377        sys_exit(-5);
     378    HEIGHT = ((int)(unsigned char)c) << 8;
     379    if(sys_read(0, &c, 1) != 1)
     380        sys_exit(-5);
     381    HEIGHT |= (int)(unsigned char)c;
     382
     383    floats = new();
     384    if(!floats)
     385        return NULL;
     386
     387    for(y = 0; y < HEIGHT; y++)
     388    for(x = 0; x < WIDTH; x++)
     389    {
     390        if(sys_read(0, &c, 1) != 1)
     391            sys_exit(-5);
     392        put(floats, x, y, (double)(unsigned char)c / 0xff);
     393    }   
     394#else
    152395    SDL_Surface *tmp, *surface;
    153396    uint32_t *pixels;
    154     float *floats;
    155     int x, y;
    156397
    157398    tmp = IMG_Load(name);
     
    161402    WIDTH = tmp->w > MAXWIDTH ? MAXWIDTH : tmp->w;
    162403    HEIGHT = tmp->h > MAXHEIGHT ? MAXHEIGHT : tmp->h;
    163     floats = malloc(WIDTH * HEIGHT * sizeof(float));
     404    floats = new();
    164405    if(!floats)
    165406        return NULL;
     
    172413
    173414    for(y = 0; y < HEIGHT; y++)
    174         for(x = 0; x < WIDTH; x++)
    175         {
    176             int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
    177             put(floats, x, y, (float)green / 0xff);
    178         }
     415    for(x = 0; x < WIDTH; x++)
     416    {
     417        int green = (pixels[y * surface->pitch / 4 + x] >> 8) & 0xff;
     418        put(floats, x, y, (double)green / 0xff);
     419    }
     420#endif
    179421
    180422    return floats;
    181423}
    182424
    183 static void save(float const *src, char const *name)
    184 {
     425static void save(double const *src, char const *name)
     426{
     427    int x, y;
     428#ifdef BYTECODE
     429    for(y = 0; y < HEIGHT; y++)
     430    for(x = 0; x < WIDTH; x++)
     431    {
     432        double p = 255 * get(src, x, y);
     433        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
     434        char c = (char)(unsigned char)i;
     435        if(sys_write(1, &c, 1) != 1)
     436            sys_exit(-6);
     437    }
     438#else
    185439    SDL_Surface *surface;
    186440    uint32_t *pixels;
    187     int x, y;
    188441
    189442    surface = SDL_CreateRGBSurface(SDL_SWSURFACE, WIDTH, HEIGHT, 32,
     
    194447    for(x = 0; x < WIDTH; x++)
    195448    {
    196         float p = 255 * get(src, x, y);
     449        double p = 255 * get(src, x, y);
    197450        uint32_t i = p < 0 ? 0 : p > 255 ? 255 : p;
    198451        pixels[surface->pitch / 4 * y + x] = (i << 16) | (i << 8) | i;
     
    200453
    201454    SDL_SaveBMP(surface, name);
    202 }
    203 
    204 static float *ostromoukhov(float const *src)
     455#endif
     456}
     457
     458static double *ostromoukhov(double const *src)
    205459{
    206460    static int const table[][3] =
     
    240494    };
    241495
    242     float *dest = new();
    243     float *tmp = copy(src);
     496    double *dest = new();
     497    double *tmp = copy(src);
    244498    int x, y;
    245499
     
    252506            int x3 = (y & 1) ? WIDTH - 1 - x - 1 : x + 1;
    253507
    254             float p = get(tmp, x2, y);
    255             float q = p > 0.5 ? 1. : 0.;
    256             float error = (p - q);
     508            double p = get(tmp, x2, y);
     509            double q = p > 0.5 ? 1. : 0.;
     510            double error = (p - q);
    257511            int i = p * 255.9999;
    258512
     
    278532    }
    279533
    280     free(tmp);
     534    del(tmp);
    281535
    282536    return dest;
     
    288542     h i j k l
    289543*/
    290 static float *ed(float const *src, int serpentine,
     544static double *ed(double const *src, int serpentine,
    291545                 int a, int b, int c, int d, int e, int f,
    292546                 int g, int h, int i, int j, int k, int l)
    293547{
    294     float *dest = new();
    295     float *tmp = copy(src);
     548    double *dest = new();
     549    double *tmp = copy(src);
    296550    int x, y, n;
    297551
     
    310564            int x4 = swap ? WIDTH - 1 - x - 2 : x + 2;
    311565
    312             float p = get(tmp, x2, y);
    313             float q = p > 0.5 ? 1. : 0.;
    314             float error = (p - q) / n;
     566            double p = get(tmp, x2, y);
     567            double q = p > 0.5 ? 1. : 0.;
     568            double error = (p - q) / n;
    315569
    316570            put(dest, x2, y, q);
     
    363617    }
    364618
    365     free(tmp);
     619    del(tmp);
    366620
    367621    return dest;
     
    369623
    370624/* XXX */
    371 static float *dbs(float const *src, float const *orig,
    372                   float sigma, float dx, float dy)
    373 {
    374     float mat[NN][NN];
    375     float *dest, *tmp, *tmp2;
    376     float error;
     625static double *dbs(double const *src, double const *orig,
     626                  double sigma, double dx, double dy)
     627{
     628    double mat[NN][NN];
     629    double *dest, *tmp, *tmp2;
     630    double error;
    377631
    378632    makegauss(mat, sigma, 0., 0.);
     
    393647        for(x = 0; x < WIDTH; x++)
    394648        {
    395             float d, d2, e, best = 0.;
     649            double d, d2, e, best = 0.;
    396650            int opx = -1, opy = -1;
    397651
     
    407661                for(i = -N; i < N + 1; i++)
    408662                {
    409                     float m, p, q1, q2;
     663                    double m, p, q1, q2;
    410664
    411665                    if(x + i < 0 || x + i >= WIDTH)
     
    446700                    for(i = -N; i < N + 1; i++)
    447701                    {
    448                         float ma, mb, p, q1, q2;
     702                        double ma, mb, p, q1, q2;
    449703                        if(x + i < 0 || x + i >= WIDTH)
    450704                            continue;
     
    481735            for(i = -N; i < N + 1; i++)
    482736            {
    483                 float m = mat[i + N][j + N];
     737                double m = mat[i + N][j + N];
    484738                if(y + j >= 0 && y + j < HEIGHT
    485739                    && x + i >= 0 && x + i < WIDTH)
    486740                {
    487                     float t = get(tmp2, x + i, y + j);
     741                    double t = get(tmp2, x + i, y + j);
    488742                    put(tmp2, x + i, y + j, t + m * (d2 - d));
    489743                }
     
    491745                                && x + opx + i >= 0 && x + opx + i < WIDTH)
    492746                {
    493                     float t = get(tmp2, x + opx + i, y + opy + j);
     747                    double t = get(tmp2, x + opx + i, y + opy + j);
    494748                    put(tmp2, x + opx + i, y + opy + j, t + m * (d - d2));
    495749                }
     
    505759    }
    506760
    507     free(tmp);
    508     free(tmp2);
     761    del(tmp);
     762    del(tmp2);
    509763
    510764    return dest;
    511765}
    512766
    513 static void study(float const *src, float const *dest,
    514                   float sigma, float precision, float fdx, float fdy)
     767static void study(double const *src, double const *dest,
     768                  double sigma, double precision, double fdx, double fdy)
    515769{
    516770#   define Z 3
    517     float mat[NN][NN], mat2[NN][NN];
    518     float *tmp, *tmp2;
    519     float e, e0, e1;
    520     float best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
     771    double mat[NN][NN], mat2[NN][NN];
     772    double *tmp, *tmp2;
     773    double e, e0, e1;
     774    double best = 1., fx = -1., fy = -1., step = 2., bfx = 0., bfy = 0.;
    521775    int dx, dy;
    522776
     
    526780    tmp2 = gauss(dest, mat);
    527781    e0 = dist(tmp, tmp2, 1.);
    528     free(tmp2);
     782    del(tmp2);
    529783
    530784    makegauss(mat2, sigma, fdx, fdy);
    531785    tmp2 = gauss(dest, mat2);
    532786    e1 = dist(tmp, tmp2, 1.);
    533     free(tmp2);
     787    del(tmp2);
    534788
    535789    while(step > precision)
     
    541795                tmp2 = gauss(dest, mat);
    542796                e = dist(tmp, tmp2, best);
    543                 free(tmp2);
     797                del(tmp2);
    544798                if(e < best)
    545799                {
     
    555809    }
    556810
    557     free(tmp);
    558 
    559     printf("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n",
    560            1000. * e0, 1000. * e1, 1000. * best, fx, fy);
    561     fflush(stdout);
    562 }
    563 
    564 static float *merge(float const *im1, float const *im2, float t)
    565 {
    566     float *dest = new();
     811    del(tmp);
     812
     813    msg("E: %g E_fast: %g E_min: %g dx: %g dy: %g\n",
     814        1000. * e0, 1000. * e1, 1000. * best, fx, fy);
     815}
     816
     817static double *merge(double const *im1, double const *im2, double t)
     818{
     819    double *dest = new();
    567820    int x, y;
    568821
     
    587840int main(int argc, char *argv[])
    588841{
    589     float *src;
     842    double *src;
    590843    int mode, i;
    591844
     
    608861        return 2;
    609862
    610     printf("### mode %i on `%s' ###\n", mode, argv[2]);
     863    msg("### mode %i on `%s' ###\n", mode, argv[2]);
    611864
    612865    switch(mode)
     
    614867        case 1:
    615868        {
    616             float *dest = ed(src, false, 7, 0,
     869            double *dest = ed(src, false, 7, 0,
    617870                                0, 3, 5, 1, 0,
    618871                                0, 0, 0, 0, 0);
    619872            study(src, dest, 1.2, 0.001, .16, .28);
    620             free(dest);
    621             free(src);
     873            del(dest);
     874            del(src);
    622875        }
    623876        break;
     
    625878        case 2:
    626879        {
    627             float *src2, *dest, *tmp;
     880            double *src2, *dest, *tmp;
    628881
    629882            if(argc < 3)
     
    635888            for(i = 0; i <= 100; i++)
    636889            {
    637                 tmp = merge(src, src2, (float)i / 100.);
     890                tmp = merge(src, src2, (double)i / 100.);
    638891                dest = ed(tmp, false, 7, 0,
    639892                             0, 3, 5, 1, 0,
    640893                             0, 0, 0, 0, 0);
    641894                study(tmp, dest, 1.2, 0.001, .16, .28);
    642                 free(dest);
    643                 free(tmp);
    644             }
    645             free(src2);
    646             free(src);
     895                del(dest);
     896                del(tmp);
     897            }
     898            del(src2);
     899            del(src);
    647900        }
    648901        break;
     
    651904        case 4:
    652905        {
    653             float mat[NN][NN];
    654             float *dest, *tmp, *tmp2;
     906            double mat[NN][NN];
     907            double *dest, *tmp, *tmp2;
    655908            int a, b, c, d, e;
    656909
     
    689942                for(i = 1; i <= 2; i++)
    690943                {
    691                     //printf("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2);
    692                     printf("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2);
     944                    //msg("[%i] K: %d,%d,%d,%d,%d ", i, a2, b2, c2, d2, e2);
     945                    msg("[%i] K: %d,%d,%d,%d ", i, a2, c2, d2, e2);
    693946
    694947                    dest = ed(src, i == 2, a2, 0,
     
    704957                        tmp = gauss(src, mat);
    705958                        tmp2 = gauss(dest, mat);
    706                         printf("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.));
    707                         free(tmp);
    708                         free(tmp2);
     959                        msg("E: %.5g\n", 1000. * dist(tmp, tmp2, 1.));
     960                        del(tmp);
     961                        del(tmp2);
    709962                    }
    710963                    fflush(stdout);
    711                     free(dest);
     964                    del(dest);
    712965                }
    713966            }
    714967
    715             free(src);
     968            del(src);
    716969        }
    717970        break;
     
    719972        case 5:
    720973        {
    721             float *dest;
     974            double *dest;
    722975
    723976            dest = ed(src, false, 7, 0,
    724977                         0, 3, 5, 1, 0,
    725978                         0, 0, 0, 0, 0);
    726             printf("[1] ");
     979            msg("[1] ");
    727980            study(src, dest, 1.2, 0.001, .16, .28);
    728             free(dest);
     981            del(dest);
    729982
    730983            dest = ed(src, false, 7, 5,
    731984                         3, 5, 7, 5, 3,
    732985                         1, 3, 5, 3, 1);
    733             printf("[2] ");
     986            msg("[2] ");
    734987            study(src, dest, 1.2, 0.001, .26, .76);
    735             free(dest);
     988            del(dest);
    736989
    737990            dest = ostromoukhov(src);
    738             printf("[3] ");
     991            msg("[3] ");
    739992            study(src, dest, 1.2, 0.001, .0, .19);
    740             free(dest);
     993            del(dest);
    741994
    742995            dest = ed(src, true, 2911, 0,
    743996                  0, 1373, 3457, 2258, 0,
    744997                  0,    0,    0,    0, 0);
    745             printf("[4] ");
     998            msg("[4] ");
    746999            study(src, dest, 1.2, 0.001, .0, .34);
     1000        }
     1001        break;
     1002
     1003        case 6:
     1004        case 7:
     1005        {
     1006            double mat[NN][NN];
     1007            double *dest;
     1008            int a, ax, b, bx, c, cx, d, dx;
     1009
     1010            if(mode == 6)
     1011                a = 7, b = 3, c = 5, d = 1;
     1012            else
     1013                a = 7, b = 5, c = 4, d = 0;
     1014
     1015#undef TOTAL
     1016#define TOTAL (a+b+c+d)
     1017            makegauss(mat, 1.2, 0, 0);
     1018            for(ax = -2; ax <= 2; ax++)
     1019            for(bx = -2; bx <= 2; bx++)
     1020            for(cx = -2; cx <= 2; cx++)
     1021            for(dx = -2; dx <= 3; dx++)
     1022            {
     1023                int a2 = a + ax, b2 = b + bx, c2 = c + cx, d2 = d + dx;
     1024
     1025                if(a2 < 0 || a2 > TOTAL || b2 < 0 || b2 > TOTAL ||
     1026                     c2 < 0 || c2 > TOTAL || d2 < 0 || d2 > TOTAL)
     1027                    continue;
     1028
     1029                if(a2 + b2 + c2 + d2 != TOTAL)
     1030                    continue;
     1031
     1032                msg("K: %d,%d,%d,%d ", a2, b2, c2, d2);
     1033
     1034                dest = ed(src, mode == 7, a2, 0,
     1035                               0, b2, c2, d2, 0,
     1036                               0,  0,  0,  0, 0);
     1037                /* XXX: E_fast is meaningless, ignore it */
     1038                study(src, dest, 1.2, 0.001, 0., 0.);
     1039                fflush(stdout);
     1040                del(dest);
     1041            }
     1042
     1043            del(src);
    7471044        }
    7481045        break;
     
    7541051    //dest = dbs(src, tmp, 0.158718, 0.283089);
    7551052    //dest = copy(tmp);
    756     free(tmp);
     1053    del(tmp);
    7571054    study(src, dest, 1.2, 0.00001);
    7581055    save(dest, "output.bmp");
    759     free(dest);
     1056    del(dest);
    7601057#endif
    7611058
     
    7721069    dest = ostromoukhov(src);
    7731070    //dest = ed(src, 7, 0, 0, 4, 5, 0, 0, 0, 0, 0, 0, 0);
    774     //printf("%s: ", argv[1]);
     1071    //msg("%s: ", argv[1]);
    7751072    //study(src, dest, 1.2, 0.001);
    7761073    save(dest, "output.bmp");
    777     free(dest);
     1074    del(dest);
    7781075#endif
    7791076
     
    7871084        for(dx = 0; dx < STEP; dx++)
    7881085        {
    789             float fy = 2. / STEP * (dy - STEP / 2.);
    790             float fx = 2. / STEP * (dx - STEP / 2.);
     1086            double fy = 2. / STEP * (dy - STEP / 2.);
     1087            double fx = 2. / STEP * (dx - STEP / 2.);
    7911088
    7921089            makegauss(mat, 1.2, fx, fy);
    7931090            tmp2 = gauss(dest, mat);
    794             printf("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
     1091            msg("%g %g %g\n", fy, fx, 1000. * dist(tmp, tmp2, 1.));
    7951092            fflush(stdout);
    796             free(tmp2);
    797         }
    798         printf("\n");
     1093            del(tmp2);
     1094        }
     1095        msg("\n");
    7991096    }
    8001097
     
    8131110        dest = ed(src, a, 0, 0, b, c, d, 0, 0, 0, 0, 0, 0);
    8141111        tmp2 = gauss(dest, 0., 0.);
    815         printf("%.5g: (%02d %02d %02d %02d)\n",
    816                1000. * dist(tmp, tmp2, 1.), a, b, c, d);
    817         free(dest);
    818         free(tmp2);
     1112        msg("%.5g: (%02d %02d %02d %02d)\n",
     1113            1000. * dist(tmp, tmp2, 1.), a, b, c, d);
     1114        del(dest);
     1115        del(tmp2);
    8191116    }
    8201117
     
    8231120
    8241121#if 0
    825     printf("%s\n", argv[1]);
     1122    msg("%s\n", argv[1]);
    8261123
    8271124    dest = ed(src, false, 7, 0, 0, 3, 5, 1, 0, 0, 0, 0, 0, 0);
    8281125    study(src, dest, 1.2, 0.01);
    829     free(dest);
     1126    del(dest);
    8301127
    8311128    dest = ed(src, false, 7, 5, 3, 5, 7, 5, 3, 1, 3, 5, 3, 1);
    8321129    study(src, dest, 1.2, 0.01);
    833     free(dest);
     1130    del(dest);
    8341131
    8351132    dest = ostromoukhov(src);
    8361133    study(src, dest, 1.2, 0.01);
    837     free(dest);
     1134    del(dest);
    8381135
    8391136    dest = ed(src, true, 2911, 0, 0, 1373, 3457, 2258, 0, 0, 0, 0, 0, 0);
    8401137    study(src, dest, 1.2, 0.01);
    841     free(dest);
    842 
    843     printf("\n");
     1138    del(dest);
     1139
     1140    msg("\n");
    8441141#endif
    8451142
     
    8521149    //for(sigma = 0.8; sigma < 2.; sigma *= 1.01)
    8531150    {
    854         printf("%g ", sigma);
     1151        msg("%g ", sigma);
    8551152        study(src, dest, sigma, 0.01);
    8561153    }
     
    8651162    {
    8661163        char buf[1024];
    867         printf("%i: %g\n", a, sigma);
     1164        msg("%i: %g\n", a, sigma);
    8681165        dest = dbs(src, tmp, sigma, 0., 0.);
    8691166        sprintf(buf, "output-dbs-%i.bmp", a++);
  • research/2008-displacement/study-4.sh

    r2281 r2286  
    2727
    2828    touch "$OUTPUT"
    29     ./main -4 $i | tee -a "$OUTPUT" | while read line; do
     29    ./main -4 "$i" | tee -a "$OUTPUT" | while read line; do
    3030        echo -n .
    3131    done
  • research/2008-displacement/study-5.sh

    r2279 r2286  
    2626        continue
    2727    fi
    28     if TMP="$(./main -5 $i)"; then
     28    if TMP="$(./main -5 "$i")"; then
    2929        echo "$TMP" >> "$OUTPUT"
    3030        echo "$TMP"
  • research/2008-displacement/study-6.sh

    r2285 r2286  
    1818
    1919find $1 -type f | while read i; do
    20     OUTPUT="$(echo "$i" | tr -d / | sed 's/^[.]*//').txt"
     20    OUTPUT="$(echo "$i" | tr -d / | sed 's/^[.]*//')-raster.txt"
    2121    echo "$0: outputting to $OUTPUT"
    2222
     
    2727
    2828    touch "$OUTPUT"
    29     ./main -4 $i | tee -a "$OUTPUT" | while read line; do
     29    ./main -6 "$i" | tee -a "$OUTPUT" | while read line; do
    3030        echo -n .
    3131    done
  • research/2008-displacement/study-7.sh

    r2285 r2286  
    1818
    1919find $1 -type f | while read i; do
    20     OUTPUT="$(echo "$i" | tr -d / | sed 's/^[.]*//').txt"
     20    OUTPUT="$(echo "$i" | tr -d / | sed 's/^[.]*//')-serp.txt"
    2121    echo "$0: outputting to $OUTPUT"
    2222
     
    2727
    2828    touch "$OUTPUT"
    29     ./main -4 $i | tee -a "$OUTPUT" | while read line; do
     29    ./main -7 "$i" | tee -a "$OUTPUT" | while read line; do
    3030        echo -n .
    3131    done
  • research/2008-displacement/vote.c

    r2284 r2286  
    6666    int n, i, j;
    6767
     68    for(i = 0; i < nitems; i++)
     69        values[i] = 9999.;
     70
    6871    f = fopen(file, "r");
    6972    if(!f)
Note: See TracChangeset for help on using the changeset viewer.