#include <bits/stdc++.h>
using namespace std;

struct pardes{ // parallel design
        int q,k,r,cost; // q is the number of blocks per class, k is the block size and r is the number of classes
        vector < vector<int> > classes; // each vector stores the "permutation" for a paralel class
        vector < vector<int> > a; // matrix of satisfied pairs
        pardes(int qp,int kp,int rp){ // this just sets up the size of the "classes" vectors
                // in which each class is the same "continuous" blocks
                q = qp;
                k = kp;
                r = rp;
                vector <int> v (q*k);
                for(int i=0;i<q*k;i++){
                        v[i] = i;
                }
                for(int i=0;i<r;i++){
                        classes.push_back(v);
                }
        }
        void builda(){ // builds the a matrix from scratch
                a = vector < vector <int> >  (q*k, vector <int> (q*k) ); // matrix of fullfiled pairs;
                for(int i=0;i<r;i++){
                        for(int j=0;j<q;j++){
                                for(int p1=0;p1<k;p1++){
                                        for(int p2=0;p2<k;p2++){
                                                if( p1 == p2) continue;
                                                a[classes[i][j*k+p1]][classes[i][j*k+p2]] ++;
                                        }
                                }
                        }
                }

        }
        void setcost(){
                // the number of uncovered pairs, as described in the article
                // this assumes a has been built and is correct
                cost = 0;
                for(int i=0;i<k*q;i++){
                        for(int j=i+1;j<k*q;j++){
                                if( a[i][j] == 0) cost ++;
                        }
                }
        }
        void  modify( int i,int j,int d){
                int old = a[i][j];
                a[i][j] += d;
                if( i<j && old == 0 && a[i][j] != 0 ) cost--;
                if( i<j && old != 0 && a[i][j] == 0) cost ++;
        }
        void transpose(int i,int p1,int p2){
                // this cheaply transposes two positions while keeping all the stuff correct
                int b1 = p1/k;
                int b2 = p2/k;
                int v1 = classes[i][p1];
                int v2 = classes[i][p2];
                p1 = p1%k;
                p2 = p2%k;
                for(int p=0;p<k;p++){
                        if( p == p1) continue;
                        int v = classes[i][b1*k+p];
                        modify(v1,v,-1);
                        modify(v,v1,-1);
                        modify(v2,v,1);
                        modify(v,v2,1);
                }
                for(int p=0;p<k;p++){
                        if( p == p2) continue;
                        int v = classes[i][b2*k+p];
                        modify(v2,v,-1);
                        modify(v,v2,-1);
                        modify(v1,v,1);
                        modify(v,v1,1);
                }
                swap(classes[i][b1*k+p1],classes[i][b2*k+p2]);
        }
        void impri(){
                int comps = 0;
                vector < vector<int> > e(q*k);
                vector < vector<int> > cs(q*k);
                vector <int> c(q*k);
                for(int i=0;i<q*k;i++){
                        for(int j=0;j<q*k;j++){
                                if(a[i][j] == 0) e[i].push_back(j);
                        }
                }
                for(int i=0;i<q*k;i++){
                        if( c[i] == 0){
                                comps++;
                                deque <int> B;
                                B.push_back(i);
                                cs[comps].push_back(i);
                                c[i] = comps;
                                while( !B.empty() ){
                                        int x = B.front();
                                        B.pop_front();
                                        for(int y: e[x]){
                                                if( c[y] ) continue;
                                                c[y] = comps;
                                                cs[comps].push_back(y);
                                                B.push_back(y);
                                        }
                                }
                        }
                }
                for(int i=1;i<=comps;i++){
                        if( cs[i].size() > 5) return;
                }
                cout << cost << "\n";
                for(int i=0;i<r;i++){
                        cout << "Class " << i << ": ";
                        for(int j=0;j<q;j++){
                                for(int p=0;p<k;p++){
                                        cout <<  classes[i][j*k+p] << " ";
                                }
                                cout  << "| ";
                        }
                        cout << "\n";
                }
                cout << "Components formed by missing edges\n";
                for(int i=1;i<=comps;i++){
                        if( cs[i].size() == 1) continue;
                        cout << "component : ";
                        for(int x: cs[i]){
                                cout << x << " ";
                        }
                        cout << "\n";
                }
        }

};

double randomdouble(){
        return rand()/ (RAND_MAX +1.);
}

int flip(int delta,int i){ // does a coin flip with probability e^(delta/c_i)
        double ci = 0.5* pow(0.99,pow(10,-4)*i);
        double p = exp(delta/ci);
        if( randomdouble() < p ) return 1;
        return 0;
}



int nextsiman( pardes& D, int i){ // this function gives us the next term in the search
        int q = D.q, k = D.k, r = D.r;
        // first select a random class
        int x = rand()%r;
        // next select two points in different blocks
        int p1,p2;
        while( true ){
                p1 = rand()%(k*q);
                p2 = rand()%(k*q);
                if( p1/k != p2/k) break;
        }
        int oldcost = D.cost;
        // now consider the design in which those blocks are swapped
        D.transpose(x,p1,p2);
        int newcost =  D.cost;
        double delta = (oldcost - newcost);
        if( delta >= 0 ) return 1;
        if( flip(delta,i) ) return 1;
        else D.transpose(x,p1,p2);
        return 0;
}

pardes tryfind(int q, int k, int r){
        // this method tries to find a covering design resolvable in r classes
        // utilizing the simmulated annealing method discussed in the article
        vector < pardes > starters; // these are the designs from which the system starts out
        int T = 1;
        // our first approach starts out with a set of T random coverings
        for(int i=0;i<T;i++){ // we are creating our starting designs
                pardes D(q,k,r);
                for(int j=0;j<r;j++){
                        random_shuffle(D.classes[j].begin(),D.classes[j].end()); // randomly shuffle each class
                }
                D.builda();
                D.setcost();
                starters.push_back(D);
        }

        int depth = 800000; // how deep we are willing to go in the search

        pardes Best = starters[0];

        vector < pardes> temp;
        for(  pardes D : starters ){
                for(int i=0;i< depth;){
                        i += nextsiman(D,i);
                        if( D.cost < Best.cost ){
                                Best = D;
                        }
                        if( D.cost < 10 ){
                                D.impri();
                        }
                }
                temp.push_back(D);
        }
        starters = temp;
        return Best;

}

int main(){
        srand(time(0));
        int k = 5;
        int q = 6;
        int r ;
        cout << "Insert r\n";
        cin >> r;
        pardes D = tryfind(q,k,r);
        cout << D.cost << endl;
        D.impri();
}