calc_raddist.cpp

#include <iostream>
#include <cassert>
#include <string>
#include <vector>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include "mvector3.hpp"

using namespace std;

const double cutof_leng = 1.4;
const double radial_div = 0.1;
const int    n_rad_div  = (int)(cutof_leng / radial_div);

enum{
  CLUSTER_N = 2,
  MOL_NUMB  = 4,
};

double3 grid_leng(cutof_leng);
int grid_numb[3] = { -1 };
int all_grid_numb( -1 );
double3 L( -1.0 );

struct prtcl{
  double3 r;
  int prop, chem;
  int l_idx, unit;
  int prtcl_in_cell;
  bool outer;
}; 

struct ByPropPred{
  bool operator()(const prtcl& p1, const prtcl& p2) const {
    return p1.prop < p2.prop;
  }
};

struct ByLipidPred{
  bool operator()(const prtcl& p1, const prtcl& p2) const {
    return (MOL_NUMB * p1.l_idx + p1.unit) < (MOL_NUMB * p2.l_idx + p2.unit);
  }
};

void LoadParticleDat(ifstream& fin, vector<prtcl>& Prtcl, const int sys_size){
  prtcl buf_p;
  double3 buf_d3;
  buf_p.outer = false;
  
  for(int i=0; i<sys_size; i++){
    fin >> buf_p.r.x >> buf_p.r.y >> buf_p.r.z >> buf_d3.x >> buf_d3.y >> buf_d3.z >> buf_p.prop >> buf_p.chem >> buf_p.unit >> buf_p.l_idx;
    if(buf_p.prop == 0){
      Prtcl.push_back(buf_p);
    }else if(buf_p.chem && (buf_p.unit == 0)){
      Prtcl.push_back(buf_p);
    }else if(buf_p.chem && (buf_p.unit == 1)){
      Prtcl.push_back(buf_p);
    }
  }
}

void SortParticleDat(vector<prtcl>& Prtcl, const int wN){
  sort(Prtcl.begin(), Prtcl.end(), ByPropPred());
  sort(Prtcl.begin() + wN, Prtcl.end(), ByLipidPred());
}

void LoadScale(ifstream& fin){
  double buf[6] = {0.0};
  fin >> L.x >> L.y >> L.z >> buf[0] >> buf[1] >> buf[2] >> buf[3] >> buf[4] >> buf[5];
}

void SetScale(){
  grid_numb[0] = static_cast<int>(L.x / cutof_leng);
  grid_numb[1] = static_cast<int>(L.y / cutof_leng);
  grid_numb[2] = static_cast<int>(L.z / cutof_leng);
  all_grid_numb = grid_numb[0] * grid_numb[1] * grid_numb[2];
  grid_leng.x = L.x / grid_numb[0];
  grid_leng.y = L.y / grid_numb[1];
  grid_leng.z = L.z / grid_numb[2];
}

int GenHash(const double3& r){
  int i = (int)(r.x / grid_leng.x);
  int j = (int)(r.y / grid_leng.y);
  int k = (int)(r.z / grid_leng.z);
 
  if(i == grid_numb[0]) i--;
  if(j == grid_numb[1]) j--;
  if(k == grid_numb[2]) k--;
  
  const int hash = i + grid_numb[0] * (j + k * grid_numb[1]);
  assert(hash < all_grid_numb);
  assert(hash >= 0);
  return hash;
}

int GenHash(int* q){
  if(q[0] < 0 || q[0] > grid_numb[0]) q[0] -= (int)q[0]*grid_numb[0]/fabs(q[0]);
  if(q[1] < 0 || q[1] > grid_numb[1]) q[1] -= (int)q[1]*grid_numb[1]/fabs(q[1]);
  if(q[2] < 0 || q[2] > grid_numb[2]) q[2] -= (int)q[2]*grid_numb[2]/fabs(q[2]);
  const int hash = q[0] + grid_numb[0]*(q[1] + q[2]*grid_numb[1]);
  return hash;
}

void MinImage(double3& dr){
  dr.x -= L.x * round(dr.x / L.x);
  dr.y -= L.y * round(dr.y / L.y);
  dr.z -= L.z * round(dr.z / L.z);  
}

void GetNearImage(double& dr, double& rj, double L){
  if(dr < -0.5 * L){
    rj += L; dr += L;
  }else if(dr > 0.5 * L){
    rj -= L; dr -= L;
  }
}

void MinImage(double3& dr, double3& rj){
  GetNearImage(dr.x, rj.x, L.x);
  GetNearImage(dr.y, rj.y, L.y);
  GetNearImage(dr.z, rj.z, L.z);
}

void CreateNearlist(vector<vector<int> > &ilist){
  ilist.resize(all_grid_numb);
  for(int i=0; i<all_grid_numb; i++){
    ilist[i].resize(27);
  }
  
  int i_num=0;
  for(int iz=0; iz<grid_numb[2]; iz++){
    for(int iy=0; iy<grid_numb[1]; iy++){
      for(int ix=0; ix<grid_numb[0]; ix++){
	int j_num=0;
	for(int jz=-1; jz<2; jz++){
	  for(int jy=-1; jy<2; jy++){
	    for(int jx=-1; jx<2; jx++){
	      int q[3] = {ix+jx, iy+jy, iz+jz};
	      ilist[i_num][j_num]=GenHash(q);
	      j_num++;
	    }
	  }
	}
	i_num++;
      }
    }
  }
}

void RegistPrtclIdx(vector<int> &prtcl_idx, vector<prtcl> &Prtcl, vector<int> &elem_in_grid, vector<int> &adrs_grid){
  elem_in_grid.assign(elem_in_grid.size(), 0);
  adrs_grid[0] = 0;
  for(size_t i=0; i<Prtcl.size(); i++){
    const int hash = GenHash(Prtcl[i].r);
    Prtcl[i].prtcl_in_cell = hash;
    elem_in_grid[hash]++;
  }
  for(int i=0; i<all_grid_numb; i++){
    adrs_grid[i+1] = adrs_grid[i] + elem_in_grid[i];
  }
  for(size_t i=0; i<Prtcl.size(); i++){
    const int hash = Prtcl[i].prtcl_in_cell;
    const int temp_idx = adrs_grid[hash];
    prtcl_idx[temp_idx] = i;
    adrs_grid[hash]++;
  }
  for(int i=0; i<all_grid_numb; i++){
    adrs_grid[i] = adrs_grid[i] - elem_in_grid[i];
  }
}

void CheckReg(const vector<int>& adrs_grid, const vector<int> &prtcl_idx, const vector<prtcl> &Prtcl){
  for(int tar=0; tar<all_grid_numb; tar++){
    for(int pi = adrs_grid[tar]; pi<adrs_grid[tar+1]; pi++){
      const int hash = GenHash(Prtcl[prtcl_idx[pi]].r);
      assert(hash == tar);
    }
  }
}

void GenerateGraph(vector<prtcl>& Prtcl, vector<int>& prtcl_idx, vector<int>& adrs_grid, vector<vector<int> >& n_prtcl_id, vector<vector<int> >& ilist)
{
  const size_t size = Prtcl.size();
  for(size_t idx=0; idx<size; idx++){
    const int pi     = prtcl_idx[idx];
    const int prop_i = Prtcl[pi].prop;
    const int hash   = Prtcl[pi].prtcl_in_cell;
    const double3 ri = Prtcl[pi].r;
    for(int i=0; i<27; i++){
      const int intr_grid = ilist[hash][i];
      const int beg_id    = adrs_grid[intr_grid    ];
      const int end_id    = adrs_grid[intr_grid + 1];
      for(int id=beg_id; id < end_id; id++){
	const int pj     = prtcl_idx[id];
	const int prop_j = Prtcl[pj].prop;
	double3 rj = Prtcl[pj].r;
	double3 dr = rj - ri;
	MinImage(dr);
	const double dist = dr.norm2();
	if((dist < cutof_leng) && (pi != pj) && (prop_i == 1) && (prop_j == 1)){
	  n_prtcl_id[pi].push_back(pj);
	}
      }
    }
  }
}


void Dfs(int pid, vector<bool>& reg_flag, vector<vector<int> >& n_prtcl_id, vector<int>& patch_id, vector<prtcl>& Prtcl, int cur_patch){
  if(reg_flag[pid] || (Prtcl[pid].prop != 1)) return;
  reg_flag[pid] = true;
  patch_id[pid] = cur_patch;
  const int n_prtcl_size = n_prtcl_id[pid].size();
  for(int i=0; i<n_prtcl_size; i++){
    const int	pjd = n_prtcl_id[pid][i];
    double3	dr  = Prtcl[pjd].r - Prtcl[pid].r;
    MinImage(dr, Prtcl[pjd].r);
    Dfs(pjd, reg_flag, n_prtcl_id, patch_id, Prtcl, cur_patch);
  }
}

int CountRegN(const vector<bool>& reg_flag){
  int n=0;
  for(size_t i=0; i<reg_flag.size(); i++) n += reg_flag[i];
  return n;
}

int SearchNextPid(vector<bool>& reg_flag, const int size, vector<prtcl>& Prtcl){
  bool flag = true;
  int ret = 0;
  while(flag){
    int pid = rand() % size;
    if((!reg_flag[pid]) && (Prtcl[pid].prop == 1)){
      ret = pid;
      flag = false;
    }
  }
  return ret;
}

void SearchPatch(vector<bool>& reg_flag, vector<vector<int> >& n_prtcl_id, vector<int>& patch_id, vector<prtcl>& Prtcl, int& cur_patch, const int lN)
{
  int    regist_n = 0;
  int    pid      = 0;
  const int size  = reg_flag.size();
  while(regist_n != lN){
    pid = SearchNextPid(reg_flag, size, Prtcl);
    Dfs(pid, reg_flag, n_prtcl_id, patch_id, Prtcl, cur_patch);
    regist_n = CountRegN(reg_flag);
    cout << regist_n << endl;
    cur_patch++;
  }
}

double3 CalcCMpos(const vector<prtcl>& Prtcl){
  double3 cmpos(0.0); int elem = 0;
  for(size_t i=0; i<Prtcl.size(); i++){
    if(Prtcl[i].prop == 1){
      cmpos += Prtcl[i].r;
      elem++;
    }
  }
  cmpos /= elem;
  return cmpos;
}

void DumpInOutElem(const vector<int>& patch_id, ofstream& fout, const int cur_patch, const int cur_time){
  fout << cur_time << " ";
  int ret[CLUSTER_N] = {0};
  for(int i=0; i<cur_patch; i++) ret[i] = count(patch_id.begin(), patch_id.end(), i);
  sort(ret, ret + CLUSTER_N);
  for(int i=0; i<CLUSTER_N; i++) fout << ret[i] << " ";
  fout << endl;
}

void CalcRadialDist(const vector<prtcl>& Prtcl, vector<int>& prtcl_idx, vector<int>& adrs_grid, vector<vector<int> > ilist,
		    const int wN, vector<int>& wn_in_shell, vector<int>& hn_in_shell)
{
  const int size = Prtcl.size();
  for(int i=wN; i<size; i += 2){
    if(Prtcl[i].outer){
      const int hash = Prtcl[i+1].prtcl_in_cell;
      const double3 ri = Prtcl[i+1].r;
      for(int cell=0; cell<27; cell++){
	const int intr_grid = ilist[hash][cell];
	const int beg_id = adrs_grid[intr_grid  ];
	const int end_id = adrs_grid[intr_grid+1];
	for(int id=beg_id; id<end_id; id++){
	  const int pj     = prtcl_idx[id];
	  const int prop_j = Prtcl[pj].prop;
	  double3 rj = Prtcl[pj].r;
	  double3 dr = rj - ri;
	  MinImage(dr);
	  const double dist = dr.norm2();
	  if(dist < cutof_leng){
	    if(prop_j == 0){
	      const int shell_id = (int)(dist / radial_div);
	      wn_in_shell[shell_id]++;
	    }else if(prop_j == 1 && Prtcl[pj].outer){
	      const int shell_id = (int)(dist / radial_div);
	      hn_in_shell[shell_id]++;
	    }
	  }
	}
      }
    }
  }
}

void LabelOutOrIn(const double3 cmpos, vector<prtcl>& Prtcl, const vector<int>& patch_id, const int wN)
{
  const int size = Prtcl.size();
  int in_out_elem[2] = {0};
  for(int i=0; i<size; i++){
    const int id = patch_id[i];
    in_out_elem[id]++;
  }
  
  int tar_patch_id = -1;
  if(in_out_elem[0] < in_out_elem[1]){
    tar_patch_id = 1;
  }else{
    tar_patch_id = 0;    
  }
  
  for(int i=wN; i<size; i++){
    if(patch_id[i] == tar_patch_id){
      Prtcl[i].outer = true;
    }
  }
}

void DumpRadialDiv(ofstream& fout, vector<int>& wn_in_shell, vector<int>& hn_in_shell, const int sample){
  for(size_t i=0; i<wn_in_shell.size(); i++){
    const double Rin  = i     * radial_div;
    const double Rout = (i+1) * radial_div;
    const double shell_vol = 4.0 * M_PI / 3.0 * (Rout*Rout*Rout - Rin*Rin*Rin);
    fout << (i + 0.5) * radial_div << " " <<  (double)wn_in_shell[i] / (double)sample / shell_vol << " " << (double)hn_in_shell[i] / (double)sample / shell_vol << endl;
  }
}

void AssertFileOpen(const istream& ist, const string& str){
  if(!ist){
    cerr << "Cannot open " << str << endl;
    assert(ist);
  }
}

int main(int argc,char* argv[]){
  assert(argc == 3);
  srand((unsigned)time(NULL));
  const int beg_time = atoi(argv[2]);
  string s_cur  = argv[1];
  const string fname[6] = {
    s_cur + "/particle_data.txt",
    s_cur + "/macro_data.txt",
    s_cur + "/inout_elem.txt",
    s_cur + "/rad_water_dense.txt",
    s_cur + "/mean_rad.txt",
    s_cur + "/radial_dist.txt"
  };
  
  ifstream f_p(fname[0].c_str());
  ifstream f_m(fname[1].c_str());
  ofstream fout0(fname[2].c_str());
  ofstream fout1(fname[3].c_str());
  ofstream fout2(fname[4].c_str());
  ofstream fout3(fname[5].c_str());

  AssertFileOpen(f_p, fname[0]);
  AssertFileOpen(f_m, fname[1]);
  
  int wN, hbN, all_time, time_step;
  double buf;
  f_m >> wN >> hbN >> L.x >> buf >> all_time >> time_step;
  const int sys_size = wN + hbN;
  const int lN       = hbN / MOL_NUMB;
  L.y = L.x;
  L.z = L.x;

  SetScale();

  //nearest cell
  vector<vector<int> > ilist;
  
  vector<int> elem_in_grid(2 * all_grid_numb, 0);
  vector<int> adrs_grid(2 * all_grid_numb + 1, 0);
  
  vector<prtcl> Prtcl;
  vector<vector<int> > n_prtcl_id;
  vector<int> patch_id;
  vector<int> prtcl_idx;
  vector<bool> reg_flag;
  vector<int> wn_in_shell(n_rad_div, 0), hn_in_shell(n_rad_div, 0);

  CreateNearlist(ilist);
  
  for(int time=0; time < all_time; time += time_step){
    cout << time << endl;
    LoadParticleDat(f_p, Prtcl, sys_size);
    SortParticleDat(Prtcl, wN);
    
    if(time >= beg_time){
      const size_t p_size = Prtcl.size();
      prtcl_idx.resize(p_size,-1);
      patch_id.resize(p_size,-1);
      n_prtcl_id.resize(p_size);
      
      RegistPrtclIdx(prtcl_idx, Prtcl, elem_in_grid, adrs_grid);
      CheckReg(adrs_grid, prtcl_idx, Prtcl);
      
      GenerateGraph(Prtcl, prtcl_idx, adrs_grid, n_prtcl_id, ilist);
      
      reg_flag.resize(p_size, false);
      int   cur_patch = 0;
      SearchPatch(reg_flag, n_prtcl_id, patch_id, Prtcl, cur_patch, lN);
      
      if(cur_patch == CLUSTER_N){
	DumpInOutElem(patch_id, fout0, cur_patch, time);
	const double3 cm_pos = CalcCMpos(Prtcl);
	LabelOutOrIn(cm_pos, Prtcl, patch_id, wN);
	CalcRadialDist(Prtcl, prtcl_idx, adrs_grid, ilist, wN, wn_in_shell, hn_in_shell);	
      }
    }
    Prtcl.clear();
    prtcl_idx.clear();
    reg_flag.clear();
    for(size_t i=0; i<n_prtcl_id.size(); i++) n_prtcl_id[i].clear();
  }
  const int samplen = (all_time - beg_time) / time_step;
  DumpRadialDiv(fout3, wn_in_shell, hn_in_shell, samplen);
  
  return 0;
}