ga_grow.c

来自「药物开发中的基于结构的从头设计代码」· C语言 代码 · 共 879 行 · 第 1/2 页

	// copy these elites to the next population

	for(i=0;i<count;i++) member[i]=pool[i];

	for(i=count;i<num_member;i++) member[i].valid=0;	

	num_member=count;

	printf("%d elite molecules are selected\n",count);

	delete [] index;
	delete [] pool;

	return;
}

void Population::Select_Parents(int num_wanted)	
{
	extern Parameter *parm;
	int i,j,tmp,mark,id;
	int total,num_candidate,num_left;
	float sim;
	struct Wheel
		{
		 int valid;
		 int id;
		 int min;
		 int max;
		} *wheel;

	if(num_member==0) {num_parent=0; return;}

	wheel=new Wheel[num_member];
	if(wheel==NULL) Memory_Allocation_Error();

	// initialize the wheel, according to the wins

	 total=0; num_candidate=0; 

	 for(i=0;i<num_member;i++)
		{
		 if(member[i].valid==0) continue;
		 else
			{
			 wheel[num_candidate].valid=1;
			 wheel[num_candidate].id=i+1;
		 	 wheel[num_candidate].min=total;
		 	 total+=member[i].wins;
		 	 wheel[num_candidate].max=total;
			 num_candidate++;
			}
		}

	if(num_candidate<=num_wanted) 	// not enough molecules for selection
		{
		 num_parent=0;

		 for(i=0;i<num_member;i++)
			{
			 if(member[i].valid==0) continue;
			 else
				{
				 parent[num_parent]=member[i];
				 num_parent++;
				}
			}

		 delete [] wheel;

		 printf("%d parent molecules are selected\n", num_parent);
		
		 return;
		}

	// select the parents into the mating pool

	num_parent=0; num_left=num_candidate;

	do
	{
	 tmp=(int)(drand48()*total);

	 for(i=0;i<num_candidate;i++)
		{
		 if(wheel[i].valid==0) continue;
		 else if(tmp<wheel[i].min) continue;
		 else if(tmp>=wheel[i].max) continue;
		 else 
			{
			 id=wheel[i].id;

			 // first, similarity check, i.e. niching algorithm

			 mark=FALSE;

			 for(j=0;j<num_parent;j++)
			{
			 sim=Cal_Molecular_Similarity(member[id-1],parent[j]);
			 if(sim<parm->similarity_cutoff) continue;
			 else {mark=TRUE; break;}
			}

			 if(mark==FALSE)	// record the selected molecule
				{
				 parent[num_parent]=member[id-1];
				 num_parent++;
				}
		
			 // delete the current piece and re-arrange the wheel
		
			 wheel[i].valid=0; num_left--;

			 total-=member[id-1].wins;

			 if(i==(num_candidate-1)) break;

			 for(j=i+1;j<num_candidate;j++)
				{
				 wheel[j].min-=member[id-1].wins;
				 wheel[j].max-=member[id-1].wins;
				}
			 break;
			}
		}

	} while((num_parent<num_wanted)&&(num_left>0));

	delete [] wheel;

	printf("%d parent molecules are selected\n", num_parent);
 
	return;
}

void Population::Initialize(Ligand core, int num_generation)
{
	extern Parameter *parm;
	extern FragLibrary *fraglib;
	int i,j,k,tmp,mark,count;
	int max_record,num_core_seed,num_frag_seed,num_frag,num_new;
	Ligand lig,frag,result[24];
	int *core_seed_list,*frag_seed_list,*frag_list,*tmp_seed_list;
	Growing_Record tmp_record;

	core_seed_list=new int[core.num_atom];
	if(core_seed_list==NULL) Memory_Allocation_Error();

	frag_seed_list=new int[500];
	if(frag_seed_list==NULL) Memory_Allocation_Error();

	tmp_seed_list=new int[500];
	if(tmp_seed_list==NULL) Memory_Allocation_Error();

	frag_list=new int[fraglib->num_frag];
	if(frag_list==NULL) Memory_Allocation_Error();

	// first, get the seed hydrogens on the core structure

	num_core_seed=core.Make_Seed_List(core_seed_list);

	// growing on each seed hydrogen on the core structure

	max_record=parm->max_population*fraglib->num_frag;

        growing_record=new Growing_Record[max_record];
        if(growing_record==NULL) Memory_Allocation_Error();

	num_member=num_record=0;

	lig=core;

	while(1)
	{
	 num_new=0;

	 // first, estimate the number of fragments which will be
	 // grown on each growing site and make the fragment list.

	 num_frag=(parm->max_population-num_member)/num_core_seed;

	 for(i=0;i<num_core_seed;i++)
		{
	 	 num_frag=fraglib->Make_Fragment_List(num_frag,frag_list);

		 // applying each fragment in the list

		 for(j=0;j<num_frag;j++)
			{
			 mark=fraglib->Get_A_Fragment(frag_list[j],frag);
			 if(mark==FALSE) continue; // get fragment failed 

			 num_frag_seed=frag.Make_Seed_List(frag_seed_list);

			 for(k=0;k<num_frag_seed;k++)
			{
			 tmp_record.ligand_id=1;
			 tmp_record.lig_seed=core_seed_list[i];
			 tmp_record.fragment_id=frag_list[j];
			 tmp_record.frag_seed=frag_seed_list[k];

			 mark=Growing_Record_Check(tmp_record);
			 if(mark==FALSE) continue;
			 else frag_seed_list[k]=0;  // has been tried before 
			}

			 tmp=0;		 
			 for(k=0;k<num_frag_seed;k++)
			{
			 if(frag_seed_list[k]==0) continue;
			 else {tmp_seed_list[tmp]=frag_seed_list[k]; tmp++;}
			}

			 if(tmp==0) 
				{
				 // puts("no available seed H on the fragment");
				 continue;
				}
			 else k=(int)(drand48()*tmp);

                         tmp_record.ligand_id=1;
                         tmp_record.lig_seed=core_seed_list[i];
                         tmp_record.fragment_id=frag_list[j];
                         tmp_record.frag_seed=tmp_seed_list[k];

			 growing_record[num_record]=tmp_record;
			 if(num_record<(max_record-1)) num_record++;

			 mark=lig.Generate_New_Structure
			      (core_seed_list[i],frag,tmp_seed_list[k],result);

			 if(mark==FALSE) 
				{
				 // puts("growing fails");
				 continue;
				}
			 else if(mark>(parm->max_population-num_member))
				goto End_of_growing; // too many! stop growing

			 // add the new molecules

			 count=0;	

			 for(k=0;k<24;k++)
			{
			 if(result[k].valid==0) continue;
			 else if(Duplicate_Check(result[k])==TRUE) continue;
			 else
				{
			 	 member[num_member]=result[k];
			 	 sprintf(member[num_member].name,"No_%d_%d", 
				 	 num_generation, num_member+1);
			 	 member[num_member].id=num_member+1;
			 	 num_member++; count++;
				}
			}

			 num_new+=count;
			}
		}
	 if(num_new<num_core_seed) break; 
	}

	End_of_growing:

	delete [] core_seed_list;
	delete [] frag_seed_list;
	delete [] tmp_seed_list;
	delete [] frag_list;
	delete [] growing_record;
	
	return;
}

int Population::Growing_Record_Check(Growing_Record &record) const
{
	int i;

	for(i=0;i<num_record;i++)
	{
 	 if(record.ligand_id!=growing_record[i].ligand_id) continue; 
	 else if(record.lig_seed!=growing_record[i].lig_seed) continue; 
	 else if(record.fragment_id!=growing_record[i].fragment_id) continue;
	 else if(record.frag_seed!=growing_record[i].frag_seed) continue;
	 else return TRUE; 
	}

	return FALSE;
}

void Population::Growing_On_Population(int num_generation)
{
	extern Parameter *parm;
	extern FragLibrary *fraglib;
	int i,j,k,tmp,mark,count;
	int max_record,core_seed,num_frag_seed,num_frag,num_new;
	Ligand lig,frag,result[24];
	int *frag_seed_list,*frag_list,*tmp_seed_list;
	Growing_Record tmp_record;

	if(num_parent==0) return;

	frag_seed_list=new int[500];
	if(frag_seed_list==NULL) Memory_Allocation_Error();

	tmp_seed_list=new int[500];
	if(tmp_seed_list==NULL) Memory_Allocation_Error();

	frag_list=new int[fraglib->num_frag];
	if(frag_list==NULL) Memory_Allocation_Error();

	// growing on each parent molecule in the mating pool.

        max_record=parm->max_population*fraglib->num_frag;

        growing_record=new Growing_Record[max_record];
        if(growing_record==NULL) Memory_Allocation_Error();

	num_record=0;

	while(1)
	{
	 num_new=0;

         num_frag=(parm->max_population-num_member)/num_parent;
	 if(num_frag<1) num_frag=1;

	 for(i=0;i<num_parent;i++)
		{
		 lig=parent[i];

	 	 // first, get the fragment list 

                 num_frag=fraglib->Make_Fragment_List(num_frag,frag_list);

		 // applying each fragment in the list

		 for(j=0;j<num_frag;j++)
			{
		 	 core_seed=parent[i].Select_A_Seed_Hydrogen();
		 	 if(core_seed==FALSE) continue;

			 mark=fraglib->Get_A_Fragment(frag_list[j],frag);
			 if(mark==FALSE) continue; // get fragment failed 

			 num_frag_seed=frag.Make_Seed_List(frag_seed_list);

			 // check whether the seed H has been tried before

			 for(k=0;k<num_frag_seed;k++)
			{
			 tmp_record.ligand_id=i+1;
			 tmp_record.lig_seed=core_seed;
			 tmp_record.fragment_id=frag_list[j];
			 tmp_record.frag_seed=frag_seed_list[k];

			 mark=Growing_Record_Check(tmp_record);
			 if(mark==FALSE) continue;
			 else frag_seed_list[k]=0; 
			}

			 // make a new list containing the available seed Hs

			 tmp=0;		 
			 for(k=0;k<num_frag_seed;k++)
			{
			 if(frag_seed_list[k]==0) continue;
			 else {tmp_seed_list[tmp]=frag_seed_list[k]; tmp++;}
			}

			// select a seed H on the fragment and growing

			 if(tmp==0) continue; // no available seed 
			 else k=(int)(drand48()*tmp);

			 tmp_record.ligand_id=i+1;
			 tmp_record.lig_seed=core_seed;
			 tmp_record.fragment_id=frag_list[j];
			 tmp_record.frag_seed=tmp_seed_list[k];	

			 growing_record[num_record]=tmp_record;
			 if(num_record<(max_record-1)) num_record++;

			 mark=lig.Generate_New_Structure
			      (core_seed,frag,tmp_seed_list[k],result);

			 if(mark==FALSE) continue; // growing fails 
			 else if(mark>(parm->max_population-num_member))
				goto End_of_growing; // too many! stop growing

			 // check growing records and add the new molecules

			 count=0;	

			 for(k=0;k<24;k++)
			{
			 if(result[k].valid==0) continue;
			 member[num_member]=result[k];
			 sprintf(member[num_member].name,"No_%d_%d", 
				 num_generation, num_member+1);
			 num_member++; count++;
			}

			 num_new+=count;
			}
		}

	 if(num_new<num_parent) break; 
	}

	End_of_growing:

	delete [] frag_seed_list;
	delete [] tmp_seed_list;
	delete [] frag_list;
	delete [] growing_record;
	
	return;
}

int Population::Duplicate_Check(Ligand &lig) const
{
	int i,mark;

	mark=FALSE;

	for(i=0;i<num_member;i++)
		{
		 if(member[i].valid==0) continue;
		 else if(Molecule_Duplicate_Check(lig,member[i])==TRUE)
			{
			 mark=TRUE; break;
			}
		 else continue;
		}

	return mark;
}