//  WOLFF.C
//  Wolff cluster algorithm
//  © 2021 Peter J. Meyer

#include "iss.h"

#define STACK_POINTER_CHECK true
#define DEBUG false

unsigned int num_spins_visited_this_run;    //  set to 0 in EQUIL.C at start of each run

char wolff_map_filepath[128];               //  pathname of map file

short int w_i, w_j, w_k, w_l, w_r, r0;
//  not static, since accessed in DISP_LAT.C

static int w_spin_value, w_new_spin_value;

static int perform_wolff_update(void);
static void grow_and_flip_wolff_cluster(void);
static void flip_wolff_spin(void);
static int possible_new_wolff_spin(void);

//  This implementation is slow for T considerably larger than
//  the critical temperature (e.g. T = 4.0), but is about the
//  same speed as the Swendsen-Wang algorithm at the critical temperature.

/*--------------------------*/
void perform_wolff_sweep(void)
{
unsigned int target_num, target_num_spins_visited_this_run;

if ( model_is_ising )
    target_num_spins_visited_this_run = ( mcs_num+1 )*((unsigned int)num_spins/2);
else    //  model is q-state Potts
    target_num_spins_visited_this_run = ( mcs_num+1 )*(((unsigned int)num_spins*(q-1))/q);

if (  target_num_spins_visited_this_run < num_spins_visited_this_run )
    target_num = 0;
else
    target_num = target_num_spins_visited_this_run - num_spins_visited_this_run;

do  {
    num_spins_visited_this_sweep += perform_wolff_update();
    num_clusters_traced_this_sweep++;
    if ( map_file_requested )
        {
        strcpy(wolff_map_filepath,map_filepath);
        display_map(wolff_map_filepath,4);
        }
    } while ( num_spins_visited_this_sweep < target_num );

num_spins_visited_this_run += num_spins_visited_this_sweep;
}

/*---------------------------------*/
static int perform_wolff_update(void)
{
mark_all_sites_unvisited();

switch ( dimensionality )
    {
    case 2:
    //  Select spin at random.
    do  {
        w_i = (int)(PRNG1()*size);
        w_j = (int)(PRNG1()*size);
        } while ( !SITE_IS_OCCUPIED_2(w_i,w_j) );

    w_spin_value = spin_at_occupied_site(w_i,w_j);
    break;

    case 3:
    do  {
        w_i = (int)(PRNG1()*size);
        w_j = (int)(PRNG1()*size);
        w_k = (int)(PRNG1()*size);
        } while ( !SITE_IS_OCCUPIED_3(w_i,w_j,w_k) );

    w_spin_value = spin_at_occupied_site(w_i,w_j,w_k);
    break;

    case 4:
    do  {
        w_i = (int)(PRNG1()*size);
        w_j = (int)(PRNG1()*size);
        w_k = (int)(PRNG1()*size);
        w_l = (int)(PRNG1()*size);
        } while ( !SITE_IS_OCCUPIED_4(w_i,w_j,w_k,w_l) );

    w_spin_value = spin_at_occupied_site(w_i,w_j,w_k,w_l);
    break;
    }

if ( model_is_q_state_potts )
    {
    //  Pick a new q-state Potts spin value
    //  randomly selected from values other than the old.
    do  {
        w_new_spin_value = (int)(PRNG1()*q) + 1;
        } while ( w_new_spin_value == w_spin_value );
    }

clear_stack();

ct_size_of_cluster = 0;

grow_and_flip_wolff_cluster();   

if ( ct_size_of_cluster > max_cluster_size )
    max_cluster_size = ct_size_of_cluster;

#if STACK_POINTER_CHECK
check_stack_is_empty(__FILE__);
#endif

return ( ct_size_of_cluster );
}

/*
Grow and flip a Wolff cluster (call with 1st spin at i,j):

    Flip spin.
    Put r0 = -1

    Loop
        While there is a direction r > r0 such that there is a nn spin in direction r
        which has the same value as the seed spin and is unvisited:
            If PRNG < probability
                Push i,j,r on stack
                Put i,j = location of new spin
                Put r0 = -1
                Flip spin
            End If
        End While

        If stack is empty
            finish
        End If

        Pop i,j,r0 from stack
    End Loop

Flip spin:

    Flip the spin
    Mark spin as visited
    Increase cluster count        
 */

/*-----------------------------------------*/
static void grow_and_flip_wolff_cluster(void)
{
flip_wolff_spin();
r0 = -1;

switch ( dimensionality )
    {
    case 2:
    if ( map_file_requested )
        //  Erase v_bonds2
        memset(&v_bonds2[0][0],0,num_sites);

    loop
        {
        while ( possible_new_wolff_spin() )
            {
            //  direction of possible new spin is w_r
            //  and location of possible new spin is (nn_i,nn_j).
            if ( PRNG1() < v_bond_probability )
                {
                if ( push_onto_stack_with_r(w_r,w_i,w_j) )
                    {
                    w_i = nn_i;
                    w_j = nn_j;
                    r0 = -1;
                    flip_wolff_spin();
                    #if DEBUG
                    if ( map_file_requested )
                        {
                        strcpy(wolff_map_filepath,map_filepath);
                        display_map(wolff_map_filepath,4);
                        }
                    #endif
                    }
                }
            else
                r0 = w_r;
            }

        if ( stack_empty() )
            return;

        pop_from_stack_with_r(&r0,&w_i,&w_j);
        }
    break;

    case 3:
    if ( map_file_requested )
        //  Erase v_bonds3
        memset(&v_bonds3[0][0][0],0,num_sites);

    loop
        {
        while ( possible_new_wolff_spin() )
            {
            //  direction of possible new spin is w_r
            //  and location of possible new spin is (nn_i,nn_j,nn_k).
            if ( PRNG1() < v_bond_probability )
                {
                if ( push_onto_stack_with_r(w_r,w_i,w_j,w_k) )
                    {
                    w_i = nn_i;
                    w_j = nn_j;
                    w_k = nn_k;
                    r0 = -1;
                    flip_wolff_spin();
                    #if DEBUG
                    if ( map_file_requested )
                        {
                        strcpy(wolff_map_filepath,map_filepath);
                        display_map(wolff_map_filepath,4);
                        }
                    #endif
                    }
                }
            else
                r0 = w_r;
            }

        if ( stack_empty() )
            return;

        pop_from_stack_with_r(&r0,&w_i,&w_j,&w_k);
        }
    break;

    case 4:
    if ( map_file_requested )
        //  Erase v_bonds3
        memset(&v_bonds4[0][0][0][0],0,num_sites);

    loop
        {
        while ( possible_new_wolff_spin() )
            {
            //  direction of possible new spin is w_r
            //  and location of possible new spin is (nn_i,nn_j,nn_k,nn_l).
            if ( PRNG1() < v_bond_probability )
                {
                if ( push_onto_stack_with_r(w_r,w_i,w_j,w_k,w_l) )
                    {
                    w_i = nn_i;
                    w_j = nn_j;
                    w_k = nn_k;
                    w_l = nn_l;
                    r0 = -1;
                    flip_wolff_spin();
                    #if DEBUG
                    if ( map_file_requested )
                        {
                        strcpy(wolff_map_filepath,map_filepath);
                        display_map(wolff_map_filepath,4);
                        }
                    #endif
                    }
                }
            else
                r0 = w_r;
            }

        if ( stack_empty() )
            return;

        pop_from_stack_with_r(&r0,&w_i,&w_j,&w_k,&w_l);
        }
    break;
    }
}

/*-----------------------------*/
static void flip_wolff_spin(void)
{
ct_size_of_cluster++;

switch ( dimensionality )
    {
    case 2:
    //  Flip the spin.
    if ( model_is_ising )
        do_ising_spin_flip(w_spin_value,w_i,w_j); 
    else //  model is q-state Potts
        do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j); 

    MARK_SITE_AS_VISITED_2(w_i,w_j);
    break;

    case 3:
    if ( model_is_ising )
        do_ising_spin_flip(w_spin_value,w_i,w_j,w_k); 
    else //  model is q-state Potts
        do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j,w_k); 

    MARK_SITE_AS_VISITED_3(w_i,w_j,w_k);
    break;

    case 4:
    if ( model_is_ising )
        do_ising_spin_flip(w_spin_value,w_i,w_j,w_k,w_l); 
    else //  model is q-state Potts
        do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j,w_k,w_l); 

    MARK_SITE_AS_VISITED_4(w_i,w_j,w_k,w_l);
    break;
    }
}

//  This returns true if there is a direction r > r0 such that
//  there is a nn spin in direction r which has the same spin value
//  as the seed spin and which is unvisited.
//  In this case the direction of the possible new spin on return is in w_r
//  and the location of the possible new spin is (nn_i,nn_j,...).
/*------------------------------------*/
static int possible_new_wolff_spin(void)
{
if ( r0 == coord_num - 1 )
    return ( false );

switch ( dimensionality )
    {
    case 2:
    if ( !bonds2[w_i][w_j] )
        return ( false );

    //  There is at least one open bond at this site.

    get_nn_sites(w_i,w_j);

    for ( w_r=r0+1; w_r<coord_num; w_r++ )
        {
        if ( bonds2[w_i][w_j] & bit[w_r] )  
        //  if an open bond exists in direction w_r
            {
            nn_i = nn_sites[w_r][0];
            nn_j = nn_sites[w_r][1];
            if ( !SITE_IS_VISITED_2(nn_i,nn_j) )
            //  if the nn spin has not been visited (=flipped)
                {
                //  This spin has not been visited, so it has not been flipped
                //  so it has the same value as when we began to grow the cluster.
                if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j) )
                //  if this spin has the same spin value as the seed spin
                    {
                    if ( map_file_requested )
                        {
                        //  Open the virtual bond to show that 
                        //  this bond has been considered.
                        v_bonds2[w_i][w_j] |= bit[w_r];
                        //  Mark the virtual bond as open at the other site.
                        v_bonds2[nn_i][nn_j] 
                            |= bit[inverse_direction((w_i+w_j)%2,w_r)];                     
                        }
                    return ( true );
                    }
                }
            }
        }
    break;

    case 3:
    if ( !bonds3[w_i][w_j][w_k] )
        return ( false );

    //  There is at least one open bond at this site.

    get_nn_sites(w_i,w_j,w_k);

    for ( w_r=r0+1; w_r<coord_num; w_r++ )
        {
        if ( bonds3[w_i][w_j][w_k] & bit[w_r] )  
        //  if an open bond exists in direction w_r
            {
            nn_i = nn_sites[w_r][0];
            nn_j = nn_sites[w_r][1];
            nn_k = nn_sites[w_r][2];
            if ( !SITE_IS_VISITED_3(nn_i,nn_j,nn_k) )
            //  if the nn spin has not been visited (=flipped)
                {
                //  This spin has not been visited, so it has not been flipped
                //  so it has the same value as when we began to grow the cluster.
                if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j,nn_k) )
                //  if this spin has the same spin value as the seed spin
                    {
                    if ( map_file_requested )
                        {
                        //  Open the virtual bond to show that 
                        //  this bond has been considered.
                        v_bonds3[w_i][w_j][w_k] |= bit[w_r];
                        //  Mark the virtual bond as open at the other site.
                        v_bonds3[nn_i][nn_j][nn_k] 
                            |= bit[inverse_direction((w_i+w_j+w_k)%2,w_r)];                     
                        }
                    return ( true );
                    }
                }
            }
        }
    break;

    case 4:
    if ( !bonds4[w_i][w_j][w_k][w_l] )
        return ( false );

    //  There is at least one open bond at this site.

    get_nn_sites(w_i,w_j,w_k,w_l);

    for ( w_r=r0+1; w_r<coord_num; w_r++ )
        {
        if ( bonds4[w_i][w_j][w_k][w_l] & bit[w_r] )  
        //  if an open bond exists in direction w_r
            {
            nn_i = nn_sites[w_r][0];
            nn_j = nn_sites[w_r][1];
            nn_k = nn_sites[w_r][2];
            nn_l = nn_sites[w_r][3];
            if ( !SITE_IS_VISITED_4(nn_i,nn_j,nn_k,nn_l) )
            //  if the nn spin has not been visited (=flipped)
                {
                //  This spin has not been visited, so it has not been flipped
                //  so it has the same value as when we began to grow the cluster.
                if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j,nn_k,nn_l) )
                //  if this spin has the same spin value as the seed spin
                    {
                    if ( map_file_requested )
                        {
                        //  Open the virtual bond to show that 
                        //  this bond has been considered.
                        v_bonds4[w_i][w_j][w_k][w_l] |= bit[w_r];
                        //  Mark the virtual bond as open at the other site.
                        v_bonds4[nn_i][nn_j][nn_k][nn_l]
                            |= bit[inverse_direction((w_i+w_j+w_k+w_l)%2,w_r)];                     
                        }
                    return ( true );
                    }
                }
            }
        }
    break;
    }

return ( false );
}