/* hamheuristic.c  */

#define USAGE "hamheuristic [-sgu] [-vq] [-L#] [-t#] [infile [outfile]]"

#define HELPTEXT \
" Apply a heuristic for finding hamiltonian cycles.\n\
  Output those which are unsuccessful.\n\
\n\
    -s  force output to sparse6 format\n\
    -g  force output to graph6 format\n\
        If neither -s or -g are given, the output format is\n\
        determined by the header or, if there is none, by the\n\
        format of the first input graph.\n\
    -u  Suppress output to outfile, give statistics instead.\n\
\n\
    The output file will have a header if and only if the input file does.\n\
\n\
    -p  Be content with a hamiltonian path\n\
    -v  Give a cycle or path if one is found.\n\
    -L# Limit number of sideways steps (default 1000+5*n)\n\
    -t# Try # times (default 1)\n\
\n\
    -q  suppress auxiliary information\n"

#define DEBUG(X) X

/*************************************************************************/

#include "gtools.h"

#define NO 0
#define YES 1
#define TIMEOUT 2

static long seed = 314159265;

/**************************************************************************/

static int
hamheur(sparsegraph *sg, boolean pathok, unsigned long limit, int *cyc)
/* Try up to limit, fill cyc if YES and cyc!=NULL */
/* For pathok = TRUE, return when a hamiltonian path is found */
{
    DYNALLSTAT(int,path,path_sz);
    DYNALLSTAT(int,work,work_sz);
    DYNALLSTAT(int,pos,pos_sz);  /* position on path or -1 */
    int end0,end1,v0,v1;
    size_t *v;
    int n,*e,*d,len,d0,*e0,d1,*e1,dx,*ex;
    int i,ix,x,j,w,vext,exts;
    unsigned long count;
    boolean left,cycle;
    long ran;
 
    SG_VDE(sg,v,d,e);
    n = sg->nv;
    if (n < 3) return NO;

    DYNALLOC1(int,path,path_sz,2*n+4,"malloc");
    DYNALLOC1(int,pos,pos_sz,n,"malloc");
    DYNALLOC1(int,work,work_sz,n,"malloc");

    for (i = 0; i < n; ++i)
    {
        if (d[i] < 2) return NO;
        pos[i] = -1;
    }
    count = 0;
 
    v0 = KRAN(n);
    v1 = e[v[v0]+KRAN(d[v0])];
    end0 = n+1;
    end1 = n+2;
    path[end0] = v0; pos[v0] = end0;
    path[end1] = v1; pos[v1] = end1;
    len = 2;

    for (; count < limit;)
    {
     /* First look for an extension, but note if cycle */

        if (pathok && len == n)
        {
            if (cyc)
            {
                j = 0;
                for (i = end0; i <= end1; ++i)
                    cyc[j++] = path[i];
            }
            return YES;
        }

        v0 = path[end0];
        v1 = path[end1];

        exts = 0;
        d0 = d[v0];
        e0 = e + v[v0];
        cycle = FALSE;
        ran = NEXTRAN;
        for (i = 0; i < d0; ++i)
        {
            w = e0[i];
            if (pos[w] < 0)
            {
                ++exts;
                if (ran%exts == 0) {left = TRUE; vext = w;}
            }
            else if (w == v1) cycle = TRUE;
        }
        d1 = d[v1];
        e1 = e + v[v1];
        ran = NEXTRAN;
        for (i = 0; i < d1; ++i)
        {
            w = e1[i];
            if (pos[w] < 0) 
            {
                ++exts;
                if (ran%exts == 0) {left = FALSE; vext = w;}
            }
        }

        if (exts > 0)
        {
            if (left)
            {
                --end0;
                path[end0] = vext;
                pos[vext] = end0;
            }
            else
            {
                ++end1;
                path[end1] = vext;
                pos[vext] = end1;
            }
            ++len;
            continue;
        }

     /* Can't extend but a cycle was found */

        if (cycle)
        {
            if (len == n)
            {
                if (cyc)
                {
                    j = 0;
                    for (i = pos[0]; i <= end1; ++i) cyc[j++] = path[i];
                    for (i = end0; i < pos[0]; ++i) cyc[j++] = path[i];
                }
                return YES;
            }

            ix = end0 + KRAN(len);
            for (i = 0; i < len; ++i)      
            {
                x = path[ix];

                exts = 0;
                dx = d[x];
                ex = e + v[x];
                ran = NEXTRAN;
                for (j = 0; j < dx; ++j)
                {
                    w = ex[j];
                    if (pos[w] < 0) 
                    {
                        ++exts;
                        if (ran%exts == 0) vext = w;
                    }
                }
                if (exts > 0) break;
                if (ix == end1) ix = end0; else ++ix;
            }
            if (i == len) return NO;    /* isolated component */

            work[0] = vext;
            j = 1;
            for (i = ix; i <= end1; ++i) work[j++] = path[i];
            for (i = end0; i < ix; ++i) work[j++] = path[i];

            for (i = 0; i < j; ++i)
            {
                path[end0+i] = work[i];
                pos[work[i]] = end0+i;
            }
            ++end1;
            ++len;
            continue;
        }

     /* In the last resort, do a sideways move. */
    
        exts = 0;
        d0 = d[v0];
        e0 = e + v[v0];
        ran = NEXTRAN;
        for (i = 0; i < d0; ++i)
        {
            w = e0[i];
            if (pos[w] >= 0 && w != path[end0+1])
            {
                ++exts;
                if (ran%exts == 0) {left = TRUE; vext = w;}
            }
        }
        d1 = d[v1];
        e1 = e + v[v1];
        ran = NEXTRAN;
        for (i = 0; i < d1; ++i)
        {
            w = e1[i];
            if (pos[w] >= 0 && w != path[end1-1]) 
            {
                ++exts;
                if (ran%exts == 0) {left = FALSE; vext = w;}
            }
        }

        if (left)
        {
            i = end0;
            j = pos[vext]-1;
        }
        else
        {
            i = pos[vext]+1;
            j = end1;
        }
        for (; i < j; ++i, --j)
        {
            w = path[i];
            path[i] = path[j];
            path[j] = w;
            pos[path[i]] = i;
            pos[path[j]] = j;
        }

        ++count;
    }

    return TIMEOUT;
}

/**************************************************************************/

int
main(int argc, char *argv[])
{
    sparsegraph sg;
    int n,codetype;
    int argnum,i,j,outcode,tvalue;
    char *arg,sw;
    boolean badargs;
    boolean pswitch,sswitch,gswitch,qswitch,Lswitch,tswitch,vswitch,uswitch;
    long Lvalue;
    double t;
    char *infilename,*outfilename;
    FILE *infile,*outfile;
    nauty_counter nin,nout,nNO,nYES,nTIMEOUT;
    int status;
    DYNALLSTAT(int,cyc,cyc_sz);

    HELP; PUTVERSION;

    INITRANBYTIME;

    uswitch = sswitch = gswitch = Lswitch = qswitch = FALSE;
    pswitch = vswitch = tswitch = FALSE;
    infilename = outfilename = NULL;

    argnum = 0;
    badargs = FALSE;
    for (j = 1; !badargs && j < argc; ++j)
    {
        arg = argv[j];
        if (arg[0] == '-' && arg[1] != '\0')
        {
            ++arg;
            while (*arg != '\0')
            {
                sw = *arg++;
                     SWBOOLEAN('s',sswitch)
                else SWBOOLEAN('g',gswitch)
                else SWBOOLEAN('u',uswitch)
                else SWBOOLEAN('p',pswitch)
                else SWBOOLEAN('q',qswitch)
                else SWBOOLEAN('v',vswitch)
                else SWLONG('L',Lswitch,Lvalue,"-L")
                else SWINT('t',tswitch,tvalue,"-t")
                else badargs = TRUE;
            }
        }
        else
        {
            ++argnum;
            if      (argnum == 1) infilename = arg;
            else if (argnum == 2) outfilename = arg;
            else                  badargs = TRUE;
        }
    }

    if ((sswitch != 0) + (gswitch != 0) + (uswitch != 0) > 1)
        gt_abort(">E hamheuristic: -s, -g and -u are incompatible\n");

    if (!tswitch) tvalue = 1;

    if (badargs || argnum > 2)
    {
        fprintf(stderr,">E Usage: %s\n",USAGE);
        GETHELP;
        exit(1);
    }

    if (!qswitch)
    {
        fprintf(stderr,">A hamheuristic");
        if (pswitch || sswitch || gswitch || vswitch || tswitch || Lswitch)
            fprintf(stderr," -");
        if (sswitch) fprintf(stderr,"s");
        if (gswitch) fprintf(stderr,"g");
        if (uswitch) fprintf(stderr,"u");
        if (vswitch) fprintf(stderr,"v");
        if (pswitch) fprintf(stderr,"p");
        if (Lswitch) fprintf(stderr,"L%ld",Lvalue);
        if (tswitch) fprintf(stderr,"t%d",tvalue);
        if (argnum > 0) fprintf(stderr," %s",infilename);
        if (argnum > 1) fprintf(stderr," %s",outfilename);
        fprintf(stderr,"\n");
        fflush(stderr);
    }

    if (infilename && infilename[0] == '-') infilename = NULL;
    infile = opengraphfile(infilename,&codetype,FALSE,1);
    if (!infile) exit(1);
    if (!infilename) infilename = "stdin";

    NODIGRAPHSYET(codetype);

    if (!uswitch)
    {
        if (!outfilename || outfilename[0] == '-')
        {
            outfilename = "stdout";
            outfile = stdout;
        }
        else if ((outfile = fopen(outfilename,"w")) == NULL)
        {
            fprintf(stderr,"Can't open output file %s\n",outfilename);
            gt_abort(NULL);
        }

        if (sswitch || (!gswitch && (codetype&SPARSE6)))
            outcode = SPARSE6;
        else 
            outcode = GRAPH6;

        if (codetype&HAS_HEADER)
        {
            if (outcode == SPARSE6) writeline(outfile,SPARSE6_HEADER);
            else                    writeline(outfile,GRAPH6_HEADER);
        }
    }

    nin = nout = 0;
    nNO = nYES = nTIMEOUT = 0;
    t = CPUTIME;
    SG_INIT(sg);

    while (read_sg(infile,&sg) != NULL)
    {
        ++nin;

        n = sg.nv;

        if (vswitch) { DYNALLOC1(int,cyc,cyc_sz,n,"malloc"); }
        else         cyc = NULL;

        status = TIMEOUT;
        for (i = 0; i < tvalue; ++i)
        {
            status = hamheur(&sg,pswitch,(Lswitch?Lvalue:1000+5*n),cyc);
            if (status != TIMEOUT) break;
        }

        if (status == NO) ++nNO;
        else if (status == YES) ++nYES;
        else ++nTIMEOUT;

        if (status != YES && !uswitch)
        {
            if (outcode == SPARSE6) writes6_sg(outfile,&sg);
            else if (outcode == GRAPH6) writeg6_sg(outfile,&sg);
            ++nout;
        }

        if (vswitch)
        {
            fprintf(stderr,">H " COUNTER_FMT ":",nin);
            if (status == TIMEOUT)
                fprintf(stderr," timed out\n");
            else if (status == NO)
                    fprintf(stderr," disconnected or low degree\n");
            else
            {
                for (i = 0; i < n; ++i) fprintf(stderr," %d",cyc[i]);
                fprintf(stderr,"\n");
            }
        }
    }
    t = CPUTIME - t;

    if (uswitch)
    {
        fprintf(stderr,">Z " COUNTER_FMT " graphs read from %s; "
           COUNTER_FMT " hamiltonian, " COUNTER_FMT " not, "
           COUNTER_FMT " timed out; %3.2f sec.\n",
           nin,infilename,nYES,nNO,nTIMEOUT,t);
    }
    else if (!qswitch)
    {
        fprintf(stderr,
                ">Z " COUNTER_FMT " graphs read from %s, "
                COUNTER_FMT " written to %s; %3.2f sec.\n",
                nin,infilename,nout,outfilename,t);
    }

    exit(0);
}