//important parameters that slow it down: the number of rows (significantly) and the plethysm coefficient
#include "straighten.h"
#include "plethysm.h"
#include <getopt.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <dirent.h>
#include <time.h>

typedef long long int i64;
#define FINITEFIELDSIZE 1000003LL
//define FINITEFIELDSIZE 1009LL
typedef enum { false, true } bool;
#define NUMTHREADS 2
#define MAXCOLUMNLENGTH 10
#define MAXBOXES 75
#define MAXNUMBEROFLETTERS 15
//#define MAXPREPROCESSEDDETERMINANTSPERPOINT 32768 //2^15
#define MAXPREPROCESSEDDETERMINANTSPERPOINT 1024 //2^10
#define MAXPLETHYSM 100
//higher MAXPLETHYSM will take a very long time to fill the MAXPLETHYSM x MAXPLETHYSM matrix with evaluations
#define MAXOUTER 10
//higher MAXOUTER is slow, because there will be partitions with that many rows

//<partitions and tableaux>
typedef struct partition_t { int x[MAXBOXES]; } partition_t;

void debug_printpartition(partition_t p) {
  int i;
  for (i=0;i<MAXBOXES;i++) {
    if (p.x[i]==0) {
      printf("\n");
      return;
    }
    if (i!=0) printf(",");
    printf("%d",p.x[i]);
  }
}

typedef struct {
  partition_t shape;
  partition_t transposed_shape;
  int data_col_row[MAXBOXES][MAXCOLUMNLENGTH];
} tableau_t;

int numberofboxes(partition_t shape) {
  int ret = 0;
  {int row;for (row=0;row<MAXCOLUMNLENGTH;row++) {
    ret+=shape.x[row];
  }}
  if (ret>MAXBOXES) {
    printf("error: numberofboxes>MAXBOXES\n");
    exit(0);
  }
  return ret;
}

void debug_printtableau(tableau_t *t) {
  {int row;for (row=0;row<t->transposed_shape.x[0];row++) {
    {int col;for (col=0;col<t->shape.x[row];col++) {
      printf("%3i",(t->data_col_row[col][row])+1);
    }}
    printf("\n");
  }}
}


void init_randomizetableau(struct sstd_data_c * sstd_data, struct shape_data_c * s_data, partition_t p, int outer, int inner) {
    struct sstd_data_c_options options = {};

    uint8_t max_filling_value = outer;
    uint8_t * content = (uint8_t*) calloc(max_filling_value, sizeof(uint8_t));
    int c;
    for(c = 0; c < max_filling_value; c++) {
      content[c] = inner;
    }

    uint32_t shape_length = 0;
    for(c = 0; c < MAXBOXES; c++) {
      if(p.x[c] != 0) {
        shape_length++;
      }
    }
    uint32_t * shape = (uint32_t*) calloc(shape_length, sizeof(uint32_t));
    for(c = 0; c < shape_length; c++) {
      shape[c] = p.x[c];
    }
    //uint32_t shape[5] = {18,7,6,4,1};  
    
    //uint8_t max_filling_value = 6;//max outer
    //uint8_t content[6] = {6,6,6,6,6,6};

    // set the options to just create the semistandard tableau
    options.set_sstd = 1;
    //options.gen_tableau_threaded = 1;
    // if we are straightening then the semistandard tableau array must be sorted in the row word order
    // however, if we are only interested in a random tableau then they do not need to be sorted
    options.do_not_sort_sstd = 1;

    // this constructs the shape_data_c struct that is required for all operations
    construct_shape_data_c(s_data, shape, shape_length);

    // this constructs the sstd_data_c struct that contains the array of semistandard tableau and the D-basis
    construct_sstd_data_c(sstd_data, content, max_filling_value, s_data, &options);

    free(shape);
    free(content);
}

void finalize_randomizetableau(struct sstd_data_c * sstd_data, struct shape_data_c * s_data) {
    // free the data
    destruct_shape_data_c(s_data);
    destruct_sstd_data_c(sstd_data);
}

void randomizetableau(tableau_t *tp, int inner, struct sstd_data_c * sstd_data, struct shape_data_c * s_data) {
  struct tableau *rt = random_sstd(sstd_data);
  /*
  int count=0;
  {int row;for (row=0;row<tp->transposed_shape.x[0];row++) {
    {int col;for (col=0;col<tp->shape.x[0];col++) {
      uint32_t* box = &s_data[0].row_col_to_box[count];
      tp->data_col_row[col][row]=t->entries[*box]-1;
      count++;
    }}
  }}
  */
  {int row;for (row=0;row<s_data[0].conjugate[0];row++) {
    {int col;for (col=0;col<s_data[0].shape[row];col++) {
      tp->data_col_row[col][row] = rt->entries[s_data[0].row_col_to_box[(row * s_data[0].shape[0]) + col]] - 1;
    }}
  }}
}


/*
void randomizetableau(tableau_t *tp, int inner) {
  int numbox = numberofboxes((*tp).shape);
  int outer = numbox/inner;
  if (outer>MAXOUTER) {
    printf("error: outer>MAXOUTER\n");
    exit(0);
  }
  
  outermost:
  {
  
    {int row;for (row=0;row<MAXCOLUMNLENGTH;row++) {
      {int col;for (col=0;col<MAXBOXES;col++) {
        (*tp).data_col_row[col][row]=-2;
      }}
    }}
    {int row;for (row=0;row<(*tp).transposed_shape.x[0];row++) {
      {int col;for (col=0;col<(*tp).shape.x[row];col++) {
        (*tp).data_col_row[col][row]=-1;
      }}
    }}
    {int out;for (out=0;out<outer;out++) {
      {int in;for (in=0;in<inner;in++) {
        bool found = false;
        int tries=0;
        while (!found) {
          tries++;
          if(tries==1000) {
            goto outermost;
          }
          int row,col;
          row = rand()%((*tp).transposed_shape.x[0]);
          col = rand()%((*tp).shape.x[0]);
          if ((*tp).data_col_row[col][row] == -1) {
            bool repeat = false;
            {int row2;for (row2=0;row2<(*tp).transposed_shape.x[col];row2++) {
              if ((*tp).data_col_row[col][row2]==out) {
                repeat = true;
              }
            }}
            if (!repeat) {
              (*tp).data_col_row[col][row] = out;
              found = true;
            }
          }
        }
      }}
    }}
  }
}
*/

partition_t transpose(partition_t partition) {
  partition_t ret;
  {int i;for (i=0;i<MAXBOXES;i++){
    ret.x[i]=0;
  }}
  {int index;for (index=0;index<MAXBOXES;index++) {
    int count=0;
    {int row;for (row=0;row<MAXBOXES;row++) {
      if (partition.x[row]>index) {
        count++;
      } else {
        break;
      }
    }}
    ret.x[index]=count;
  }}
  return ret;
}

//</partitions and tableaux>



i64 preprocesseddeterminants[MAXPLETHYSM][MAXPREPROCESSEDDETERMINANTSPERPOINT];




//==============================================================================

//euclidean algorithm implemented by Jesko H\"uttenhain
//typedef long long int i64;

typedef struct {
    i64 s;
    i64 t;
} bezout_t;

bezout_t exeucl(i64 a, i64 b) { /* invert a mod b */
    bezout_t B = { 1,0 };
    if (!b) return B;
    else B = exeucl(b, a%b);
    i64 t = B.s - (a/b)*B.t;
    B.s = B.t;
    B.t = t;
    return B;
}

i64 invert(i64 a) {
    bezout_t B = exeucl(FINITEFIELDSIZE,a);
        while (B.t<0) B.t+=FINITEFIELDSIZE;
    return B.t;
}

void debug_printmatrix(int dimension, i64 workmatrix[MAXPLETHYSM][MAXPLETHYSM]) {
        {int row;for(row=0;row<dimension;row++){
                {int col;for(col=0;col<dimension;col++){
                        printf("%lld,",workmatrix[row][col]);
                }}
                printf("\n");
        }}
        printf("\n");
}

i64 gaussalgor_det(int dimension, i64 *workmatrix[MAXPLETHYSM]) {
        if (dimension==0) return 0;
	i64 product = 1;
	{int row;for(row=0;row<dimension;row++){
		i64 e = workmatrix[row][row];
		if (e==0) {
			bool nonzeropivot = false;
			{int row2;for(row2=row+1;row2<dimension;row2++){
				if (workmatrix[row2][row]!=0) {
					{int col;for(col=row;col<dimension;col++){
						i64 buffer = workmatrix[row][col];
						workmatrix[row][col] = workmatrix[row2][col];
						workmatrix[row2][col] = buffer;
					}}
					e=workmatrix[row][row];
					nonzeropivot=true;
					product = FINITEFIELDSIZE-product;//pivot switch gives rise to a sign change.
					break;
				}
			}}
			if (!nonzeropivot) return 0;
		}
		{int row2;for(row2=row+1;row2<dimension;row2++){
			i64 scale = workmatrix[row2][row]*invert(workmatrix[row][row])%FINITEFIELDSIZE;
			scale = FINITEFIELDSIZE-scale;
			{int col;for(col=row+1;col<dimension;col++){
				workmatrix[row2][col] += ((workmatrix[row][col]*scale)%FINITEFIELDSIZE);
				workmatrix[row2][col] %= FINITEFIELDSIZE;
			}}
		}}
		product*=workmatrix[row][row];
		product%=FINITEFIELDSIZE;
	}}
	return product;
}



//https://www.geeksforgeeks.org/count-set-bits-in-an-integer/
/* Function to get no of set bits in binary 
   representation of passed binary no. */
unsigned int countSetBits(i64 n) 
{ 
    unsigned int count = 0; 
    while (n) 
    { 
      n &= (n-1LL) ; 
      count++; 
    } 
    return count;
}

//precache all determinants:
void precachedeterminants(int maxcolumnlength, int maxnumberofletters, int pointindex) {
  //i64 workmatrix[MAXPLETHYSM][MAXPLETHYSM];
  i64* workmatrix[MAXPLETHYSM];//[MAXPLETHYSM];
  for(int i=0; i<MAXPLETHYSM; i++) {
     workmatrix[i] = (i64*)malloc(MAXPLETHYSM * sizeof(i64));
  }
  srand(7*pointindex+100);//note that srand(0) and srand(1) have the same effect
  i64 v[maxnumberofletters][maxcolumnlength];
  {
    {int i,ii;for(i=0;i<maxnumberofletters;i++)for(ii=0;ii<maxcolumnlength;ii++)v[i][ii]=(((i64)rand())%FINITEFIELDSIZE);}
  }
  
  i64 bound = 2LL<<(maxnumberofletters-1);
  {int set;for (set=0LL;set<bound;set++) {
    int dim = countSetBits(set);
    if (dim>maxcolumnlength) continue;
    int letters[dim];
    int localset = set;
    {
      int found = 0;
      int index = 0;
      while (found<dim) {
        if ((localset%2)==1) {
          letters[found]=index;
          found++;
        }
        index++;
        localset >>= 1;
      }
    }
    {int row;for (row=0;row<dim;row++) {
      {int col;for (col=0;col<dim;col++) {
        workmatrix[row][col] = v[letters[col]][row];
      }}
    }}
    preprocesseddeterminants[pointindex][set]=gaussalgor_det(dim, workmatrix);
  }}
  srand(time(0));
  for(int i=0; i<MAXPLETHYSM; i++) {
     free(workmatrix[i]);
  }
}


void precachealldeterminants() {
  {int i;for (i=0;i<MAXPLETHYSM;i++) {precachedeterminants(MAXCOLUMNLENGTH,MAXOUTER,i);}}
}



#define MAXNUMBEROFHYPEREDGES MAXBOXES
#define MAXNUMBEROFVERTICES MAXNUMBEROFLETTERS
#define MAXNUMBEROFVERTICESPERHYPEREDGE MAXCOLUMNLENGTH
#define MAXNUMBEROFHYPEREDGESPERVERTEX MAXOUTER
#define MAXNUMBEROFDIFFERENTCOLORSTOBEPLACED MAXOUTER

struct hypergraph {
        int numberofvertices;
        int numberofhyperedges;
        int numberofverticesofhyperedge[MAXNUMBEROFHYPEREDGES];
        int verticesofhyperedge[MAXNUMBEROFHYPEREDGES][MAXNUMBEROFVERTICESPERHYPEREDGE];
        int numberofhyperedgesofvertex[MAXNUMBEROFVERTICES];
        int hyperedgesofvertex[MAXNUMBEROFVERTICES][MAXNUMBEROFHYPEREDGESPERVERTEX];
};

#ifndef max
  #define max(a,b) \
    ({ __typeof__ (a) _a = (a); \
        __typeof__ (b) _b = (b); \
      _a > _b ? _a : _b; })
#endif

void initializehypergraph(tableau_t *tableau, struct hypergraph * hypergraph){
        int numberofrows=tableau->transposed_shape.x[0];
        int *shape = tableau->shape.x;
        int data[numberofboxes(tableau->shape)];
        {
          int index=0;
          {int row;for (row=0;row<numberofrows;row++) {
            {int col;for (col=0;col<tableau->shape.x[row];col++) {
              data[index] = tableau->data_col_row[col][row]+1;
              index++;
            }}
          }}
        }
        
        int numberofboxesintableau = 0;
        {int row;for (row=0;row<numberofrows;row++){
                numberofboxesintableau+=shape[row];
        }}
        
        hypergraph->numberofvertices = 0;
        {int v;for (v=0;v<numberofboxesintableau;v++){
                hypergraph->numberofvertices=max(hypergraph->numberofvertices,data[v]);
        }}
        hypergraph->numberofhyperedges = shape[0];
        {int v;for (v=0;v<hypergraph->numberofvertices;v++){
                hypergraph->numberofhyperedgesofvertex[v]=0;
        }}
        {int e;for (e=0;e<hypergraph->numberofhyperedges;e++){
                hypergraph->numberofverticesofhyperedge[e]=0;
        }}

        {int row=0;
        int col=0;
        int v;for (v=0;v<numberofboxesintableau;v++){
                hypergraph->hyperedgesofvertex[data[v]-1][hypergraph->numberofhyperedgesofvertex[data[v]-1]] = col;
                hypergraph->numberofhyperedgesofvertex[data[v]-1]++;
                
                hypergraph->verticesofhyperedge[col][hypergraph->numberofverticesofhyperedge[col]] = data[v]-1;
                hypergraph->numberofverticesofhyperedge[col]++;
                
                col++;
                if (col==shape[row]) {
                        col=0;row++;
                }
        }}



//#############################
//debug output:
//#############################
/*
debug_printtableau(tableau);
     printf("#hyperedges=%d\n",hypergraph->numberofhyperedges);
     {int e;for (e=0;e<hypergraph->numberofhyperedges;e++){
             printf(" #vertices=%d\n",hypergraph->numberofverticesofhyperedge[e]);
             {int v;for (v=0;v<hypergraph->numberofverticesofhyperedge[e];v++){
                     printf("  %d\n",hypergraph->verticesofhyperedge[e][v]);
             }}
     }}
     printf("#vertices=%d\n",hypergraph->numberofvertices);
     {int v;for (v=0;v<hypergraph->numberofvertices;v++){
             printf(" #hyperedges=%d\n",hypergraph->numberofhyperedgesofvertex[v]);
             {int e;for(e=0;e<hypergraph->numberofhyperedgesofvertex[v];e++){
                     printf("  %d\n",hypergraph->hyperedgesofvertex[v][e]);
             }}
     }}
     
     exit(0);
*/
}


bool canbeplacedhereingeneral(int color, int vertex) {
  return true;
}
bool canbeplacedofhyperedge(int hyperedge, int color, i64 partialassignments[MAXNUMBEROFVERTICES], struct hypergraph * hypergraph) {
  {int i;for(i=0;i<hypergraph->numberofverticesofhyperedge[hyperedge];i++){
    if (partialassignments[hypergraph->verticesofhyperedge[hyperedge][i]]==color) {
      return false;
    }
  }}
  return true;
}

/**
 bubble sort. Speeding this up could be beneficial. This is performed very often.
**/
int sortcolumn_sign(int *column, int collength) {
        int ret=1;
        {int i;for (i=0;i<collength-1;i++) {
                {int j;for (j=1;j<collength-i;j++) {
                        if (column[j-1]>column[j]) {
                                int temp = column[j-1];
                                column[j-1]=column[j];
                                column[j]=temp;
                                ret*=-1;
                        }
                }}
        }}
        return ret;
}


i64 sumup_recursive(int numberofcolors, int pointindex, i64 partialassignments[MAXNUMBEROFVERTICES], struct hypergraph * hypergraph) {
        int best_vertex = -1;
        int best_numberofpossibilities = 100000;
        i64 best_possibilities[MAXNUMBEROFDIFFERENTCOLORSTOBEPLACED];
        {int vertex;for (vertex=0;vertex<hypergraph->numberofvertices;vertex++){
                if (partialassignments[vertex]==-1) {
                        int numberofpossibilities = 0;
                        i64 possibilities[MAXNUMBEROFDIFFERENTCOLORSTOBEPLACED];
                        {int color;for (color=0;color<numberofcolors;color++){
                                if (canbeplacedhereingeneral(color, vertex)) {
                                        bool placingok = true;
                                        {int hyperedge;for (hyperedge=0;hyperedge<hypergraph->numberofhyperedgesofvertex[vertex];hyperedge++){
                                                if (!canbeplacedofhyperedge(hypergraph->hyperedgesofvertex[vertex][hyperedge], color, partialassignments, hypergraph)) {
                                                        placingok=false;break;
                                                }
                                        }}
//                                        if (!canbeplacedbecauseofweight(color)) placingok=false;
                                        if (placingok) {
                                                possibilities[numberofpossibilities]=color;
                                                numberofpossibilities++;
                                        }
                                }
                        }}
                        if (numberofpossibilities==0) {
                                return 0;//pruning
                        }
                        if (numberofpossibilities < best_numberofpossibilities) {
                                best_numberofpossibilities=numberofpossibilities;
                                {int i;for (i=0;i<best_numberofpossibilities;i++)best_possibilities[i]=possibilities[i];}
                                best_vertex=vertex;
                        }
                }
        }}
        if (best_vertex==-1) {
                i64 product = 1;
                {int hyperedge;for(hyperedge=0;hyperedge<hypergraph->numberofhyperedges;hyperedge++){
                        {
                                int column[hypergraph->numberofverticesofhyperedge[hyperedge]];
                                {int i;for(i=0;i<hypergraph->numberofverticesofhyperedge[hyperedge];i++){
                                        column[i]=partialassignments[hypergraph->verticesofhyperedge[hyperedge][i]];
                                }}
                                int sign = sortcolumn_sign(column,hypergraph->numberofverticesofhyperedge[hyperedge]);//"column" is sorted after this
                                i64 rank = 0;
                                {int i;for (i=0;i<hypergraph->numberofverticesofhyperedge[hyperedge];i++) {
                                  rank += 1LL<<column[i];//2^(return_sortcolumn.column[i])
                                }}
//                                {int i;for (i=0;i<hypergraph->numberofverticesofhyperedge[hyperedge];i++) {
//                                  printf("%i,",column[i]);
//                                }}
//                                printf("  %lld\n",rank);
                                if (pointindex >= MAXPLETHYSM) {
                                  printf("error: pointindex >= MAXPLETHYSM\n");
                                  exit(0);
                                }
                                i64 factor = preprocesseddeterminants[pointindex][rank];
                                product*=factor;
                                product%=FINITEFIELDSIZE;
                                if (sign==-1) {
                                  product = FINITEFIELDSIZE-product;
                                }
                        }
                }}
                return product;
        }
        {
          i64 ret = 0;
          {int i;for(i=0;i<best_numberofpossibilities;i++){
                  partialassignments[best_vertex]=best_possibilities[i];
                  ret+=sumup_recursive(numberofcolors,pointindex, partialassignments, hypergraph);
                  ret %= FINITEFIELDSIZE;
          }}
          partialassignments[best_vertex]=-1;
          return ret;
        }
}

i64 sumup(int numberofcolors, int pointindex, i64 partialassignments[MAXNUMBEROFVERTICES], struct hypergraph * hypergraph) {
  {int i;for(i=0;i<MAXNUMBEROFVERTICES;i++){
          partialassignments[i]=-1;
  }}
  return sumup_recursive(numberofcolors,pointindex, partialassignments, hypergraph);
}

i64 evalTableau(tableau_t *t, int waringrank, int pointindex, i64 partialassignments[MAXNUMBEROFVERTICES], struct hypergraph * hypergraph) {
  initializehypergraph(t, hypergraph);
  return sumup(waringrank, pointindex, partialassignments, hypergraph);
}











































void tableaubasis(partition_t lambda, int inner, int plethysmcoefficient, i64 *M[MAXPLETHYSM], struct sstd_data_c * sstd_data, struct shape_data_c * s_data) {
  if (plethysmcoefficient>MAXPLETHYSM) {
    printf("error: plethysmcoefficient>MAXPLETHYSM\n");
    exit(0);
  }

  if (numberofboxes(lambda) % inner != 0) {
    printf("Error: The inner plethysm parameter does not divide the number of boxes!\n");
    exit(0);
  }

  i64 partialassignments[MAXNUMBEROFVERTICES];
  struct hypergraph hypergraph = {};
  i64* workmatrix[MAXPLETHYSM];//[MAXPLETHYSM];
  for(int i=0; i<MAXPLETHYSM; i++) {
     workmatrix[i] = (i64*)malloc(MAXPLETHYSM * sizeof(i64));
  }
  //tableau_t foundlist[MAXPLETHYSM];//foundlist[foundcount] is being tested for linear independence
  tableau_t * foundlist = (tableau_t *) calloc(MAXPLETHYSM, sizeof(tableau_t));

  int outer = numberofboxes(lambda)/inner;
  {partition_t lambda_transpose = transpose(lambda);
    int i;for (i=0;i<plethysmcoefficient;i++) {
    foundlist[i].shape=lambda;
    foundlist[i].transposed_shape=lambda_transpose;
  }}
  unsigned int foundcount = 0;
  int waringranks[MAXPLETHYSM];
  {int i;for(i=0;i<plethysmcoefficient;i++){
    waringranks[i]=foundlist[0].transposed_shape.x[0];//length(lambda);
//    waringranks[i]=outer;//XXX
  }}
  unsigned int failcount = 0;
  while (foundcount<plethysmcoefficient) {
    //a rectangular (foundcount+1)x(foundcount) matrix is stored in M
    randomizetableau(&(foundlist[foundcount]), inner, sstd_data, s_data);
    {int row;for (row=0;row<foundcount+1;row++) {
      M[row][foundcount] = evalTableau(&foundlist[foundcount],waringranks[row],row, partialassignments, &hypergraph);
    }}
    {int row;for (row=0;row<foundcount+1;row++) {
      {int col;for (col=0;col<foundcount+1;col++) {
        workmatrix[row][col]=M[row][col];
      }}
    }}
    i64 determ = gaussalgor_det(foundcount+1, workmatrix);
    if (determ!=0) {
      //printf("Tableau %d",(foundcount+1));printf(" out of %d:\n",plethysmcoefficient);
      //debug_printtableau(&(foundlist[foundcount]));
      foundcount++;//this makes the last entry in foundcount valid
      if (foundcount==plethysmcoefficient) {
        // store the results
        struct tableau * basis = (struct tableau *) calloc(plethysmcoefficient, sizeof(struct tableau));
        set_tableau_bulk(basis, plethysmcoefficient, s_data);
        int tab;
        for(tab=0;tab < plethysmcoefficient; tab++) {
          {int row;for (row=0;row<s_data[0].conjugate[0];row++) {
            {int col;for (col=0;col<s_data[0].shape[row];col++) {
              basis[tab].entries[s_data[0].row_col_to_box[(row * s_data[0].shape[0]) + col]] = foundlist[tab].data_col_row[col][row] + 1;
            }}
          }}
        }
        store_isotypic_full(outer, inner, s_data, plethysmcoefficient, basis);
        
        free(basis[0].entries);
        free(basis);
        for(int i=0; i<MAXPLETHYSM; i++) {
           free(workmatrix[i]);
        }
        free(foundlist);
        return;
      }
      {int col;for (col=0;col<foundcount;col++) {
        M[foundcount][col] = evalTableau(&foundlist[col],waringranks[foundcount],foundcount, partialassignments, &hypergraph);
      }}
      //again, a rectangular (foundcount+1)x(foundcount) matrix is stored in M
      failcount = 0;
    } else {
      failcount++;
      if (failcount>50) {
//         if (waringranks[plethysmcoefficient-1]<outer) {
//           printf("Waring rank increases!\n");
//         }
        failcount = 0;
        {int i;for (i=foundcount;i<plethysmcoefficient;i++) {
          waringranks[i]++;
          if (waringranks[i]>outer) {
            waringranks[i]=outer;
          }
        }}
        {int col;for (col=0;col<foundcount;col++) {
          M[foundcount][col] = evalTableau(&foundlist[col],waringranks[foundcount],foundcount, partialassignments, &hypergraph);
        }}
      }
    }
  }
}

int isotypic_already_exists(uint32_t outer, uint32_t inner, struct shape_data_c * s_data) {
  char filename[128];
  int index = 0;
  index += sprintf(filename, "HWV_Isotypic/PlethysmDB%u/Inner%u/[", outer, inner);
  for(int entry = 0; entry < s_data[0].conjugate[0]; entry++) {
    index += sprintf(&filename[index], "%u,", s_data[0].shape[entry]);
  }
  index--; //removes the trailing comma
  index += sprintf(&filename[index], "].txt");
  filename[index] = '\0';

  FILE *in_file = fopen(filename, "r");
  if(in_file == NULL) {
    return -1;
  }
  else {
    fclose(in_file);
    return 1;
  }
}


int main(int argc, char *argv[]) {
  srand(time(0));
  setbuf(stdout, NULL);

  //int inner = 6;
  //int outer = 8;
  //p.x[0]=31;p.x[1]=9;p.x[2]=7;p.x[3]=1;
  //int plethysmcoefficient = 178+10;

  //int inner = 6;
  //int outer = 6;
  //p.x[0]=18;p.x[1]=7;p.x[2]=6;p.x[3]=4;p.x[4]=1;//coeff = 39
  //int plethysmcoefficient = 39;//add 10 so that we see that it does not find more linearly independent ones than it should

  //int inner = 6;
  //int outer = 10;
  //p.x[0]=19;p.x[1]=14;p.x[2]=8;p.x[3]=3;p.x[4]=3;p.x[5]=3;p.x[6]=3;p.x[7]=3;p.x[8]=3;p.x[9]=1;//coeff = 1
  //int plethysmcoefficient = 2;

//  int inner = 4;
//  p.x[0]=15;p.x[1]=10;p.x[2]=7;//coeff is 14
//  int plethysmcoefficient = 15;//if this is set to 15, then it will only find 14 linearly independent tableaux

  
  // this will be a 2d array of shape - num_partitions x outer
  uint32_t * partitions;
  // this will store the plethysm coefficients associated with the partitions in the previous arry
  uint32_t * plethysms;
  uint32_t outer = atoi(argv[1]);
  uint32_t inner = atoi(argv[2]);

  if(outer < 2 || inner < 2 || outer > 25 || inner > 25) {
    straighten_log(STRAIGHTEN_FATAL, "The outer and inner plethysm values must be in the range of 2-25.");
    exit(EXIT_FAILURE);
  }
  uint32_t num_partitions = load_foulkes(outer, inner, &partitions, &plethysms);
  int plethysmcoefficient=0;
  straighten_set_log_level(1);
  precachealldeterminants();
  straighten_log(STRAIGHTEN_INFO, "Precaching done.");
  
  
  // start up the threads
  #pragma omp parallel num_threads(NUMTHREADS) private(plethysmcoefficient)
    {
      // these are the variables that used to be globals, they are now declared for each thread and passed around to subfunctions as need
      i64 *M[MAXPLETHYSM];//[MAXPLETHYSM];
      for(int i=0; i<MAXPLETHYSM; i++) {
         M[i] = (i64*)malloc(MAXPLETHYSM * sizeof(i64));
      } 
      struct shape_data_c s_data = {};
      struct sstd_data_c sstd_data = {};
      partition_t p;{int i;for (i=0;i<MAXBOXES;i++) p.x[i]=0;}

      // this for loop is sent off to the various threads
      #pragma omp for schedule(dynamic)
      for(int par = 0; par < num_partitions; par++) {
        for(int i=0; i<MAXPLETHYSM; i++) {
          memset(M[i], 0, MAXPLETHYSM * sizeof(i64));
        }
        for(int entry = 0; entry < outer; entry++) {
          p.x[entry] = partitions[(par * outer) + entry];
        }
        plethysmcoefficient = plethysms[par];

        init_randomizetableau(&sstd_data, &s_data, p, outer, inner);
        if(isotypic_already_exists(outer, inner, &s_data) == -1) {
          if(plethysmcoefficient > 0) {
            straighten_log(STRAIGHTEN_WARNING, "Starting partition number %d. (dimension:%d)", par, plethysmcoefficient);

            tableaubasis(p,inner,plethysmcoefficient, M, &sstd_data, &s_data);
      
            straighten_log(STRAIGHTEN_INFO, "Success for partition number %d.", par);
          }
          else {            
            store_isotypic_full(outer, inner, &s_data, plethysmcoefficient, NULL);            
          }
        }
        finalize_randomizetableau(&sstd_data, &s_data);
        
    }
  }

  straighten_log(STRAIGHTEN_INFO, "Successfully freed memory and exiting.");
}
