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

  Copyright
  Alessandro MIRONE
  mirone@esrf.fr

  Copyright 2002  by European Synchrotron Radiation Facility, Grenoble, 
                  France

                               ----------
 
                           All Rights Reserved
 
                               ----------

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that copyright notice and this permission notice appear in
supporting documentation, and that the names of European Synchrotron
Radiation Facility or ESRF or SCISOFT not be used in advertising or 
publicity pertaining to distribution of the software without specific, 
written prior permission.

EUROPEAN SYNCHROTRON RADIATION FACILITY DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL EUROPEAN SYNCHROTRON
RADIATION FACILITY OR ESRF BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, 
DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

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

 
#include<stdio.h>
#include<stdlib.h>
#include<iso646.h>
#include<math.h>
#include<string.h>

#include<iostream>

#include"Hilbertxx.hh"

#define FORMAT_BINARY
#define min(a,b) ((a)<(b))?(a):(b)


#define DEBUG(a)


void  Hxx_TypeForBits::print_values() {
  int n = sizeof(Hxx_TypeForBits)*8;
  int k;
  unsigned long   bits[n];

  Hxx_TypeForBits x = *this ;
  for(int k=0; k<n; k++) {
    bits[k] = x.lastBit() ; 
    x = x >> 1;
  }
  for(int k=0; k<n; k++) {
    printf("%u",bits[n-k-1]);
  }
  printf("\n");

}

void Hxx_TypeForBits::return_string(char* value) {
  int n = sizeof(Hxx_TypeForBits)*8;
  int k;
  unsigned long   bits[n];

  Hxx_TypeForBits x = *this ;
  for(int k=0; k<n; k++) {
    bits[k] = x.lastBit() ; 
    x = x >> 1;
  }
  for(int k=0; k<n; k++) {
    value[k] = bits[n-k-1] + '0';
  }
  value[n]=0;

}

int Hxx_TypeForBits::countOne() {
  int N_one;
  int n = sizeof(Hxx_TypeForBits)*8;

  N_one=0;
  Hxx_TypeForBits x = *this ;
  //for(int k=0; k<n; k++) {
  while (! x.isZero()) {
    if (x.firstBit()) N_one+=1 ;
    x = x << 1;
  }
  return N_one;
}

void Hxx_Filter::addCondition (Hxx_TypeForBits &mask, char logic_op, double factor) {
  int i;
  i=nop;
  filtreOperation[i].mask = mask;
  filtreOperation[i].operation = logic_op;
  filtreOperation[i].factor = factor;
  nop=++i;
};

Hxx_Filter::Hxx_Filter (int nop_max) {
  filtreOperation = new struct opFilter[nop_max];
  nop=0;
  max=0;
  min=0;
}

Hxx_Filter::~Hxx_Filter () {
  delete [] filtreOperation;
}

int Hxx_Filter::filtre(Hxx_TypeForBits &state) {
  int i;
  double res ;
  Hxx_TypeForBits tmp;
  res=0;


  if (nop==0) return 1;
  for(i=0;i<nop;i++){
    switch (filtreOperation[i].operation){
    case '|':
      tmp = filtreOperation[i].mask | state;
      break;
    case '&':
      tmp = filtreOperation[i].mask & state;
      break;
    case '^':
      tmp = filtreOperation[i].mask ^ state;
      break;
    default :
      printf("wrong operator in filtre");
      exit(1);
    }
    res += tmp.countOne() *filtreOperation[i].factor;
  }

  return this->validate(res);
}

int Hxx_Filter::validate(double res){
  if (res <= max and res >= min) {
    return 1;
  } else {
    return 0;
  }
}


Hxx_Chunker::Hxx_Chunker(int stride, int chunck_size ) {
  this->stride       = stride;
  this->chunck_size  = chunck_size ;

  this->chunck_pos= chunck_size ; // to trigger allocation at first use

  this->actual_chunk = this->top_chunk = new Hxx_Chunck ; 
  this->top_chunk->next=NULL;
  this->actual_chunk->ptr = NULL;
};

void Hxx_Chunker::free_space() {
  Hxx_Chunck *ptr, *next;
  ptr = this->top_chunk;
  while(ptr!=NULL) {
    next = ptr->next;
    if(ptr->ptr ) delete [] ptr->ptr ; 
    delete ptr;
    ptr=next;
  }

  this->chunck_pos= chunck_size ; // to trigger allocation at first use
  this->actual_chunk = this->top_chunk = new Hxx_Chunck ; 
  this->top_chunk->next=NULL;
  this->actual_chunk->ptr = NULL;

}

// ----------------------------------------------------------------------

void Hxx_Tree::Free_Chunkes() {
  this->chunker->free_space();
  this->left = NULL ; 
  this->right = NULL ; 
  this->previous = NULL ; 
  this->Leaf = NULL;
  this->n_items=0;

 
}



Hxx_Tree::Hxx_Tree(Hxx_Generic_Comparator * Less , Hxx_Object * user_data ,int  chunck_size  ) {


  this->Less    = Less    ; 

  this->left = NULL ; 
  this->right = NULL ; 
  this->previous = NULL ; 
  this->Leaf = NULL;

  this->chunker = new Hxx_Chunker(sizeof(Hxx_Tree_Node), chunck_size );

  this->user_data = user_data ; 
  
  this->n_items=0;

}


Hxx_Object *  Hxx_Tree::Add_Leaf(Hxx_Object * leaf, int &depth){
  Hxx_Tree_Node * ptr, *newLeaf;
  int left_right;
  ptr = this;

  depth=0;

  if( this->left == NULL  && this->right == NULL && this->Leaf == NULL) {
    // first leaf
    this-> Leaf = leaf;
    this->n_items++;
    return NULL;
  }

  while(1) {
    depth++;
    if( this->Less(leaf , ptr->Leaf  , this->user_data   ) ) {
      if( ptr->left ) {
	ptr = ptr->left;
      } else {
	left_right = -1;
	break;
      }
    } else if( this->Less( ptr->Leaf ,  leaf , this->user_data   ) )  {
      if( ptr->right ) {
	ptr = ptr->right;
      } else {
	left_right = 1;
	break;
      }
    } else {
      return  (Hxx_Object *) ptr->Leaf;
    }
  }

  // newLeaf = new  Hxx_Tree_Node ; 
  newLeaf = (Hxx_Tree_Node *) this->chunker->get_space() ; 
  this->n_items++;


  newLeaf->previous = ptr;
  newLeaf->right=newLeaf->left  = NULL ; 
  newLeaf->Leaf = leaf ; 
  
  if ( left_right  == -1 ) {
    ptr -> left  = newLeaf ;
  } else {
    ptr -> right = newLeaf ;
  }
  return NULL;
}

Hxx_Object  * Hxx_Tree::getLowest_and_remove(int &depth) {
  this->ActualNode = this;
  depth=1;
  return this->Again_getLowest_and_remove(depth);
}


Hxx_Object  * Hxx_Tree::Again_getLowest_and_remove(int &depth) {
  Hxx_Tree_Node * ptr;
  Hxx_Object * res;

  // depth--;

  ptr = this->ActualNode;
  while(ptr->left) {
    ptr = ptr->left;
    depth++;
  }
  res = ptr->Leaf ; 
  if( ptr != this ) {
    ptr->previous->left = ptr->right  ;
    if(ptr->right) {
      ptr->right->previous = ptr->previous;
    }
    this->ActualNode = ptr->previous;
    // delete ptr;  freeing is made by freeing chunkes
  } else {
    if(this->right) {
      this->Leaf = this->right->Leaf ; 
      this->left = this->right->left;
      this->right = this->right->right;
      if(this->right) this->right->previous = this;
      if(this->left) this->left->previous = this;
    } else {
      this->Leaf = NULL;
    }
    depth=0;
  }
  return res  ; 
}

Hxx_Object  * Hxx_Tree::getLowest(int &depth) {
  Hxx_Tree_Node * ptr;
  Hxx_Object * res;

  depth=0;

  ptr = this;
  while(ptr->left) {
    ptr = ptr->left;
    depth++;
  }
  res = ptr->Leaf ; 
  return res  ; 
}
// ---------------------------------------------------------------------------------------------------------------------------------


Hxx_Normal_Operators_Collection::Hxx_Normal_Operators_Collection ( int Ncreation, int Ndestruction,
								   int chunck_size) {

  this->Ncreation = this->user_data[0]=Ncreation ;
  this->Ndestruction = this->user_data[1]=Ndestruction;

  this->n_items=0;
  this->op_order = Ncreation+Ndestruction;

  this->orderer =    new Hxx_Tree ( (Hxx_Generic_Comparator *) Hxx_Normal_Operators_Collection::less , (Hxx_Object *)  this->user_data  );
 


  this->size_of_creators = (this->Ncreation)* sizeof(int);
  this->size_of_destructors = (this->Ndestruction)* sizeof(int);
  this -> size_of_teeth =  this->size_of_creators+ this->size_of_destructors  + sizeof(Hxx_FLOAT) ; 

  this->chunker = new Hxx_Chunker(this -> size_of_teeth ,   chunck_size );

  this->average_depth  = 0;
  this->average_depth2 = 0;
}




int  Hxx_Normal_Operators_Collection ::less(int  *a, int  *b , int *user_data ) {

  int op_ord = user_data[0]+user_data[1];
  for(int j = op_ord-1 ; j>=0 ; j--) {
   if(Hxx_Utilities::SCRAMBLE_FOR_BITS[a[j]]<Hxx_Utilities::SCRAMBLE_FOR_BITS[b[j]]) {
      return 1;
    } else if (Hxx_Utilities::SCRAMBLE_FOR_BITS[a[j]]>Hxx_Utilities::SCRAMBLE_FOR_BITS[b[j]]) {
      return 0;
    }
  }
  return 0;

}


void Hxx_Normal_Operators_Collection ::add_product(int * creators, int *destructors, Hxx_FLOAT coeff ) {
  int depth;
  Hxx_Object * ptr;
  Hxx_FLOAT *pf;

  char * new_space = this->chunker->get_space();
  int *prova =(int*) new_space;

  memcpy(new_space                                                      , creators    , this->size_of_creators);
  memcpy(new_space + this->size_of_creators                             , destructors , this->size_of_destructors);
  memcpy(new_space + this->size_of_creators+this->size_of_destructors   , &coeff      , sizeof(Hxx_FLOAT) );


  ptr = this -> orderer -> Add_Leaf( (Hxx_Object * ) new_space , depth);

  if(ptr==0) {
    this->average_depth  += depth;
    this->average_depth2  += depth*depth;
    this->n_items++;
  } else {
    pf = (Hxx_FLOAT*) ( ptr +this->size_of_creators+this->size_of_destructors );
    *pf += coeff; 
  }
  
};

void Hxx_Normal_Operators_Collection ::clean_space() {
  this->chunker->free_space();
}

void Hxx_Normal_Operators_Collection::add_memorized(int *newTuple, int count) {
  for(int i=0; i< this->op_order ; i++) {
    this->indexes[ count * this->op_order   + i  ] = newTuple[i]; 
  }
  this->coeffs[ count ] = *((Hxx_FLOAT * )  (& newTuple[this->op_order]));
}

void Hxx_Normal_Operators_Collection::reorder() {
  char *newTuple;

  this->indexes      = new int [ (this->op_order) *  (this->n_items) * sizeof(int)];
  this->next_indexes = new int [ (this->op_order) *  (this->n_items) * sizeof(int)];
  this->coeffs       = new Hxx_FLOAT [   (this->n_items) * sizeof(int)];
  
  printf(" PROFONDITA MEDIA = %e   VARIANZA = %e \n", this->average_depth / this->n_items ,
	 sqrt(-this->average_depth*this->average_depth+ this->n_items* this->average_depth2)/this->n_items    );


  int depth;
  int count = 0;
  newTuple = (char*) (this->orderer->getLowest_and_remove(depth));
  if(newTuple) {
    this->add_memorized((int *) newTuple, count);
    count++;
  }
  while(newTuple ) {
    newTuple = (char *) (this->orderer->Again_getLowest_and_remove(depth));
    if(newTuple) {
      this->add_memorized((int *)newTuple, count);
      count++;
    }
  }
  this->orderer->Free_Chunkes();
  if(count != this->n_items) {
    printf(" BIG BIG BIG PROBLEMS \n");
    exit(0);
  }

  this->clean_space();

  this->set_next_indexes();

}

void Hxx_Normal_Operators_Collection::set_next_indexes() {
  int next;
  int previous;
  int nop;
  for(int index = 0; index< this->op_order; index++) {
    next = this->n_items;
    previous=-1;
    for(int j= this->n_items -1; j>=0; j--) {
      nop = this->indexes[ j * this->op_order +index] ; 
      if(  nop  != previous) {
	previous =  nop;
	next = j+1;
      } 
      this->next_indexes[ j * this->op_order +index]  = next;
      for(int done =0; done < index; done++) {
	this->next_indexes[ j * this->op_order +done]  = min(next, this->next_indexes[ j * this->op_order +done] ) ;
      }
    }
  }
};

void Hxx_Normal_Operators_Collection::print() {
  for(int j=0; j<this->n_items; j++) {
    for(int index = 0; index< this->op_order; index++) {
      printf("%d %d  ", this->indexes[ j * this->op_order +index] , 
	     this->next_indexes[ j * this->op_order +index] );
    }
    printf(" == %e  ", this->coeffs[ j   ]);
    
    printf("\n");
  }
}




Hxx_FLOAT Hxx_Normal_Operators_Collection::operate( int i, Hxx_FermionicState* orig, Hxx_FermionicState* result, int & next) {
  int nop;
  if(i>= this->n_items) {
    printf(" i too big in Hxx_Normal_Operators_Collection::operate\n");
    exit(0);
  }
  int sign=1;

  memcpy(result, orig, sizeof( Hxx_FermionicState ));

  for(int index =this->op_order-1 ; index>= 0 ; index-- ) {
    nop = this->indexes[ i * this->op_order +index] ; 
    if(index < this->Ncreation) {
      Hxx_1p_FermionicOperators::operators[nop].Creation(result, result, sign );
    } else {
      Hxx_1p_FermionicOperators::operators[nop].Destruction(result, result, sign );
    }
    if(sign==0) {
      next = this->next_indexes[ i * this->op_order +index] ; 
    return Hxx_FLOAT(0.0);
    }
  }
  next = i+1 ; 
  return sign*this->coeffs[i];
}




void  Hxx_FermionicState::print_values(int mask ) {
  int n = sizeof(Hxx_TypeForBits)*8;
  int result;
  int bits[n];
  {

    Hxx_TypeForBits x = this->occupancies;
    for(int k=0; k<n; k++) {
      bits[k] = x.lastBit() ; 
      x = x >> 1;
    }
    for(int k=0; k<n; k++) {
      printf("%d",bits[n-k-1]);
    }
    printf("\n");
  }
  if(mask) {
    printf("MASK          ");
    Hxx_TypeForBits x = this->signes;
    for(int k=0; k<n; k++) {
      bits[k] = x.lastBit() ; 
      x = x >> 1;
    }
    for(int k=0; k<n; k++) {
      printf("%d",bits[n-k-1]);
    }
    printf("\n");
  }
}

int  Hxx_FermionicState::n_particles() {
  int n = sizeof(Hxx_TypeForBits)*8;
  int result;
  result=0;
  Hxx_TypeForBits x = this->occupancies;
  for(int k=0; k<n; k++) {
    result += x.lastBit() ; 
    x = x >> 1;
  }
  return result;
}

int  Hxx_FermionicState::n_particles(int no, int *ops) {

  int n = sizeof(Hxx_TypeForBits)*8;
  int result;
  result=0;
  Hxx_TypeForBits x = this->occupancies;
  int count=no-1;
  for(int k=0; k<n; k++) {

    if ( (n-k-1) == ops[count] ) {
      result += x.lastBit() ; 
      count--;
      if(count==-1) break;
    }
    x = x >> 1;
  }

  return result;


//   int n = sizeof(Hxx_TypeForBits)*8;
//   int result;
//   result=0;
//   Hxx_TypeForBits x = this->occupancies;
//   int start;
//   int end;
//   for(int ope =0 ; ope < no; ope++) {


//     if(ope==0) { 
//       start=0;
//     } else {
//       start = ops[ope-1]; 
//     }
    
//     end = ops[ope];
//     if(start>end) {
//       printf("start>end in routine Hxx_FermionicState::n_particles\n");
//       exit(0);
//     }
    
//     for(int k=start; k<end; k++) {
//       result += x.lastBit() ; 
//       x = x >> 1;
//     }
//   }
//   return result;
}



int filter(Hxx_FermionicState* s, void *clientdata) {
  return 1;
}


void Hxx_1p_FermionicOperator_::initialize(int pos) {
  int nop = sizeof(Hxx_TypeForBits)*8;

  this-> c_d_mask.setToOne();
  this-> c_d_mask  = this-> c_d_mask << (nop-pos -1) ; 
  this-> sign_mask.setTo0() ;
  this-> sign_mask = ~ this-> sign_mask ;
  this->sign_mask = this -> sign_mask >> pos + 1 ;
}

Hxx_1p_FermionicOperator_ * Hxx_1p_FermionicOperators ::Initialise() {
  int nop = sizeof(Hxx_TypeForBits)*8;

  Hxx_1p_FermionicOperator_ * operators  = new Hxx_1p_FermionicOperator_[nop];

  for(int i=0; i<nop; i++) {
    operators[i].initialize(i);
  }
  return operators;
}

int Hxx_N_counter::N(Hxx_FermionicState* orig ) {
  return orig->n_particles(this->n, this->ops)- this->reference_n;
}

Hxx_basis::Hxx_basis(int initialdimension, int chunck_size) {
  Hxx_TypeForBits zero;
  zero.setTo0();
  this->nstates=0;
  this->basis = new Hxx_FermionicState[initialdimension] ;
  for(int kb=0; kb< initialdimension; kb++) {
    this->basis[kb].initialise(zero);
  }
  this->states_chunker = new Hxx_Chunker(sizeof(Hxx_FermionicState) ,   chunck_size );
  this->orderer = new Hxx_Tree( (Hxx_Generic_Comparator *)  Hxx_FermionicState::less   );
  
  this->average_depth=0;
  this->average_depth2=0;
  this->set_filter(filter , 0);
}  

void Hxx_basis::add_state(Hxx_FermionicState * state) {
  memcpy( this->basis + this->nstates ,  state, sizeof(Hxx_FermionicState) ) ; 
  this->nstates++ ; 
}

int  Hxx_basis::operate_andputin_Tree( Hxx_Normal_Operators_Collection * collection , int Ncolls, Hxx_Filter* Hxfiltre) {
  
  int n_operators ;
  int i_operator ; 
  int nextOperator;
  int validation;
  Hxx_FermionicState * operandum ; 
  Hxx_FLOAT  result;
  int depth;
  int * scrambler ;
  Hxx_Object * found;
  Hxx_FermionicState *newState;
  
  n_operators = collection->get_n_items() ; 
  scrambler = Hxx_Utilities::create_scrambler(this->nstates);
  
  this->Ncreated = 0;


  for(int k=0; k< this->nstates; k++) {

    operandum = this->basis + scrambler[k];    
    i_operator =0 ; 

    newState = (Hxx_FermionicState * ) this->states_chunker->get_space();      
    memcpy(newState, operandum, sizeof( Hxx_FermionicState ));

    found = this->orderer->Add_Leaf( (Hxx_Object*) newState , depth);
    if(found == 0 ) {
      this->Ncreated ++;
      newState = (Hxx_FermionicState * ) this->states_chunker->get_space();
    } else {
      ((Hxx_FermionicState*) found)->tag=operandum->tag+1;
    }
  

    if(operandum->tag<Ncolls) {
      operandum->tag++;
	
      while( i_operator < n_operators ) {
	// printf("OP %d\n",i_operator);
	result  =  collection->operate( i_operator , operandum , newState, nextOperator);
	newState->tag=0;
	// printf(" result = %e \n", result);
	if (Hxfiltre==0){
	  validation = (filtre)(newState, clientdata);
	} else {
	  validation = Hxfiltre->filtre(newState->occupancies);
	}
	if(result!=Hxx_FLOAT(0.0) && validation ) {
	  found = this->orderer->Add_Leaf( (Hxx_Object*) newState , depth);
	  if(found == 0 ) {
	    this->Ncreated ++;
	    newState = (Hxx_FermionicState * ) this->states_chunker->get_space();
	    this->average_depth += depth;
	    this->average_depth2 += depth*depth;
	  }
	}
  	i_operator = nextOperator;
      }
    }
  }
  delete [] scrambler;
  return this->Ncreated;
}

//  int  Hxx_basis::operate_andputin_Tree( Hxx_Normal_Operators_Collection * collection , int Ncolls) {
  
//    int n_operators ;
//    int i_operator ; 
//    int nextOperator;
//    Hxx_FermionicState * operandum ; 
//    Hxx_FLOAT  result;
//    int depth;
//    int * scrambler ;
//    Hxx_Object * found;
//    Hxx_FermionicState *newState;
  
//    n_operators = collection->get_n_items() ; 
//    scrambler = Hxx_Utilities::create_scrambler(this->nstates);
  
//    this->Ncreated = 0;


//    for(int k=0; k< this->nstates; k++) {

//      operandum = this->basis + scrambler[k];    
//      i_operator =0 ; 

//      newState = (Hxx_FermionicState * ) this->states_chunker->get_space();      
//      memcpy(newState, operandum, sizeof( Hxx_FermionicState ));

//      found = this->orderer->Add_Leaf( (Hxx_Object*) newState , depth);
//      if(found == 0 ) {
//        this->Ncreated ++;
//        newState = (Hxx_FermionicState * ) this->states_chunker->get_space();
//      } else {
//        ((Hxx_FermionicState*) found)->tag=operandum->tag+1;
//      }
  

//      if(operandum->tag<Ncolls) {
//        operandum->tag++;
	
//        while( i_operator < n_operators ) {
//  	// printf("OP %d\n",i_operator);
//  	result  =  collection->operate( i_operator , operandum , newState, nextOperator);
//  	newState->tag=0;
//  	// printf(" result = %e \n", result);
//  	if(result!=Hxx_FLOAT(0.0) && (filtre)(newState, clientdata) ) {
//  	  found = this->orderer->Add_Leaf( (Hxx_Object*) newState , depth);
//  	  if(found == 0 ) {
//  	    this->Ncreated ++;
//  	    newState = (Hxx_FermionicState * ) this->states_chunker->get_space();
//  	    this->average_depth += depth;
//  	    this->average_depth2 += depth*depth;
//  	  }
//  	}
//    	i_operator = nextOperator;
//        }
//      }
//    }
//    delete [] scrambler;
//    return this->Ncreated;
//  }


void Hxx_basis::set_filter(Callback func, void *cdata){
    filtre = func ;
    clientdata = cdata;
  }

void Hxx_basis::reconstruct_basis(Hxx_Filter* Hxfiltre){
  Hxx_FermionicState * popup;
  Hxx_TypeForBits zero;
  zero.setTo0();
  Hxx_FermionicState * newbasis = new Hxx_FermionicState[this->Ncreated];
  for(int kb=0; kb<this->Ncreated ; kb++) {
    newbasis[kb].initialise(zero);
  }
  
  int newnstates=0;

  for(int k=0; k<this->nstates; k++) {
    popup = this->basis + k; 
    if( Hxfiltre->filtre(popup->occupancies) !=0) { 
        memcpy(newbasis  +newnstates   ,popup , sizeof(Hxx_FermionicState));
        newnstates++;
    } 
  }
  delete [] this->basis;
  this->basis = newbasis;
  this->nstates = newnstates;
}

void Hxx_basis::reconstruct_basis(CounterForHxxBasisFilter *filter){
  Hxx_FermionicState * popup;
  Hxx_TypeForBits zero;
  zero.setTo0();
  Hxx_FermionicState * newbasis = new Hxx_FermionicState[this->Ncreated];
  for(int kb=0; kb<this->Ncreated ; kb++) {
    newbasis[kb].initialise(zero);
  }
  
  int newnstates=0;

  for(int k=0; k<this->nstates; k++) {
    popup = this->basis + k; 
    if( filter(popup) ==0) { 
        memcpy(newbasis  +newnstates   ,popup , sizeof(Hxx_FermionicState));
        newnstates++;
    } 
  }
  delete [] this->basis;
  this->basis = newbasis;
  this->nstates = newnstates;
}

void Hxx_basis::reconstruct_basis() {
  Hxx_FermionicState * popup;
  int depth;
  delete [] this->basis;
  Hxx_TypeForBits zero;
  zero.setTo0();
  this->basis = new Hxx_FermionicState[this->Ncreated];

  for(int kb=0; kb<this->Ncreated ; kb++) {
    this->basis[kb].initialise(zero);
  }

  this->nstates=this->Ncreated;

  for(int k=0; k<this->Ncreated; k++) {
    if(k==0) {
      popup = (Hxx_FermionicState * )  this->orderer->getLowest_and_remove(depth);
    } else {
      popup = (Hxx_FermionicState * )  this->orderer->Again_getLowest_and_remove(depth);
    }
    memcpy(this->basis  + k  ,popup , sizeof(Hxx_FermionicState));
  }
  this->states_chunker->free_space();
  this->orderer->Free_Chunkes();

  printf( " averaged depth is %e  ; variance is %e \n", this->average_depth/this->Ncreated,
	  sqrt(- this->average_depth*this->average_depth + this->Ncreated*this->average_depth2)/ this->Ncreated );

  this->average_depth=0;
  this->average_depth2=0;

}

int Hxx_basis::find_state(Hxx_FermionicState * state  ) {
  int h= this->nstates;
  int l=-1;
  int m;
  while( h>l+1) {
    m=(l+h)>>1;
    if(Hxx_FermionicState::less(this->basis+m, state) ) {
      l=m;
    } else {
      h=m;
    }
  }
  if( h<0 || h>=this->nstates) {
    return -1;
  }
  if( (state->get_value()) != ( this->basis+h)->get_value() ){
    return -1;
  }
  return h;
}

#undef RAND_MAX
#define RAND_MAX        2147483647


int * Hxx_Utilities::create_scrambler(int nstates)
{
  std::cout << " sono in create scrambler " << std::endl;
  int *scrambler = new int [nstates];
  for(int i=0; i<nstates; i++) scrambler[i]=i;
  for(int i=0; i<4*nstates; i++) {
    int i1;
    int i2;
    int swap;
    i1 = (int) (  random()*1.0*nstates*1.0/RAND_MAX     )  ;
    i2 = (int) (  random()*1.0*nstates*1.0/RAND_MAX     )  ;
    swap = scrambler[i1];
    scrambler[i1] = scrambler[i2] ; 
    scrambler[i2] = swap ; 
  }
  return scrambler;
}
  

Hxx_TransitionMatrix::Hxx_TransitionMatrix() {

  this->chunker = new Hxx_Chunker(sizeof(Hxx_transition), 100000 ) ;
  this->n_basis_from = 0;
  this->n_basis_to   = 0;
  this->Transitions_Orderer = new Hxx_Tree( (Hxx_Generic_Comparator *) Hxx_transition::less, NULL, 100000);


  aux_total_transitions=0;
  aux_i_from=NULL;
  aux_i_to=NULL;
  aux_coeff=NULL;


}

void Hxx_TransitionMatrix::free_space() {

  chunker->free_space();

  Transitions_Orderer->Free_Chunkes();
  if(aux_total_transitions){
    delete aux_i_from;
    delete aux_i_to;
    delete aux_coeff;
    
    aux_i_from=NULL;
    aux_i_to=NULL;
    aux_coeff=NULL;
    aux_total_transitions=0;
  }
}

Hxx_TransitionMatrix::~Hxx_TransitionMatrix(){

  this->free_space();
  delete chunker;
  delete Transitions_Orderer;
}


void Hxx_TransitionMatrix::add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
					    int n_counters, Hxx_N_counter  ** diag_pro,
					    Hxx_FLOAT  coeff) {
  
  Hxx_transition *newTrans;
  Hxx_transition *foundTrans;

  int n_Abasis;

  int * A_Scrambler;

  int i_a_state, a_state, b_state;
  
  Hxx_FermionicState * fs_a;
  
  Hxx_FLOAT amplitude;
  
  int depth;

  double average_depth=0.0;
  double average_depth2=0.0;

  this->n_basis_from =  a_basis->get_dimension () ;
  this->n_basis_to   =  this->n_basis_from  ;
   
  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  A_Scrambler = Hxx_Utilities::create_scrambler(this->n_basis_from);
  
  int total_Transitions=0;
  
  for(i_a_state =0 ; i_a_state<this->n_basis_from ; i_a_state++) {

    a_state = A_Scrambler [i_a_state];
    
    fs_a = a_basis->get_basis() + a_state;
    
    amplitude = coeff;

    for(int  i_operator =0 ; i_operator < n_counters ; i_operator++ ) {
      
      amplitude  =  amplitude *  ( (*(diag_pro + i_operator ))->N(  fs_a ) ) ;

    }
      

    if(amplitude!=Hxx_FLOAT(0.0) ) {
      
      b_state = a_state;
      
      newTrans->initialise(a_state, b_state, amplitude);
      
      foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);
      
      if(foundTrans==0 ) {
	
	average_depth += depth;
	average_depth2 += depth*depth;
	total_Transitions++;
	
	newTrans = (Hxx_transition * ) (this->chunker->get_space());
	
      } else {
	
	foundTrans -> add_coeff(amplitude);
	
      }
    }
     
  }
  printf( " Hxx_TransitionMatrix::add_contribution  for   counters \n");
  printf( " ALL TRANSITIONS (%d)HAVE BEEN CALCULATED AND STORED IN A BINARY TREE\n",  total_Transitions);
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / total_Transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(total_Transitions*average_depth2 -
			  average_depth* average_depth)/ total_Transitions );
  

  delete A_Scrambler;  


}



void    Hxx_TransitionMatrix::add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			    Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff, int complete_basis) {

  Hxx_transition *newTrans;
  Hxx_transition *foundTrans;

  int n_Abasis;

  int * A_Scrambler;

  int i_a_state, a_state, b_state;
  
  int n_operators;
  int i_operator ; 
  int nextOperator;
  
  Hxx_FermionicState * fs_a;
  Hxx_FermionicState * fs_b;
  
  Hxx_FLOAT amplitude;
  Hxx_TypeForBits zero;
  zero.setTo0();
  int depth;
  double average_depth=0.0;
  double average_depth2=0.0;

  this->n_basis_from =  a_basis->get_dimension () ;
  this->n_basis_to   =  b_basis->get_dimension () ;
   
  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  A_Scrambler = Hxx_Utilities::create_scrambler(this->n_basis_from);
  
  n_operators = Norm_Op->get_n_items() ; 
  
  fs_b = new Hxx_FermionicState(zero );

  int total_Transitions=0;
  
  for(i_a_state =0 ; i_a_state<this->n_basis_from ; i_a_state++) {

    DEBUG( printf(" i_a_state %d \n", i_a_state); )
    a_state = A_Scrambler [i_a_state];
    
    fs_a = a_basis->get_basis() + a_state;
    
    i_operator =0 ; 
    
    while( i_operator < n_operators ) {
      DEBUG( printf(" i_operator  %d  n_operators %d\n", i_operator,n_operators );)
      
      amplitude  = coeff* Norm_Op->operate( i_operator , fs_a ,fs_b , nextOperator);
      
      DEBUG( printf(" amplitude    %e nextOperator  %d\n", amplitude,nextOperator);)
      
      if(amplitude!=Hxx_FLOAT(0.0) ) {
	
	b_state = b_basis->find_state(fs_b);
	
	if(b_state==-1 ) {
	  if( !complete_basis ) {
	    i_operator = nextOperator;
	    continue;
	  } else {
	    printf("MO  e per a %d .....\n",  a_basis->find_state(fs_a) );
	    printf(" b_state==-1 in routine for transitions \n");
	    printf(" valore cercato\n");
	    fs_b->print_values();
	    printf(" valore di partenza \n");
	    (a_basis->get_basis() + a_state)->print_values();
	    printf("################################################################################\n");
	    printf("Forse devi allentare la condizione sull esistenza dei nuovi stati , argomento complete_basis di add_contribution\n");
	    for(int i=0; i< this->n_basis_from; i++) {
	      /*if((a_basis->get_basis() + i+1 )->get_value()<
		(a_basis->get_basis() + i )->get_value() )
		{
		printf(" non corretto \n");
		exit(0);
		}
	      */
	      if((a_basis->get_basis() + i )->get_value() == fs_b->get_value())
		{
		  printf(" trovato \n");
		  exit(0);
		}
	      // (a_basis->get_basis() + i )->print_values();
	    }
	    exit(0);
	  }
	}
	
	newTrans->initialise(a_state, b_state, amplitude);
	
	foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);

        DEBUG( printf(" n_items %d\n", this->Transitions_Orderer->get_n_items() );)
	
	if(foundTrans==0 ) {
	  average_depth += depth;
	  average_depth2 += depth*depth;
	  total_Transitions++;
	  DEBUG( printf(" total_Transitions  %d\n", total_Transitions  );)
	  newTrans = (Hxx_transition * ) (this->chunker->get_space());
	} else {
	  foundTrans -> add_coeff(amplitude);
	}
      }
      i_operator = nextOperator;
    }
    
  }
  
  printf( " ALL TRANSITIONS (%d)HAVE BEEN CALCULATED AND STORED IN A BINARY TREE\n",  total_Transitions);
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / total_Transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(total_Transitions*average_depth2 -
			  average_depth* average_depth)/ total_Transitions );
  

  delete fs_b;
  delete A_Scrambler;

}


void    Hxx_TransitionMatrix::add_contribution_debug(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			    Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff, int complete_basis) {

  Hxx_transition *newTrans;
  Hxx_transition *foundTrans;

  int n_Abasis;

  int * A_Scrambler;

  int i_a_state, a_state, b_state;
  
  int n_operators;
  int i_operator ; 
  int nextOperator;
  
  Hxx_FermionicState * fs_a;
  Hxx_FermionicState * fs_b;
  
  Hxx_FLOAT amplitude;
  Hxx_TypeForBits zero;
  zero.setTo0();
  int depth;
  double average_depth=0.0;
  double average_depth2=0.0;

  this->n_basis_from =  a_basis->get_dimension () ;
  this->n_basis_to   =  b_basis->get_dimension () ;
   
  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  A_Scrambler = Hxx_Utilities::create_scrambler(this->n_basis_from);
  
  n_operators = Norm_Op->get_n_items() ; 
  
  fs_b = new Hxx_FermionicState(zero );

  int total_Transitions=0;
  
  for(i_a_state =0 ; i_a_state<this->n_basis_from ; i_a_state++) {

    printf(" i_a_state %d \n", i_a_state); 
    a_state = A_Scrambler [i_a_state];
    
    fs_a = a_basis->get_basis() + a_state;
    
    i_operator =0 ; 
    printf("fs_a\n");
    fs_a->print_values();
    while( i_operator < n_operators ) {
      printf(" i_operator  %d  n_operators %d\n", i_operator,n_operators );
      
      amplitude  = coeff* Norm_Op->operate( i_operator , fs_a ,fs_b , nextOperator);
      
      printf(" amplitude    %e  i_operator %d  nextOperator  %d\n", amplitude,i_operator, nextOperator);
      fs_b->print_values();
      
      if(amplitude!=Hxx_FLOAT(0.0) ) {
	printf(" aggiungo stato \n");
	b_state = b_basis->find_state(fs_b);
	
	if(b_state==-1 && complete_basis) {
	  printf("MO  e per a %d .....\n",  a_basis->find_state(fs_a) );
	  printf(" b_state==-1 in routine for transitions \n");
	  printf(" valore cercato\n");
	  fs_b->print_values();
	  (a_basis->get_basis() + a_state)->print_values();
	  printf("################################################################################\n");
	  for(int i=0; i< this->n_basis_from-1; i++) {
	    if((a_basis->get_basis() + i+1 )->get_value()<
	       (a_basis->get_basis() + i )->get_value() )
	      {
		printf(" non corretto \n");
		exit(0);
	      }
	    if((a_basis->get_basis() + i+1 )->get_value() == fs_b->get_value())
	      {
		printf(" trovato \n");
		exit(0);
	      }
	    (a_basis->get_basis() + i )->print_values();
	  }
	  exit(0);
	}
	
	newTrans->initialise(a_state, b_state, amplitude);
	
	foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);

        printf(" n_items %d\n", this->Transitions_Orderer->get_n_items() );
	
	if(foundTrans==0 ) {
	  average_depth += depth;
	  average_depth2 += depth*depth;
	  total_Transitions++;
	  printf(" total_Transitions  %d\n", total_Transitions  );
	  newTrans = (Hxx_transition * ) (this->chunker->get_space());
	} else {
	  foundTrans -> add_coeff(amplitude);
	}
      }
      i_operator = nextOperator;
    }
    
  }
  
  printf( " ALL TRANSITIONS (%d)HAVE BEEN CALCULATED AND STORED IN A BINARY TREE\n",  total_Transitions);
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / total_Transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(total_Transitions*average_depth2 -
			  average_depth* average_depth)/ total_Transitions );
  

  delete fs_b;
  delete A_Scrambler;

}
void  Hxx_TransitionMatrix::BuildAux(){
  aux_total_transitions =    Transitions_Orderer->get_n_items();
  aux_i_from = new int       [aux_total_transitions];
  aux_i_to   = new int       [aux_total_transitions];
  aux_coeff = new Hxx_FLOAT[aux_total_transitions];

  Hxx_transition * newTrans; 
  int depth;
  for(int k=0; k<aux_total_transitions; k++) {
    if(k==0) {
      newTrans  = (Hxx_transition * )  Transitions_Orderer->getLowest_and_remove(depth);
    } else {
      newTrans = (Hxx_transition * )  Transitions_Orderer->Again_getLowest_and_remove(depth);
    }
    aux_i_from [k] = newTrans->get_from() ;
    aux_i_to   [k] = newTrans->get_to()   ;
    aux_coeff[k] = newTrans->get_coeff(); 
  }
  chunker->free_space();
  Transitions_Orderer->Free_Chunkes();
};

void Hxx_TransitionMatrix::AddFrom(Hxx_TransitionMatrix * term, double factor, int dosquare){


  if( term->aux_coeff == NULL ) {
    std::cout << " ERROR term->aux_coeff == NULL in routine void Hxx_TransitionMatrix::AddFrom(Hxx_TransitionMatrix * term, double f_coeff)\n";
    exit(0);
  }
  
  int *scrambler = Hxx_Utilities::create_scrambler(term->aux_total_transitions);

  Hxx_transition * newTrans, *foundTrans;
  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  double average_depth=0.0;
  double average_depth2=0.0;
  int i_from;
  int i_to;
  int depth;
  Hxx_FLOAT f_coeff;
  // std::cout << "il totale delle trans est "<<term->aux_total_transitions << std::endl ; 
  int i;
  int nuove=0;
  for(int iseq=0; iseq<term->aux_total_transitions; iseq++) {
    i = scrambler[iseq];
    // std::cout << i << std::endl ; 
    
    i_from=term->aux_i_from[i];
    i_to= term->aux_i_to[i]  ;
    if(dosquare==0) {
      f_coeff=term->aux_coeff[i] * factor ;
    } else {
      f_coeff=term->aux_coeff[i] * term->aux_coeff[i] *  factor ;
    }

    newTrans->initialise(i_from, i_to,f_coeff );
    foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);
    average_depth += depth;
    average_depth2 += depth*depth;

    if(foundTrans==0 ) {
      nuove++;
      newTrans = (Hxx_transition * ) (this->chunker->get_space());
    } else {
      foundTrans -> add_coeff(f_coeff);
    }
  }

  delete scrambler;


  printf( " a number of %d nuove su %d   TRANSITIONS HAVE BEEN SUMMED AND STORED IN A BINARY TREE\n",nuove,  term->aux_total_transitions);
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / term->aux_total_transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(term->aux_total_transitions*average_depth2 -
			  average_depth* average_depth)/ term->aux_total_transitions );

};

     
void  Hxx_TransitionMatrix::Read_4arrays( char * nome_file, int &Nels , int *&from, 
					  int *&to, double *&coeffs,int * & qcol) {

  
  printf( " NOW READING 4  TRANSITIONS ARRAYS  FROM FILE %s\n", nome_file);
  FILE * f_in = fopen(nome_file,"r");
  
  if(!f_in) {
    printf(" ERROR file %s not found \n", nome_file);
    exit(0);
  }
  
  fread(&Nels   , sizeof(int)      , 1    , f_in) ;
  from = new int [Nels];
  to   = new int [Nels];
  qcol = new int [Nels];
  coeffs = new double [Nels];

  fread(to   , sizeof(int)      , Nels    , f_in) ;
  fread(from , sizeof(int)      , Nels    , f_in) ;
  fread(coeffs, sizeof(Hxx_FLOAT), Nels    , f_in) ;
  fread(qcol   , sizeof(int)      , Nels    , f_in) ;
  fclose(f_in);


}


void  Hxx_TransitionMatrix::Read_arrays( char * nome_file, int &Nels , int *&from, int *&to, double *&coeffs, int binary) {

  int n_basis_from , n_basis_to;
  
  printf( " NOW READING  TRANSITIONS ARRAYS  FROM FILE %s\n", nome_file);
  FILE * f_in = fopen(nome_file,"r");
  
  if(!f_in) {
    printf(" ERROR file %s not found \n", nome_file);
    exit(0);
  }
  
  if(binary) {
    fread(&Nels   , sizeof(int)      , 1    , f_in) ;
    fread(&n_basis_from  , sizeof(int)      , 1    , f_in) ;
    fread(&n_basis_to    , sizeof(int)      , 1    , f_in) ;
  } else {
    fscanf(f_in," %d\n " , &Nels);
    fscanf(f_in," %d\n " , &n_basis_from);
    fscanf(f_in," %d\n " , &n_basis_to  );
  }
  
  printf(" adesso leggo %d transitions\n", Nels);
  
  
  from = new int [Nels];
  to = new int [Nels];
  coeffs = new double [Nels];


  int index;
  int i_from , i_to;
  double f_coeff;
  
  for(int i=0; i<Nels; i++) { 
    index = i;
    if(binary) {    
      fread(&i_to   , sizeof(int)      , 1    , f_in) ;
      fread(&i_from , sizeof(int)      , 1    , f_in) ;
      fread(&f_coeff, sizeof(Hxx_FLOAT), 1    , f_in) ;
    } else { 
      fscanf(f_in," %d %d  %le", &i_to, &i_from ,  &f_coeff );
    }
    from[index]=i_from;
    to[index] = i_to;
    coeffs[index] = f_coeff ;
  }

  
  
}



void  Hxx_TransitionMatrix::Read_add( char * nome_file, Hxx_FLOAT coeff, int binary) {

  printf( " NOW READING  TRANSITIONS FROM FILE %s\n", nome_file);
  
  FILE * f_in = fopen(nome_file,"r");
  
  if(!f_in) {
    printf(" ERROR file %s not found \n", nome_file);
    exit(0);
  }
  int total_Transitions;
  if(binary) {
    fread(&total_Transitions   , sizeof(int)      , 1    , f_in) ;
    fread(&this->n_basis_from  , sizeof(int)      , 1    , f_in) ;
    fread(&this->n_basis_to    , sizeof(int)      , 1    , f_in) ;
  } else {
    fscanf(f_in," %d\n " , &total_Transitions);
    fscanf(f_in," %d\n " , &this->n_basis_from);
    fscanf(f_in," %d\n " , &this->n_basis_to  );
  }
 
  printf(" adesso leggo %d transitions\n", total_Transitions);


  int *from = new int [total_Transitions];
  int *to = new int [total_Transitions];
  Hxx_FLOAT *amplitude = new Hxx_FLOAT [total_Transitions];
  int index;
  int i_from , i_to;
  Hxx_FLOAT f_coeff;
  int depth;


  int *scrambler = Hxx_Utilities::create_scrambler(total_Transitions);

  Hxx_transition * newTrans, *foundTrans;

  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  double average_depth=0.0;
  double average_depth2=0.0;
  
  for(int i=0; i<total_Transitions; i++) { 
    index = scrambler[i];
    if(binary) {    
      fread(&i_to   , sizeof(int)      , 1    , f_in) ;
      fread(&i_from , sizeof(int)      , 1    , f_in) ;
      fread(&f_coeff, sizeof(Hxx_FLOAT), 1    , f_in) ;

    } else { 
      fscanf(f_in," %d %d  %le", &i_to, &i_from ,  &f_coeff );
    }
    from[index]=i_from;
    to[index] = i_to;
    amplitude[index] = f_coeff*coeff;
  }

  for(int i=0; i<total_Transitions; i++) {
    i_from=from[i];
    i_to= to[i]  ;
    f_coeff=amplitude[i] ;

    newTrans->initialise(i_from, i_to,f_coeff );
    // printf(" $$$$ %d \n", this->Transitions_Orderer->get_n_items() );
    foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);

   
    if(foundTrans==0 ) {
      average_depth += depth;
      average_depth2 += depth*depth;
      newTrans = (Hxx_transition * ) (this->chunker->get_space());
    } else {
      foundTrans -> add_coeff(f_coeff);
    }
  }

  delete scrambler;
  delete from;
  delete to;
  delete amplitude;
  fclose(f_in);


  printf( " ALL TRANSITIONS HAVE BEEN READ AND STORED IN A BINARY TREE\n");
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / total_Transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(total_Transitions*average_depth2 -
			  average_depth* average_depth)/ total_Transitions );
}


void  Hxx_TransitionMatrix::Read_add_forhops( char * nome_file, Hxx_FLOAT coeff, int binary,
					      int * cc, double *hopfacts) {

  int cc1, cc2;

  printf( " NOW READING  TRANSITIONS FROM FILE %s\n", nome_file);
  
  FILE * f_in = fopen(nome_file,"r");
  
  if(!f_in) {
    printf(" ERROR file %s not found \n", nome_file);
    exit(0);
  }
  int total_Transitions;
  if(binary) {
    fread(&total_Transitions   , sizeof(int)      , 1    , f_in) ;
    fread(&this->n_basis_from  , sizeof(int)      , 1    , f_in) ;
    fread(&this->n_basis_to    , sizeof(int)      , 1    , f_in) ;
  } else {
    fscanf(f_in," %d\n " , &total_Transitions);
    fscanf(f_in," %d\n " , &this->n_basis_from);
    fscanf(f_in," %d\n " , &this->n_basis_to  );
  }
 
  printf(" adesso leggo %d transitions\n", total_Transitions);


  int *from = new int [total_Transitions];
  int *to = new int [total_Transitions];
  Hxx_FLOAT *amplitude = new Hxx_FLOAT [total_Transitions];
  int index;
  int i_from , i_to;
  Hxx_FLOAT f_coeff;
  int depth;


  int *scrambler = Hxx_Utilities::create_scrambler(total_Transitions);

  Hxx_transition * newTrans, *foundTrans;

  newTrans = (Hxx_transition * ) (this->chunker->get_space());

  double average_depth=0.0;
  double average_depth2=0.0;
  
  for(int i=0; i<total_Transitions; i++) { 
    index = scrambler[i];
    if(binary) {    
      fread(&i_to   , sizeof(int)      , 1    , f_in) ;
      fread(&i_from , sizeof(int)      , 1    , f_in) ;
      fread(&f_coeff, sizeof(Hxx_FLOAT), 1    , f_in) ;

    } else { 
      fscanf(f_in," %d %d  %le", &i_to, &i_from ,  &f_coeff );
    }
    cc1=cc[i_to];
    cc2=cc[i_from];
    if( cc2>cc1) cc1=cc2;


    from[index]=i_from;
    to[index] = i_to;
    amplitude[index] = f_coeff*coeff*hopfacts[cc1]  ;
  }

  for(int i=0; i<total_Transitions; i++) {
    i_from=from[i];
    i_to= to[i]  ;
    f_coeff=amplitude[i] ;

    newTrans->initialise(i_from, i_to,f_coeff );
    // printf(" $$$$ %d \n", this->Transitions_Orderer->get_n_items() );
    foundTrans =  (Hxx_transition *) this->Transitions_Orderer->Add_Leaf((Hxx_Object*)newTrans, depth);

   
    if(foundTrans==0 ) {
      average_depth += depth;
      average_depth2 += depth*depth;
      newTrans = (Hxx_transition * ) (this->chunker->get_space());
    } else {
      foundTrans -> add_coeff(f_coeff);
    }
  }

  delete scrambler;
  delete from;
  delete to;
  delete amplitude;
  fclose(f_in);


  printf( " ALL TRANSITIONS HAVE BEEN READ AND STORED IN A BINARY TREE\n");
  printf( " AVERAGE DEPTH IS \n");
  printf( " %e \n" , average_depth / total_Transitions );
  printf( " VARIANCE AROUND AVERAGE DEPTH IS \n");
  printf( " %e \n" , sqrt(total_Transitions*average_depth2 -
			  average_depth* average_depth)/ total_Transitions );
}




void Hxx_TransitionMatrix::create_3Array_and_Clean(int *& from, int *& to, Hxx_FLOAT *& coeff, int &ntrans) {
  
  int total_Transitions = Transitions_Orderer->get_n_items();
  
  ntrans = total_Transitions;
  
  from = new int       [total_Transitions];
  to   = new int       [total_Transitions];
  coeff = new Hxx_FLOAT[total_Transitions];

  Hxx_transition * newTrans; 
  
  int depth;

  for(int k=0; k<total_Transitions; k++) {

    if(k==0) {
      newTrans  = (Hxx_transition * )  Transitions_Orderer->getLowest_and_remove(depth);
    } else {
      newTrans = (Hxx_transition * )  Transitions_Orderer->Again_getLowest_and_remove(depth);
    }
    from [k] = newTrans->get_from() ;
    to   [k] = newTrans->get_to()   ;
    coeff[k] = newTrans->get_coeff();

  }
  this->free_space();
}  




void  Hxx_TransitionMatrix::Write_and_Clean( char * nome_file, int binary) {
   
  printf( " NOW SORTING  TRANSITIONS TO FILE %s\n", nome_file);
  
  FILE * f_out = fopen(nome_file,"w");
  
  int total_Transitions = Transitions_Orderer->get_n_items();
  
 
  Hxx_transition * newTrans; 
  
  int depth;
  if(binary==1) {
    fwrite(&total_Transitions   , sizeof(int)      , 1    , f_out) ;
    fwrite(&this->n_basis_from  , sizeof(int)      , 1    , f_out) ;
    fwrite(&this->n_basis_to    , sizeof(int)      , 1    , f_out) ;
  } else if (binary==2) {
    fprintf(f_out,"%%%%MatrixMarket matrix coordinate real general \n");
    fprintf(f_out," %d " , this->n_basis_to  );
    fprintf(f_out," %d " , this->n_basis_from);
    fprintf(f_out," %d\n " , total_Transitions);
  }

  else {
    fprintf(f_out," %d\n " , total_Transitions);
    fprintf(f_out," %d\n " , this->n_basis_from);
    fprintf(f_out," %d\n " , this->n_basis_to  );
  }
   
  
  for(int k=0; k<total_Transitions; k++) {
    if(k==0) {
      newTrans  = (Hxx_transition * )  Transitions_Orderer->getLowest_and_remove(depth);
    } else {
      newTrans = (Hxx_transition * )  Transitions_Orderer->Again_getLowest_and_remove(depth);
    }
    if(binary==1) {
      int i_to(newTrans->get_to());
      int i_from(newTrans->get_from());
      Hxx_FLOAT f_coeff(newTrans->get_coeff());
      fwrite(&i_to   , sizeof(int)      , 1    , f_out) ;
      fwrite(&i_from , sizeof(int)      , 1    , f_out) ;
      fwrite(&f_coeff, sizeof(Hxx_FLOAT), 1    , f_out) ;
    } else if (binary==2) {
      fprintf(f_out," %d %d  %20.20e\n", newTrans->get_to()+1,newTrans->get_from()+1,newTrans->get_coeff()  );
    } else {
      fprintf(f_out," %d %d  %20.20e\n", newTrans->get_to(),newTrans->get_from(),newTrans->get_coeff()  );
    }
  }
  fclose(f_out);
  this->free_space();
}  


void Inizializza_Hilbertxx(){

  std::cout << first_bit_mask_ << std::endl;
  std::cout << first_bit_mask << std::endl;
  
  if(!first_bit_mask_ ) {
    std::cout << " ma neanca gli static static mi inizializza !!!\n";
    exit(0);
  }

  first_bit_mask  =  2* first_bit_mask_;

  std::cout << " inizializzo i fermions \n";
  Hxx_1p_FermionicOperators::operators = Hxx_1p_FermionicOperators ::Initialise();
  std::cout << " inizializzo scrambler \n";
  Hxx_Utilities::SCRAMBLE_FOR_BITS = Hxx_Utilities::create_scrambler(sizeof(Hxx_TypeForBits)*8);
  std::cout << " inizializzo scrambler OK \n";
}

 
Hxx_1p_FermionicOperator_ * Hxx_1p_FermionicOperators::operators = Hxx_1p_FermionicOperators ::Initialise();

int * Hxx_Utilities::SCRAMBLE_FOR_BITS = Hxx_Utilities::create_scrambler(sizeof(Hxx_TypeForBits)*8);









//  #else ifdef FORMAT_PETER 
  
//    fprintf(f_out," %d\n " , n_Abasis);
//    fprintf(f_out," %d\n " , total_Transitions);
  
//    int from_limits[n_Abasis];
//    int from=0;
//    int new_from=0;
//    int count=0;
//    for(int k=0; k<total_Transitions; k++) {
//      if(k==0) {
//        newTrans  = (Hxx_transition * )  Transitions_Orderer->getLowest_and_remove(depth);
//      } else {
//        newTrans = (Hxx_transition * )  Transitions_Orderer->Again_getLowest_and_remove(depth);
//      }
//      new_from = newTrans->get_from();
//      if( new_from != from ) {
//        for(int j=from; j<new_from; j++) {
//  	from_limits[j]=count;
//        }
//        from = new_from;
//      }
//      #ifdef FORMAT_BINARY

//      int i_to(newTrans->get_to());
//      Hxx_FLOAT f_coeff(newTrans->get_coeff());

//      fwrite(&i_to, sizeof(int), 1    , f_out) ;
//      fwrite(&f_coeff, sizeof(Hxx_FLOAT), 1    , f_out) ;

//      #else

//      fprintf(f_out,"%d  %e\n", newTrans->get_to(),newTrans->get_coeff()  );

//      #endif



//  for(int i=0; i<n_Abasis-1; i++) {
//    fprintf(f_out,"%d\n", from_limits[i]);
//    }
//    fprintf(f_out,"%d\n",count );

//  #endif




