SciCodePile/SciCode-Domain-Code

1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Allele","Shuhua-Group/ArchaicSeeker2.0","examples.sh",".sh","162","2","./ArchaicSeeker2 -v examples/vcf.par -r examples/remap.par -m examples/model.txt -X examples/outgroup.par -p examples/pop.par -A examples/anc.par -o examples/Han3","Shell"
4"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/104.ArchaicAllelePro/run.sh",".sh","151","2","../000.scripts/archaicAllelePro ../102.ArchaicConsist/Han.Neanderthal.chr22.txt.gz ../103.AF.anc/YRI.anc.txt.gz Han.Neanderthal.archaic.allele.pro.txt5","Shell"
6"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/103.AF.anc/run.sh",".sh","193","5","zcat ../../examples/YRI.chr22.vcf.gz | ../000.scripts/getAF | gzip > YRI.af.txt.gz7 8../000.scripts/annoAnc YRI.af.txt.gz ../../examples/homo_sapiens_ancestor 22 22 YRI.anc.txt9gzip YRI.anc.txt10","Shell"
11"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/101.Seg/run.sh",".sh","165","3","../000.scripts/getSeg ind.txt Han_chr22.seg chr2212cat chr22_Neanderthal.txt | ../000.scripts/segFiltMerge | ../000.scripts/con2gvcf | bgzip > Han_Neanderthal.vcf.gz13","Shell"
14"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/102.ArchaicConsist/run.sh",".sh","165","3","../000.scripts/getArchaicCon ../../examples/Han.chr22.vcf.gz ../101.Seg/Han_Neanderthal.vcf.gz pops.txt Han Han.Neanderthal.chr22.txt15gzip Han.Neanderthal.chr22.txt16","Shell"
17"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/getSeg.cpp",".cpp","2554","106","# include <iostream>18# include <fstream>19# include <vector>20# include <sstream>21# include <map>22# include <set>23 24using namespace std;25 26class bed27{28	public:29		long start, end;30		double gstart, gend;31};32 33int main(int argc, char **argv)34{35	ifstream fpind(argv[1]);36	map<string, int> idCheck;37	vector<string> ids;38	string tid, tpop;39	while(fpind >> tid)40	{41		idCheck[tid + ""_1""] = ids.size();42		ids.push_back(tid + ""_1"");43		idCheck[tid + ""_2""] = ids.size();44		ids.push_back(tid + ""_2"");45	}46	int nind = ids.size();47	vector<map<string,vector<bed> > >data;48	data.resize(nind);49	ifstream fpi(argv[2]);50	string line;51	getline(fpi,line);52	set<string> labs;53	map<string, set<long> > pos;54	string chr;55	map<long,double> gmap;56	while(getline(fpi,line))57	{58		istringstream cl(line);59		bed t;60		string lab;61		cl >> tid >> chr >> t.start >> t.end >> t.gstart >> t.gend >> lab;62		if(lab == ""YRI_Modern_Denisova_Neanderthal"" || lab == ""Denisova_Neanderthal"")63			lab = ""Others"";64		else if(lab != ""Denisova"" && lab != ""Neanderthal"")65			continue;66		if(idCheck.count(tid) == 0)67			continue;68		labs.insert(lab);69		pos[lab].insert(t.start);70		pos[lab].insert(t.end);71		gmap[t.start] = t.gstart;72		gmap[t.end] = t.gend;73		data[idCheck[tid]][lab].push_back(t);74	}75	for(set<string>::iterator it = labs.begin(); it != labs.end(); ++it)76	{77		if(pos.count(*it) == 0)78			continue;79		vector<long> curPos(pos[*it].begin(), pos[*it].end());80		int npos = curPos.size();81		vector<vector<short> > mat;82		mat.resize(npos - 1);83		for(int i = 0 ; i < npos - 1; ++i)84			mat[i].resize(nind, 0);85		for(int i = 0 ; i < nind; ++i)86		{87			if(data[i].count(*it) == 0)88				continue;89			vector<bed>& curBed(data[i][*it]);90			int p(0);91			int nbed(curBed.size());92			for(int j = 0 ; j < npos - 1; ++j)93			{94				if(p == nbed)95					break;96				long start = curPos[j];97				long end = curPos[j + 1];98				if(curBed[p].start > end)99					continue;100				else if(curBed[p].start <= start && curBed[p].end >=end)101					mat[j][i] = 1;102				if(curBed[p].end == end)103					++p;104			}105		}106		string out(argv[3]);107		string path = out + ""_"" + *it + "".txt"";108		ofstream fpo(path.c_str());109		fpo << ""Contig\tStart(bp)\tEnd(bp)\tStart(cM)\tEnd(cM)"";110		for(int i = 0 ; i < nind ; i += 2)111			fpo << '\t' << ids[i].substr(0,ids[i].size() - 2);112		fpo << endl;113		for(int i = 0 ; i < npos - 1; ++i)114		{115			fpo << chr << '\t' << curPos[i] << '\t' << curPos[i + 1] << '\t' << gmap[curPos[i]] << '\t' << gmap[curPos[i+1]];116			for(int j = 0 ; j < nind ; j += 2)117				fpo << '\t' << mat[i][j] << ""|"" << mat[i][j + 1];118			fpo << endl;119		}120	}121}122","C++"
123"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/segFiltMerge.cpp",".cpp","417","27","# include <iostream>124# include <fstream>125# include <vector>126# include <sstream>127 128using namespace std;129 130int main()131{132	string line, tmp;133	getline(cin,line);134	cout << line << endl;135	while(getline(cin,line))136	{137		if(line[0] == 'C')138			continue;139		istringstream cl(line);140		cl >> tmp >> tmp >> tmp >> tmp >> tmp;141		while(cl >> tmp)142			if(tmp[0] == '1' || tmp[2] == '1')143			{144				cout << line << endl;145				break;146			}147	}148}149","C++"
150"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/archaicAllelePro.cpp",".cpp","2133","90","# include <iostream>151# include <fstream>152# include <vector>153# include <sstream>154# include <zlib.h>155 156using namespace std;157 158class snp159{160	public:161		long pos;162		char ref, alt, anc;163		int nafr, sumafr;164};165 166const int MAX = 4096 * 16;167 168double power(0.95);169 170int main(int argc, char **argv)171{172	gzFile fdt, fdafr;173	fdt = gzopen(argv[1], ""rb"");174	fdafr = gzopen(argv[2], ""rb"");175	ofstream fpo(argv[3]);176	vector<vector<snp> > data;177	data.resize(23);178	char buff[MAX];179	gzgets(fdafr, buff, MAX);180	snp tsnp;181	while(gzgets(fdafr, buff, MAX))182	{183		istringstream cl(buff);184		int chr;185		cl >> chr >> tsnp.pos >> tsnp.ref >> tsnp.alt >> tsnp.anc >> tsnp.nafr >> tsnp.sumafr;186		data[chr].push_back(tsnp);187	}188	fpo << ""CHR\tPOS\tREF\tALT\tANC\tAncConfid\tAfrDerFreq\tArchaicDerFreq\tModDerFreq\tArchaicAlleleProb"" << endl;189	gzgets(fdt, buff, MAX);190	vector<int> count;191	count.resize(23, 0);192	while(gzgets(fdt, buff, MAX))193	{194		istringstream cl(buff);195		int chr;196		long pos;197		char ref, alt;198		int nas, sumas, nmod, summod, sum;199		cl >> chr >> pos >> ref >> alt >> nas >> sumas >> nmod >> summod >> sum;200		if(sumas == 0 || nas == 0)201			continue;202		while(count[chr] < data[chr].size() && data[chr][count[chr]].pos < pos)203			++count[chr];204		if(data[chr][count[chr]].pos > pos || count[chr] == data[chr].size())205			continue;206		if(data[chr][count[chr]].pos == pos)207		{208			snp& cur(data[chr][count[chr]]);209			int state(0);210			if(cur.ref == ref && cur.alt == alt)211			{212				char anc(cur.anc);213				anc &= ~' ';214				if(anc == cur.ref)215					state = 0;216				else if(anc == cur.alt)217					state = 1;218				else219					state = -9;220				fpo << chr << '\t' << pos << '\t' << ref << '\t' << alt << '\t' << cur.anc << '\t' << state;221				double afrf, archf, modf, pro;222				afrf = (double) cur.nafr / cur.sumafr;223				archf = (double) nas / sumas;224				modf = (double) nmod / summod;225				if(state == 1)226				{227					afrf = 1 - afrf;228					archf = 1 - archf;229					modf = 1 - modf;230				}231				if(nmod == 0)232					modf = 0;233				pro = ( modf * ( 1 - power ) + 1 - modf ) * ( afrf * 0.003 + 1 - afrf );234				fpo << '\t' << afrf << '\t' << archf << '\t' << modf << '\t' << pro << endl;235			}236		}237	}238}239","C++"
240"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/con2gvcf.cpp",".cpp","708","32","# include <iostream>241# include <vector>242# include <fstream>243# include <sstream>244 245using namespace std;246 247int main()248{249	cout << ""#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT"";250	string line;251	getline(cin,line);252	istringstream head(line);253	head >> line >> line >> line >> line >> line;254	int nind(0);255	while(head >> line)256		cout << '\t' << line;257	cout << endl;258	while(getline(cin,line))259	{260		istringstream cl(line);261		string chr;262		long start, end;263		double gstart, gend;264		cl >> chr >> start >> end >> gstart >> gend;265		cout << chr << '\t' << start << ""\t.\t.\t.\t.\tGenetDis:"" << gstart << ""-"" << gend <<  ""\tEND="" << end - 1 << ""\tGT"";266		while(cl >> line)267			cout << '\t' << line;268		cout << endl;269	}270}271","C++"
272"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/annoAnc.cpp",".cpp","1248","54","# include <iostream>273# include <fstream>274# include <vector>275# include <sstream>276# include <zlib.h>277# include <string>278 279using namespace std;280 281const int MAX = 4096 * 16;282 283int main(int argc, char **argv)284{285	gzFile fdaf;286	fdaf = gzopen(argv[1], ""rb"");287	string pre(argv[2]);288	int chr_start = stoi(argv[3]);289	int chr_end = stoi(argv[4]);290	ofstream fpo(argv[5]);291	vector<vector<string> > anc;292	anc.resize(23);293	char buff[MAX];294	for(int i = chr_start; i <= chr_end; ++i)295	{296		ostringstream path;297		path << pre << ""_"" << i << "".fa.gz"";298		gzFile fdanc;299		fdanc = gzopen(path.str().c_str(), ""rb"");300		gzgets(fdanc, buff, MAX);301		while(gzgets(fdanc, buff, MAX))302		{303			string line(buff);304			anc[i].push_back(line);305//			cout << line.size() << endl;306		}307		gzclose(fdanc);308	}309	int n(anc[chr_start].front().size() - 1);310	fpo << ""CHR\tPOS\tREF\tALT\tANC\tN\tSum"" << endl;311	gzgets(fdaf, buff, MAX);312	while(gzgets(fdaf, buff, MAX))313	{314		istringstream cl(buff);315		int chr;316		long pos;317		string ref, alt, num, sum;318		cl >> chr >> pos >> ref >> alt >> num >> sum;319		char a(anc[chr][(pos - 1) / n][(pos - 1) % n]);320		//char b(anc[chr][pos / n][pos % n]);321 322		fpo << chr << '\t' << pos << '\t' << ref << '\t' << alt << '\t' << a << '\t' << num << '\t' << sum << endl;323	}324}325","C++"
326"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/getArchaicCon.cpp",".cpp","3295","181","/*327 * getArchaicMK.cpp328 *329 *  Created on: Mar 9, 2020330 *      Author: kyuan331 */332 333# include <iostream>334# include <fstream>335# include <sstream>336# include <map>337# include <set>338# include <vector>339# include <cstdlib>340# include <zlib.h>341# include ""boost/dynamic_bitset.hpp""342 343using namespace std;344 345class reg346{347public:348	int chr;349	long start, end;350	boost::dynamic_bitset<> snp;351};352 353class pos354{355public:356	int chr;357	long p;358	string ref, alt;359	boost::dynamic_bitset<> snp;360};361 362//-1 pos < reg; 0 pos -> reg; 1 pos > reg363int cmp(const pos& p, const reg& r)364{365	if(p.chr < r.chr)366		return -1;367	if(p.chr > r.chr)368		return 1;369	if(p.p < r.start)370		return -1;371	if(p.p > r.end)372		return 1;373	return 0;374}375 376int cmp(const pos& p1, const pos& p2)377{378	if(p1.chr < p2.chr)379		return -1;380	if(p1.chr > p2.chr)381		return 1;382	if(p1.p < p2.p)383		return -1;384	if(p1.p > p2.p)385		return 1;386	return 0;387}388 389const int MAX = 4096 * 16;390 391int main(int argc, char **argv)392{393	gzFile fptest, fpas;394	fptest = gzopen(argv[1], ""rb"");395	fpas = gzopen(argv[2], ""rb"");396	ifstream fpind(argv[3]);397	string pop(argv[4]);398	ofstream fpo(argv[5]);399 400	set<string> ids;401	string id, p;402	while(fpind >> id >> p)403		if(p == pop)404			ids.insert(id);405 406	map<string, int> check;407	vector<int> sel;408	char buff[MAX];409	gzgets(fpas, buff, MAX);410	istringstream ashead(buff);411	string line;412	for(int i = 0 ; i < 9 ; ++i)413		ashead >> line;414	int n(0);415	while(ashead >> line)416		if(ids.count(line))417		{418			check[line] = n;419			++n;420			sel.push_back(1);421		}422		else423			sel.push_back(0);424	vector<reg> asReg;425	reg treg;426	treg.snp.resize(n * 2);427	while(gzgets(fpas, buff, MAX))428	{429		istringstream cl(buff);430		treg.snp.reset();431		cl >> treg.chr >> treg.start >> line >> line >> line >> line >> line >> line;432		treg.end = atol(line.substr(4).c_str());433		cl >> line;434		int m(0);435		for(int i = 0 ; i < sel.size() ; ++i)436		{437			cl >> line;438			if(sel[i])439			{440				if(line[0] == '1')441					treg.snp.set(m * 2);442				if(line[2] == '1')443					treg.snp.set(m * 2 + 1);444				++m;445			}446		}447		asReg.push_back(treg);448	}449 450	while(gzgets(fptest, buff, MAX))451		if(buff[1] != '#')452			break;453	istringstream head(buff);454	for(int i = 0 ; i < 9 ; ++i)455		head >> line;456	vector<int> index;457	int count(0);458	while(head >> id)459	{460		if(check.count(id))461			index.push_back(check[id]);462		else463			index.push_back(-1);464	}465	pos tpos;466	tpos.snp.resize(n * 2);467	int s(0);468	fpo << ""CHR\tPos\tRef\tAlt\tAltArchaic\tNArchaic\tAltMod\tNMod\tSum"" << endl;469	while(gzgets(fptest, buff, MAX))470	{471		istringstream cl(buff);472		cl >> tpos.chr >> tpos.p >> line >> tpos.ref >> tpos.alt >> line >> line >> line >> line;473		tpos.snp.reset();474		for(int i = 0 ; i < index.size(); ++i)475		{476			cl >> line;477			if(index[i] >= 0)478			{479				if(line[0] == '1')480					tpos.snp.set(index[i] * 2);481				if(line[2] == '1')482					tpos.snp.set(index[i] * 2 + 1);483			}484		}485		int c(cmp(tpos, asReg[s]));486		if(c == -1)487			continue;488		while(c == 1)489		{490			++s;491			if(s == asReg.size())492				break;493			if(cmp(tpos, asReg[s]) <= 0)494				break;495		}496		if(s == asReg.size())497			break;498		if(c == 0)499		{500				fpo << tpos.chr<<'\t' << tpos.p << '\t' << tpos.ref << '\t' << tpos.alt << '\t' << (tpos.snp & asReg[s].snp).count() << '\t' << asReg[s].snp.count() << '\t' << (tpos.snp & (~asReg[s].snp) ).count() << '\t' << (~asReg[s].snp).count() << '\t' << n * 2 << endl;501		}502	}503}504 505 506","C++"
507"Allele","Shuhua-Group/ArchaicSeeker2.0","ArchaicAllele/000.scripts/getAF.cpp",".cpp","1365","79","# include <iostream>508# include <fstream>509# include <map>510# include <sstream>511# include <set>512# include <vector>513 514using namespace std;515 516int main(int argc, char **argv)517{518	bool all(false);519	if(argc >= 2)520		all = true;521	ifstream fpind;522	if(!all)523		fpind.open(argv[1]);524	set<string> ids;525	string id;526	vector<int> check;527	string line;528	cout << ""#CHR\tPOS\tREF\tALT\tN\tSum"" << endl;529	while(getline(cin, line))530		if(line[1] != '#')531			break;532	istringstream head(line);533 534	if(all)535	{536		while(fpind >> id)537			ids.insert(id);538		for(int i = 0 ; i < 9 ; ++i)539			head >> line;540		while(head >> id)541			if(ids.count(id))542				check.push_back(1);543			else544				check.push_back(0);545	}546	else547	{548		for(int i = 0 ; i < 9 ; ++i)549			head >> line;550		while(head >> id)551			check.push_back(1);552	}553	while(getline(cin, line))554	{555		istringstream cl(line);556		if(line[0] == '#')557			continue;558		string chr, pos, ref, alt;559		cl >> chr >> pos >> line >> ref >> alt >> line >> line >> line >> line;560		int n(0), sum(0);561		for(int i = 0 ; i < check.size(); ++i)562		{563			cl >> line;564			if(check[i])565			{566				if(line[0] == '1')567				{568					++n;569					++sum;570				}571				else if(line[0] == '0')572					++sum;573				if(line[2] == '1')574				{575					++n;576					++sum;577				}578				else if(line[2] == '0')579					++sum;580			}581		}582		cout << chr << '\t' << pos << '\t' << ref << '\t' << alt << '\t' << n << '\t' << sum << endl;583	}584}585","C++"
586"Allele","Shuhua-Group/ArchaicSeeker2.0","WaveEstimate/getAS2Seg/getAS2Seg.cpp",".cpp","3358","148","# include <iostream>587# include <vector>588# include <sstream>589# include <fstream>590# include <map>591# include <cstring>592# include <iomanip>593# include <set>594# include <zlib.h>595 596using namespace std;597 598const int MAX = 4096 * 256;599 600class gmap601{602	public:603		string gmapPath;604		gmap(string gmapPath): gmapPath(gmapPath)605		{606			pos.resize(23);607			gdis.resize(23);608			mem.resize(23);609			for(int i = 1 ; i <= 22 ; ++i)610			{611				ostringstream path;612				path << gmapPath << ""/genetic_map_chr"" << i << ""_combined_b37.txt"";613				ifstream fpi(path.str().c_str());614				vector<long> &curPos(pos[i]);615				vector<double> &curGdis(gdis[i]);616				curPos.clear();617				curGdis.clear();618				string line;619				fpi >> line >> line >> line;620				long tpos;621				double tgdis;622				while(fpi >> tpos >> line >> tgdis)623				{624					curPos.push_back(tpos);625					curGdis.push_back(tgdis);626				}627			}628		}629		double search(int chr, long p)630		{631			if(mem[chr].count(p))632				return mem[chr][p];633			const vector<long>& curPos(pos[chr]);634			const vector<double>& curGdis(gdis[chr]);635			int start(0), end(curPos.size() - 1);636			int mid;637			while(end - start > 1)638			{639				mid = (end + start) / 2;640				if(curPos[mid] > p)641					end = mid;642				else if(curPos[mid] < p)643					start = mid;644				else645				{646					mem[chr][p] = curGdis[mid];647					return curGdis[mid];648				}649			}650			while(start < curPos.size() && curPos[start] < p)651				++start;652			if(start == curPos.size())653			{654				mem[chr][p] = curGdis.back();655				return curGdis.back();656			}657			else if(start <= 1)658			{659				mem[chr][p] = curGdis.front();660				return curGdis.front();661			}662			--start;663			double r = curGdis[start] + ( p - curPos[start]	) * ( curGdis[start + 1] - curGdis[start] ) / ( curPos[start + 1] - curPos[start] );664			mem[chr][p] = r;665			return r;666		}667		vector<map<long,double> > mem;668		vector<vector<long> > pos;669		vector<vector<double> > gdis;670};671 672class bed673{674	public:675		int chr;676		long start, end;677		double gstart, gend;678		string id;679};680 681int main(int argc, char **argv)682{683	string gmapPath(argv[1]);684	gmap g(gmapPath);685	ifstream fpi(argv[2]);686	string out(argv[3]);687	set<string> check;688	for(int i = 3; i < argc ; ++i)689		check.insert(argv[i]);690	string line;691	getline(fpi,line);692	map<string, vector<vector<bed> > > segs;693	while(getline(fpi, line))694	{695		istringstream cl(line);696		string lab;697		bed tbed;698		cl >> tbed.id >> tbed.chr >> tbed.start >> tbed.end >> line >> line >> lab;699		if(!check.count(lab))700			continue;701		tbed.gstart = g.search(tbed.chr, tbed.start);702		tbed.gend = g.search(tbed.chr, tbed.end);703		segs[tbed.id].resize(23);704		segs[tbed.id][tbed.chr].push_back(tbed);705	}706	string path;707	path = out + "".seg"";708	ofstream fpo(path.c_str());709	fpo << setprecision(20);710 711	for(map<string, vector<vector<bed> > >::iterator it = segs.begin(); it != segs.end(); ++it)712	{713		for(int i = 1 ; i <= 22 ; ++i)714		{715			const vector<bed> & curData(it->second.at(i));716			if(curData.size() ==0)717				continue;718			double gstart(g.gdis[i].front());719			long start(g.pos[i].front());720			for(int j = 0 ; j < curData.size(); ++j)721			{722				fpo << gstart / 100 << '\t' << curData[j].gstart / 100 << ""\tModern"" << endl;723				fpo << curData[j].gstart / 100 << '\t' << curData[j].gend / 100 << ""\tArchaic"" << endl;724				gstart = curData[j].gend;725				start = curData[j].end;726			}727			fpo << gstart / 100 << '\t' << g.gdis[i].back() / 100 << ""\tModern"" << endl;728		}729	}730}731 732 733","C++"
734"Allele","Shuhua-Group/ArchaicSeeker2.0","WaveEstimate/MultiWaver2.1/ParamExp.hpp",".hpp","2777","124","/*735 * ParamExp.hpp736 *737 *  Created on: May 26, 2015738 *  Author: young739 */740/******************************************************************************741 * @brief Parameters of EM742 * Class ParamExp to store the parameters of EM743 * Here assume the observations are from the combination of K sets of numbers744 * which following exponential distributions, with parameters lambda[j], and the745 * mix proportions prop[j], the purpose of EM is used to estimate the parameters746 * lambda[j] and mix proportions[j], where j in [1,2,...,K]747 *****************************************************************************/748 749#ifndef PARAMEXP_HPP_750#define PARAMEXP_HPP_751 752//const double kDelta = 0.000001; //converge condition753 754# include <vector>755 756class ParamExp757{758public:759 760	/*761	 * @brief Constructor with given value K762	 * @param K number of exponential distribution763	 */764	ParamExp(int K = 1);765 766	ParamExp(int K, std::vector<double> observ);767 768    /*769	 * @brief Constructor with given value K, lambdas and proportions770	 * @param K number of exponential distribution771	 * @param lambda initial value of lambda772	 * @param prop initial value of proportions773	 */774	ParamExp(int K, double *lambda, double *prop);775 776    /*777	 * @brief Copy constructor778	 * @param rhs	old parameter779	 */780	ParamExp(const ParamExp &rhs);781 782    /*783	 * @brief overloading of operator =784	 * @param rhs	old parameter785	 */786	ParamExp & operator=(const ParamExp &rhs);787 788    /*789	 * @brief get the number of exponential distributions790	 * @return integer number791	 */792	int getK() const;793 794    /*795	 * @brief get the value of lambda with given index796	 * @param index the exponential distribution797	 * @return the value of lambda798	 */799	double getLambda(int index) const;800 801    /*802	 * @brief get the value of proportion with given index803	 * @param index the exponential distribution804	 * @return the value of proportion805	 */806	double getProp(int index) const;807//	void setK(int);808//	void setLambda(int index, double lambda);809 810    /*811	 * set the Proportion at index812	 */813	void setProp(int index, double prop);814 815    /*816	 * @brief check if the parameter is converged817	 * @param par old parameter818	 * @param epsilon epsilon to check whether a parameter converge or not819	 * @return true if is converged820	 */821	bool isConverge(const ParamExp &par, double epsilon = 0.000001);822 823    /*824	 * @brief sort the parameter by lambda in ascend order825	 */826	void sortByLambda();827 828    /*829	 * @brief print the value of parameters830	 */831	void print();832 833    /*834	 * @brief destructor835	 */836	virtual ~ParamExp();837 838private:839 840    /*841	 * number of exponential distributions842	 */843	int K;844 845    /*846	 * parameter lambdas for exponential distributions847	 */848	double *lambda;849 850    /*851	 * mix proportions for exponential distributions852	 */853	double *prop;854};855 856#endif /* PARAMEXP_HPP_ */857","Unknown"
858"Allele","Shuhua-Group/ArchaicSeeker2.0","WaveEstimate/MultiWaver2.1/EMExp.hpp",".hpp","2944","119","/*859 * EMExp.hpp860 *861 *  Created on: May 26, 2015862 *      Author: young863 */864/******************************************************************************865 * @brief A class perform EM866 * Class EMExp to perform EM algorithm in parameters estimation with hidden variable867 * The EMExp class has three attributes:868 * 1) parameters to be estimated;869 * 2) a sequence of observations, which used to estimate the parameters;870 * 3) and the likelihood corresponding to the observations and parameters.871 872 * The method iterate is used to perform EM algorithm, including two steps:873 * E-Step:874 * 1) calculate p_j=prob(z_i=j|x_i, theta_t)875 * =prob(x_i|z_i=j,theta_t)*prob(z_i=j,theta_t)/(sum_j{from 1 to K}(.)876 * denotes the numerator877 * 2) calculate p_j*xi ; j=1,...,K878 * theta_t: parameters at the t-th iteration879 * M-Step:880 * 1) update prop, refer as m, m_j=sum(p_j)/sum(p_1+p_2+...+p_K)=sum(p_j)/n881 * 2) update lambda, lambda_j=sum(p_j)/sum(p_j*x_i)882 * note: update the parameters for t+1 times iteration883 ******************************************************************************/884#ifndef EMEXP_HPP_885#define EMEXP_HPP_886 887#include <vector>888#include ""ParamExp.hpp""889 890/*891 * summation over data, with certain size892 */893double sum(double *data, int size);894 895double modernWaver(unsigned int n, const double *x, double *grad, void *data);896 897double archaicWaver(unsigned int n, const double *x, double *grad, void *data);898 899double archaicWaverConstraint(unsigned n, const double *x, double *grad, void *data);900 901class waverData902{903public:904	waverData(const std::vector<double>& _observ, double _lb, double _ub) : \905			lb(_lb), ub(_ub), observ(_observ) {};906	double lb, ub;907	const std::vector<double>& observ;908};909 910class EMExp911{912public:913 914	/*915	 * @brief constructor916	 * @param par Parameter for EM algorithm917	 * @param observ Observations of data used to estimate parameter of EM918	 */919	EMExp(const ParamExp &par, const std::vector<double> &observ);920 921    /*922	 * @brief get the log likelihood923	 * @return the log likelihood of current data and parameter924	 */925	double getLik() const;926 927    /*928	 * @brief get parameter of EM929	 * @return the reference of EM parameter930	 */931	ParamExp & getPar();932 933    /*934	 * @brief set a new parameter for EM935	 * @param new parameter to be set936	 */937	void setPar(const ParamExp &par);938 939    /*940	 * @brief update the log likelihood, usually when parameter is updated941	 */942	void updateLik();943 944    /*945	 * @brief perform EM iteration946	 * @param maxIter max number of iteration, the EM iteration terminate either947	 * parameter is converged, or reach the max number of iteration948	 */949	void iterate(int maxIter, double epsilon);950 951	void iterateClassify(int maxIter, double epsilon);952 953    /*954	 * destructor955	 */956	virtual ~EMExp();957 958private:959    /*960     * log likelihood value961     */962	double lik;963 964    /*965     * Parameter of EM966     */967	ParamExp par;968 969    /*970     * a vector of observations971     */972	std::vector<double> observ;973};974 975#endif /* EMEXP_HPP_ */976","Unknown"
977"Allele","Shuhua-Group/ArchaicSeeker2.0","WaveEstimate/MultiWaver2.1/ParamExp.cpp",".cpp","3330","191","/*978 * ParamExp.cpp979 *980 *  Created on: May 26, 2015981 *  Author: young982 */983#include <cmath>984#include <cstdlib>985#include <map>986#include <iostream>987#include <algorithm>988 989#include ""ParamExp.hpp""990 991using namespace std;992 993/* constructor994 * Initial with K exponential distributions995 * parameters lambda are randomly initial with value between 0 and 1996 * proportions are randomly initial with value between 0 and 1, and997 * the summation equals 1998 */999ParamExp::ParamExp(int K) :1000		K(K)1001{1002	lambda = new double[K];1003	prop = new double[K];1004	for (int i = 0; i < K; ++i)1005	{1006		lambda[i] = 1.0 * rand() / RAND_MAX;1007		prop[i] = 1.0 / K;1008	}1009 1010//	double tmp = 0;1011//	for (int i = 0; i < (K - 1); ++i)1012//	{1013//		prop[i] = rand() / (1.0 * K * RAND_MAX);1014//		tmp += prop[i];1015//	}1016//	//ensure sum to one1017//	prop[K - 1] = 1 - tmp;1018}1019 1020ParamExp::ParamExp(int K, std::vector<double> observ) : K(K)1021{1022	sort(observ.begin(), observ.end());1023	lambda = new double[K];1024	prop = new double[K];1025	if(K == 1)1026	{1027		lambda[0] = 1.0 * rand() / RAND_MAX;1028		prop[0] = 1.0;1029	}1030	else1031	{1032		int nobserv = observ.size();1033		for(int i = 0 ; i < K ; ++i)1034		{1035			int index = i * (double) nobserv / (K - 1);1036			if(index >= nobserv)1037				index = nobserv - 1;1038			lambda[i] = 1.0 / observ[index];1039			prop[i] = 1.0 / K;1040		}1041	}1042}1043 1044/* another constructor1045 * initialize with K, lambdas and proportions1046 */1047ParamExp::ParamExp(int K, double *l, double *p) :1048		K(K)1049{1050	lambda = new double[K];1051	prop = new double[K];1052	for (int i = 0; i < K; ++i)1053	{1054		lambda[i] = l[i];1055		prop[i] = p[i];1056	}1057}1058 1059//copy constructor1060ParamExp::ParamExp(const ParamExp &rhs)1061{1062	K = rhs.K;1063	lambda = new double[K];1064	prop = new double[K];1065	for (int i = 0; i < K; ++i)1066	{1067		lambda[i] = rhs.lambda[i];1068		prop[i] = rhs.prop[i];1069	}1070}1071 1072//overloading operator assignment1073ParamExp & ParamExp::operator=(const ParamExp &rhs)1074{1075	if (this != &rhs)1076	{1077		K = rhs.K;1078		if (lambda != NULL)1079			delete[] lambda;1080		if (prop != NULL)1081			delete[] prop;1082		lambda = new double[K];1083		prop = new double[K];1084		for (int i = 0; i < K; ++i)1085		{1086			lambda[i] = rhs.lambda[i];1087			prop[i] = rhs.prop[i];1088		}1089	}1090	return *this;1091}1092 1093int ParamExp::getK() const1094{1095	return K;1096}1097 1098double ParamExp::getLambda(int index) const1099{1100	return lambda[index];1101}1102 1103double ParamExp::getProp(int index) const1104{1105	return prop[index];1106}1107 1108void ParamExp::setProp(int index, double prop)1109{1110	this->prop[index] = prop;1111}1112 1113//test convergence1114bool ParamExp::isConverge(const ParamExp &par, double epsilon)1115{1116	bool converge = true;1117	for (int i = 0; i < K; ++i)1118	{1119		if (abs(lambda[i] - par.getLambda(i)) > epsilon || abs(prop[i] - par.getProp(i)) > epsilon)1120		{1121			converge = false;1122			break;1123		}1124	}1125	return converge;1126}1127 1128void ParamExp::sortByLambda()1129{1130	/*1131	 * the key of map are automatically sorted, therefore can be used1132	 * to sort lambda and corresponding proportion accordingly1133	 */1134	map<double, double> temp;1135	int i;1136	for (i = 0; i < K; ++i)1137	{1138		temp[lambda[i]] = prop[i];1139	}1140	i = 0;1141	for (map<double, double>::iterator it = temp.begin(); it != temp.end(); it++)1142	{1143		lambda[i] = it->first;1144		prop[i] = it->second;1145		i++;1146	}1147}1148 1149void ParamExp::print()1150{1151	cout << ""par=("";1152	for (int i = 0; i < K; ++i)1153	{1154		cout << lambda[i] << "", "" << prop[i] << ""; "";1155	}1156	cout << "")"" << endl;1157}1158 1159//1160ParamExp::~ParamExp()1161{1162	if (lambda != NULL)1163		delete[] lambda;1164	if (prop != NULL)1165		delete[] prop;1166}1167","C++"
1168"Allele","Shuhua-Group/ArchaicSeeker2.0","WaveEstimate/MultiWaver2.1/Utils.cpp",".cpp","5187","193","/*1169 * Utils.cpp1170 *1171 *  Created on: Aug 26, 20151172 *      Author: young1173 */1174#include <cmath>1175#include <iostream>1176#include ""EMExp.hpp""1177#include ""Utils.hpp""1178#include <boost/math/distributions/chi_squared.hpp>1179 1180using namespace std;1181 1182double cv_chisq(int df, double alpha)1183{1184	boost::math::chi_squared dist(df);1185	return boost::math::quantile(dist, 1 - alpha);1186}1187 1188ParamExp findOptPar(const vector<double> &observ, int K, int maxIter, double ancestryProp, \1189		double criticalValue, double epsilon, double minP, bool simple)1190{1191//	bool findOpt = false;1192//	int k = 1;1193//	ParamExp parPrev(K);1194//	EMExp em(parPrev, observ);1195//	em.iterateClassify(maxIter, epsilon);1196//	double llkPrev = em.getLik();1197//	parPrev = em.getPar();1198//	if (simple)1199//	{1200//		parPrev.sortByLambda();

Showing the first 1,200 of 156077 lines. Download the file for the rest.