1 #ifndef __PLINK2_COMMON_H__
2 #define __PLINK2_COMMON_H__
3
4 // This library is part of PLINK 2.00, copyright (C) 2005-2020 Shaun Purcell,
5 // Christopher Chang.
6 //
7 // This program is free software: you can redistribute it and/or modify it
8 // under the terms of the GNU Lesser General Public License as published by the
9 // Free Software Foundation, either version 3 of the License, or (at your
10 // option) any later version.
11 //
12 // This library is distributed in the hope that it will be useful, but WITHOUT
13 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
15 // for more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with this program. If not, see <http://www.gnu.org/licenses/>.
19
20
21 // Resources needed across a variety of plink2 modules involving
22 // plink2-specific constructs (e.g. chr_info_t, variable-length sample IDs,
23 // 2-bit genotypes). More generic library code has been moved to plink2_base
24 // and plink2_cmdline.
25
26 #include "include/pgenlib_read.h"
27 #include "plink2_decompress.h"
28
29 #ifdef __cplusplus
30 namespace plink2 {
31 #endif
32
33 #define PROG_NAME_STR "plink2"
34
35 // leave the door semi-open to 32-bit dosages (or 31? 24?)
36 // note that u31tod()/u31tof() can't be used on 32-bit dosages
37 typedef uint16_t Dosage;
38 typedef int16_t SDosage;
39 typedef uint32_t DosageProd;
40 #define kDosageMax (1U << (8 * sizeof(Dosage) - 1))
41 CONSTI32(kDosageMid, kDosageMax / 2);
42 CONSTI32(kDosage4th, kDosageMax / 4);
43 CONSTI32(kDosageMissing, kDosageMax * 2 - 1);
44 static const double kRecipDosageMax = 0.000030517578125;
45 static const double kRecipDosageMid = 0.00006103515625;
46 static const double kRecipDosageMidSq = 0.0000000037252902984619140625;
47 static const float kRecipDosageMidf = 0.00006103515625;
48 CONSTI32(kDosagePerVec, kBytesPerVec / sizeof(Dosage));
49 CONSTI32(kDosagePerCacheline, kCacheline / sizeof(Dosage));
50
DosageCtToVecCt(uintptr_t val)51 HEADER_CINLINE uintptr_t DosageCtToVecCt(uintptr_t val) {
52 return DivUp(val, kDosagePerVec);
53 }
54
DosageCtToAlignedWordCt(uintptr_t val)55 HEADER_CINLINE uintptr_t DosageCtToAlignedWordCt(uintptr_t val) {
56 return kWordsPerVec * DosageCtToVecCt(val);
57 }
58
59 // dosage_int = 0..2 value in 16384ths
60 // returns distance from 0.5 or 1.5 in 16384ths, whichever is closer
BiallelicDosageHalfdist(uint32_t dosage_int)61 HEADER_INLINE uint32_t BiallelicDosageHalfdist(uint32_t dosage_int) {
62 const uint32_t dosage_int_rem = dosage_int & (kDosageMid - 1);
63 return abs_i32(S_CAST(int32_t, dosage_int_rem) - kDosage4th);
64 }
65
HaploidDosageHalfdist(uint32_t dosage_int)66 HEADER_INLINE uint32_t HaploidDosageHalfdist(uint32_t dosage_int) {
67 return abs_i32(S_CAST(int32_t, dosage_int) - kDosage4th);
68 }
69
DosageHetdist(uint32_t dosage_int)70 HEADER_INLINE uint32_t DosageHetdist(uint32_t dosage_int) {
71 return abs_i32(S_CAST(int32_t, dosage_int) - kDosageMid);
72 }
73
DosageHomdist(uint32_t dosage_int)74 HEADER_INLINE uint32_t DosageHomdist(uint32_t dosage_int) {
75 // todo: compare vs. branchless
76 return (dosage_int > kDosageMid)? (kDosageMax - dosage_int) : dosage_int;
77 }
78
79 // on 0..32768 scale
DphaseHalfdist(uint32_t dphase_side)80 HEADER_INLINE uint32_t DphaseHalfdist(uint32_t dphase_side) {
81 return abs_i32(S_CAST(int32_t, dphase_side) - kDosageMid);
82 }
83
84 // this is a bit arbitrary
85 CONSTI32(kMaxPhenoCt, 524287);
86 #define MAX_PHENO_CT_STR "524287"
87
88 // maximum number of usable cluster computers, this is arbitrary though it
89 // shouldn't be larger than 2^32 - 1
90 // (actually, there's an overflow danger: [work units] * parallel_idx may not
91 // fit in a uint64 if parallel_tot is too high.)
92 CONSTI32(kParallelMax, 32768);
93
94 // unnecessary to use e.g. (1LLU << 0), the FLAGSET64 macros should force the
95 // integer type to 64-bit.
96 FLAGSET64_DEF_START()
97 kfMisc0,
98 kfMiscAffection01 = (1 << 0),
99 kfMiscAllowExtraChrs = (1 << 1),
100 kfMiscRealRefAlleles = (1 << 2),
101 kfMiscMajRef = (1 << 3),
102 kfMiscMajRefForce = (1 << 4),
103 kfMiscNonfounders = (1 << 5),
104 kfMiscExcludePvarFilterFail = (1 << 6),
105 kfMiscAutosomePar = (1 << 7),
106 kfMiscAutosomeOnly = (1 << 8),
107 kfMiscMergePar = (1 << 9),
108 kfMiscPhenoColNums = (1 << 10),
109 kfMiscCovarColNums = (1 << 11),
110 kfMiscHweMidp = (1 << 12),
111 kfMiscHweKeepFewhet = (1 << 13),
112 kfMiscWriteSnplistZs = (1 << 14),
113 kfMiscWriteSnplistAllowDups = (1 << 15),
114 kfMiscMafSucc = (1 << 16),
115 kfMiscGenoDosage = (1 << 17),
116 kfMiscGenoHhMissing = (1 << 18),
117 kfMiscMindDosage = (1 << 19),
118 kfMiscMindHhMissing = (1 << 20),
119 kfMiscGenotypingRateDosage = (1 << 21),
120 kfMiscSetMissingVarIds = (1 << 22),
121 kfMiscChrOverrideCmdline = (1 << 23),
122 kfMiscChrOverrideFile = (1 << 24),
123 kfMiscNewVarIdOverflowMissing = (1 << 25),
124 kfMiscNewVarIdOverflowTruncate = (1 << 26),
125 kfMiscRequirePheno = (1 << 27),
126 kfMiscRequireCovar = (1 << 28),
127 kfMiscCatPhenoFamily = (1 << 29),
128 kfMiscRefAlleleForce = (1 << 30),
129 kfMiscAlt1AlleleForce = (1U << 31),
130 kfMiscMergeX = (1LLU << 32),
131 kfMiscNoIdHeader = (1LLU << 33),
132 kfMiscNoIdHeaderIidOnly = (1LLU << 34),
133 kfMiscBiallelicOnly = (1LLU << 35),
134 kfMiscBiallelicOnlyStrict = (1LLU << 36),
135 kfMiscBiallelicOnlyList = (1LLU << 37),
136 kfMiscStrictSid0 = (1LLU << 38),
137 kfMiscAllowBadFreqs = (1LLU << 39),
138 kfMiscIidSid = (1LLU << 40),
139 kfMiscPhenoIidOnly = (1LLU << 41),
140 kfMiscCovarIidOnly = (1LLU << 42),
141 kfMiscAllowBadLd = (1LLU << 43)
142 FLAGSET64_DEF_END(MiscFlags);
143
144 FLAGSET64_DEF_START()
145 kfExportf0,
146 kfExportf01 = (1 << 0),
147 kfExportf12 = (1 << 1),
148 kfExportfSpaces = (1 << 2),
149 kfExportfRefFirst = (1 << 3),
150 kfExportf23 = (1 << 4),
151 kfExportfA = (1 << 5),
152 kfExportfATranspose = (1 << 6),
153 kfExportfAD = (1 << 7),
154 kfExportfBcf42 = (1 << 8),
155 kfExportfBcf43 = (1 << 9),
156 kfExportfBcf = kfExportfBcf42 | kfExportfBcf43,
157 kfExportfBeagle = (1 << 10),
158 kfExportfBeagleNomap = (1 << 11),
159 kfExportfBgen11 = (1 << 12),
160 kfExportfBgen12 = (1 << 13),
161 kfExportfBgen13 = (1 << 14),
162 kfExportfBimbam = (1 << 15),
163 kfExportfBimbam1chr = (1 << 16),
164 kfExportfFastphase = (1 << 17),
165 kfExportfFastphase1chr = (1 << 18),
166 kfExportfHaps = (1 << 19),
167 kfExportfHapsLegend = (1 << 20),
168 kfExportfHv = (1 << 21),
169 kfExportfHv1chr = (1 << 22),
170 kfExportfIndMajorBed = (1 << 23),
171 kfExportfLgen = (1 << 24),
172 kfExportfLgenRef = (1 << 25),
173 kfExportfList = (1 << 26),
174 kfExportfRlist = (1 << 27),
175 kfExportfOxGen = (1 << 28),
176 kfExportfPed = (1 << 29),
177 kfExportfCompound = (1 << 30),
178 kfExportfStructure = (1U << 31),
179 kfExportfTranspose = (1LLU << 32),
180 kfExportfVcf42 = (1LLU << 33),
181 kfExportfVcf43 = (1LLU << 34),
182 kfExportfVcf = kfExportfVcf42 | kfExportfVcf43,
183 kfExportfTypemask = (2LLU * kfExportfVcf43) - kfExportf23,
184 kfExportfIncludeAlt = (1LLU << 35),
185 kfExportfBgz = (1LLU << 36),
186 kfExportfOmitNonmaleY = (1LLU << 37)
187 FLAGSET64_DEF_END(ExportfFlags);
188
189 FLAGSET_DEF_START()
190 kfInfo0,
191 kfInfoPrFlagPresent = (1 << 0),
192 kfInfoPrNonflagPresent = (1 << 1),
193 kfInfoNonprPresent = (1 << 2),
194 kfInfoAll = ((kfInfoNonprPresent * 2) - kfInfoPrFlagPresent),
195 // this is set in the .bim case
196 kfInfoPrNonrefDefault = (1 << 3),
197 FLAGSET_DEF_END(InfoFlags);
198
199 FLAGSET_DEF_START()
200 kfPvarPsam0,
201 kfPvarZs = (1 << 0),
202
203 kfPvarColXheader = (1 << 1),
204 kfPvarColVcfheader = (1 << 2),
205 kfPvarColMaybequal = (1 << 3),
206 kfPvarColQual = (1 << 4),
207 kfPvarColMaybefilter = (1 << 5),
208 kfPvarColFilter = (1 << 6),
209 kfPvarColMaybeinfo = (1 << 7),
210 kfPvarColInfo = (1 << 8),
211 kfPvarColXinfo = (kfPvarColInfo * 2) - kfPvarColMaybeinfo,
212 kfPvarColMaybecm = (1 << 9),
213 kfPvarColCm = (1 << 10),
214 kfPvarColDefault = (kfPvarColXheader | kfPvarColMaybequal | kfPvarColMaybefilter | kfPvarColMaybeinfo | kfPvarColMaybecm),
215 kfPvarColAll = ((kfPvarColCm * 2) - kfPvarColXheader),
216 kfPsamColMaybefid = (1 << 11),
217 kfPsamColFid = (1 << 12),
218 kfPsamColMaybesid = (1 << 13),
219 kfPsamColSid = (1 << 14),
220 kfPsamColMaybeparents = (1 << 15),
221 kfPsamColParents = (1 << 16),
222 kfPsamColSex = (1 << 17),
223 kfPsamColPheno1 = (1 << 18),
224 kfPsamColPhenos = (1 << 19),
225 kfPsamColDefault = (kfPsamColMaybefid | kfPsamColMaybesid | kfPsamColMaybeparents | kfPsamColSex | kfPsamColPhenos),
226 kfPsamColAll = ((kfPsamColPhenos * 2) - kfPsamColMaybefid)
227 FLAGSET_DEF_END(PvarPsamFlags);
228
229 FLAGSET_DEF_START()
230 kfSampleId0,
231 kfSampleIdFidPresent = (1 << 0),
232 kfSampleIdParentsPresent = (1 << 1),
233 kfSampleIdNoIdHeader = (1 << 2),
234 kfSampleIdNoIdHeaderIidOnly = (1 << 3),
235 kfSampleIdStrictSid0 = (1 << 4)
236 FLAGSET_DEF_END(SampleIdFlags);
237
238 // These structs are small enough and ownership of the pointed-to arrays is
239 // generall clear enough that the noncopyable annotation is just intended to be
240 // a tripwire.
241 typedef struct SampleIdInfoStruct {
242 NONCOPYABLE(SampleIdInfoStruct);
243 char* sample_ids;
244 char* sids;
245 uintptr_t max_sample_id_blen;
246 uintptr_t max_sid_blen;
247 SampleIdFlags flags;
248 } SampleIdInfo;
249
250 typedef struct ParentalIdInfoStruct {
251 NONCOPYABLE(ParentalIdInfoStruct);
252 char* paternal_ids;
253 char* maternal_ids;
254 uintptr_t max_paternal_id_blen;
255 uintptr_t max_maternal_id_blen;
256 } ParentalIdInfo;
257
258 typedef struct PedigreeIdInfoStruct {
259 SampleIdInfo sii;
260 ParentalIdInfo parental_id_info;
261 } PedigreeIdInfo;
262
263 typedef struct APermStruct {
264 uint32_t min;
265 uint32_t max;
266 double alpha;
267 double beta;
268 double init_interval;
269 double interval_slope;
270 } APerm;
271
272 // (2^31 - 1000001) / 2
273 CONSTI32(kApermMax, 1073241823);
274
275 typedef struct TwoColParamsStruct {
276 NONCOPYABLE(TwoColParamsStruct);
277 uint32_t colx;
278 uint32_t colid;
279 uint32_t skip_ct;
280 char skipchar;
281 char fname[];
282 } TwoColParams;
283
284 void InitPedigreeIdInfo(MiscFlags misc_flags, PedigreeIdInfo* piip);
285
286 // no CleanupPedigreeIdInfo function since LoadPsam() allocates the arrays in
287 // bigstack.
288
289
bigstack_alloc_ac(uintptr_t ct,AlleleCode ** allele_arr_ptr)290 HEADER_INLINE BoolErr bigstack_alloc_ac(uintptr_t ct, AlleleCode** allele_arr_ptr) {
291 *allele_arr_ptr = S_CAST(AlleleCode*, bigstack_alloc(ct * sizeof(AlleleCode)));
292 return !(*allele_arr_ptr);
293 }
294
bigstack_alloc_kac(uintptr_t ct,const AlleleCode ** allele_arr_ptr)295 HEADER_INLINE BoolErr bigstack_alloc_kac(uintptr_t ct, const AlleleCode** allele_arr_ptr) {
296 *allele_arr_ptr = S_CAST(const AlleleCode*, bigstack_alloc(ct * sizeof(AlleleCode)));
297 return !(*allele_arr_ptr);
298 }
299
bigstack_alloc_dosage(uintptr_t ct,Dosage ** dosage_arr_ptr)300 HEADER_INLINE BoolErr bigstack_alloc_dosage(uintptr_t ct, Dosage** dosage_arr_ptr) {
301 *dosage_arr_ptr = S_CAST(Dosage*, bigstack_alloc(ct * sizeof(Dosage)));
302 return !(*dosage_arr_ptr);
303 }
304
bigstack_alloc_dphase(uintptr_t ct,SDosage ** dphase_arr_ptr)305 HEADER_INLINE BoolErr bigstack_alloc_dphase(uintptr_t ct, SDosage** dphase_arr_ptr) {
306 *dphase_arr_ptr = S_CAST(SDosage*, bigstack_alloc(ct * sizeof(SDosage)));
307 return !(*dphase_arr_ptr);
308 }
309
310
bigstack_alloc_acp(uintptr_t ct,AlleleCode *** acp_arr_ptr)311 HEADER_INLINE BoolErr bigstack_alloc_acp(uintptr_t ct, AlleleCode*** acp_arr_ptr) {
312 *acp_arr_ptr = S_CAST(AlleleCode**, bigstack_alloc(ct * sizeof(intptr_t)));
313 return !(*acp_arr_ptr);
314 }
315
bigstack_alloc_dosagep(uintptr_t ct,Dosage *** dosagep_arr_ptr)316 HEADER_INLINE BoolErr bigstack_alloc_dosagep(uintptr_t ct, Dosage*** dosagep_arr_ptr) {
317 *dosagep_arr_ptr = S_CAST(Dosage**, bigstack_alloc(ct * sizeof(intptr_t)));
318 return !(*dosagep_arr_ptr);
319 }
320
bigstack_alloc_dphasep(uintptr_t ct,SDosage *** dphasep_arr_ptr)321 HEADER_INLINE BoolErr bigstack_alloc_dphasep(uintptr_t ct, SDosage*** dphasep_arr_ptr) {
322 *dphasep_arr_ptr = S_CAST(SDosage**, bigstack_alloc(ct * sizeof(intptr_t)));
323 return !(*dphasep_arr_ptr);
324 }
325
326
bigstack_end_alloc_dosage(uintptr_t ct,Dosage ** dosage_arr_ptr)327 HEADER_INLINE BoolErr bigstack_end_alloc_dosage(uintptr_t ct, Dosage** dosage_arr_ptr) {
328 *dosage_arr_ptr = S_CAST(Dosage*, bigstack_end_alloc(ct * sizeof(Dosage)));
329 return !(*dosage_arr_ptr);
330 }
331
bigstack_end_alloc_dphase(uintptr_t ct,SDosage ** dphase_arr_ptr)332 HEADER_INLINE BoolErr bigstack_end_alloc_dphase(uintptr_t ct, SDosage** dphase_arr_ptr) {
333 *dphase_arr_ptr = S_CAST(SDosage*, bigstack_end_alloc(ct * sizeof(SDosage)));
334 return !(*dphase_arr_ptr);
335 }
336
337 BoolErr BigstackAllocPgv(uint32_t sample_ct, uint32_t multiallelic_needed, PgenGlobalFlags gflags, PgenVariant* pgvp);
338
339 // remainder must be in [1, 16383].
340 char* PrintDosageDecimal(uint32_t remainder, char* start);
341
342 // small_dosage must be in [0, kDosageMid * 10 - 1].
PrintSmallDosage(uint32_t small_dosage,char * start)343 HEADER_INLINE char* PrintSmallDosage(uint32_t small_dosage, char* start) {
344 *start++ = '0' + (small_dosage / kDosageMid);
345 const uint32_t remainder = small_dosage % kDosageMid;
346 if (!remainder) {
347 return start;
348 }
349 return PrintDosageDecimal(remainder, start);
350 }
351
352 // 3 decimal places. dosage on /kDosageMax rather than /kDosageMid scale
353 // (hence the extra 'd')
354 char* ddosagetoa(uint64_t dosage, char* start);
355
356 // remainder must be in [1, 32767].
357 char* PrintDdosageDecimal(uint32_t remainder, char* start);
358
359 // 5 decimal places. Only used when it is important to be able to reconstruct
360 // the exact original value.
ddosagetoa_full(uint64_t dosage,char * start)361 HEADER_INLINE char* ddosagetoa_full(uint64_t dosage, char* start) {
362 start = u32toa(dosage / kDosageMax, start);
363 const uint32_t remainder = dosage % kDosageMax;
364 if (!remainder) {
365 return start;
366 }
367 return PrintDdosageDecimal(remainder, start);
368 }
369
ZeroDosageArr(uintptr_t entry_ct,Dosage * dosage_arr)370 HEADER_INLINE void ZeroDosageArr(uintptr_t entry_ct, Dosage* dosage_arr) {
371 memset(dosage_arr, 0, entry_ct * sizeof(Dosage));
372 }
373
ZeroDphaseArr(uintptr_t entry_ct,SDosage * dphase_arr)374 HEADER_INLINE void ZeroDphaseArr(uintptr_t entry_ct, SDosage* dphase_arr) {
375 memset(dphase_arr, 0, entry_ct * sizeof(SDosage));
376 }
377
SetAllDosageArr(uintptr_t entry_ct,Dosage * dosage_arr)378 HEADER_INLINE void SetAllDosageArr(uintptr_t entry_ct, Dosage* dosage_arr) {
379 memset(dosage_arr, 255, entry_ct * sizeof(Dosage));
380 }
381
382 void PopulateDenseDosage(const uintptr_t* genoarr, const uintptr_t* dosage_present, const Dosage* dosage_main, uint32_t sample_ct, uint32_t dosage_ct, Dosage* dense_dosage);
383
384 void PopulateRescaledDosage(const uintptr_t* genoarr, const uintptr_t* dosage_present, const Dosage* dosage_main, double slope, double intercept, double missing_val, uint32_t sample_ct, uint32_t dosage_ct, double* expanded_dosages);
385
386 // assumes trailing bits of genoarr are zeroed out
AtLeastOneHetUnsafe(const uintptr_t * genoarr,uint32_t sample_ct)387 HEADER_INLINE uint32_t AtLeastOneHetUnsafe(const uintptr_t* genoarr, uint32_t sample_ct) {
388 const uint32_t sample_ctl2 = DivUp(sample_ct, kBitsPerWordD2);
389 for (uint32_t uii = 0; uii != sample_ctl2; ++uii) {
390 const uintptr_t geno_word = genoarr[uii];
391 if (Word01(geno_word)) {
392 return 1;
393 }
394 }
395 return 0;
396 }
397
398 uint32_t AtLeastOneMultiallelicHet(const PgenVariant* pgvp, uint32_t sample_ct);
399
400 void SetHetMissing(uintptr_t word_ct, uintptr_t* genovec);
401
402 void SetHetMissingCleardosage(uintptr_t word_ct, uintptr_t* genovec, uint32_t* write_dosage_ct_ptr, uintptr_t* dosagepresent, Dosage* dosage_main);
403
404 void SetHetMissingKeepdosage(uintptr_t word_ct, uintptr_t* genovec, uint32_t* write_dosage_ct_ptr, uintptr_t* dosagepresent, Dosage* dosage_main);
405
406 void GenoarrToNonmissing(const uintptr_t* genoarr, uint32_t sample_ctl2, uintptr_t* nonmissing_bitarr);
407
408 uint32_t GenoarrCountMissingUnsafe(const uintptr_t* genoarr, uint32_t sample_ct);
409
410 uint32_t GenoarrCountMissingSubset(const uintptr_t* genoarr, const uintptr_t* interleaved_vec, uint32_t sample_ct);
411
412 uint32_t GenoarrCountMissingInvsubsetUnsafe(const uintptr_t* genoarr, const uintptr_t* exclude_mask, uint32_t sample_ct);
413
414 // See also DataFidColIsRequired(), etc. in plink2_data.h, which checks whether
415 // at least one remaining value is nonzero.
FidColIsRequired(const SampleIdInfo * siip,uint32_t maybe_modifier)416 HEADER_INLINE uint32_t FidColIsRequired(const SampleIdInfo* siip, uint32_t maybe_modifier) {
417 return (maybe_modifier & 2) || ((maybe_modifier & 1) && (siip->flags & kfSampleIdFidPresent));
418 }
419
SidColIsRequired(const char * sids,uint32_t maybe_modifier)420 HEADER_INLINE uint32_t SidColIsRequired(const char* sids, uint32_t maybe_modifier) {
421 return (maybe_modifier & 2) || ((maybe_modifier & 1) && sids);
422 }
423
ParentalColsAreRequired(const PedigreeIdInfo * piip,uint32_t maybe_modifier)424 HEADER_INLINE uint32_t ParentalColsAreRequired(const PedigreeIdInfo* piip, uint32_t maybe_modifier) {
425 return (maybe_modifier & 2) || ((maybe_modifier & 1) && (piip->sii.flags & kfSampleIdParentsPresent));
426 }
427
428 void CollapsedSampleFmtidInit(const uintptr_t* sample_include, const SampleIdInfo* siip, uint32_t sample_ct, uint32_t include_fid, uint32_t include_sid, uintptr_t max_sample_fmtid_blen, char* collapsed_sample_fmtids_iter);
429
GetMaxSampleFmtidBlen(const SampleIdInfo * siip,uint32_t include_fid,uint32_t include_sid)430 HEADER_INLINE uintptr_t GetMaxSampleFmtidBlen(const SampleIdInfo* siip, uint32_t include_fid, uint32_t include_sid) {
431 uintptr_t max_sid_blen = 0;
432 if (include_sid) {
433 max_sid_blen = siip->sids? siip->max_sid_blen : 2;
434 }
435 return siip->max_sample_id_blen + max_sid_blen - 2 * (!include_fid);
436 }
437
438 BoolErr CollapsedSampleFmtidInitAlloc(const uintptr_t* sample_include, const SampleIdInfo* siip, uint32_t sample_ct, uint32_t include_fid, uint32_t include_sid, char** collapsed_sample_fmtids_ptr, uintptr_t* max_sample_fmtid_blen_ptr);
439
440 // Assumes sample_ct > 0.
441 uint32_t OnlyOneFid(const uintptr_t* sample_include, const SampleIdInfo* siip, uint32_t sample_ct);
442
443 uint32_t GetMajIdxMulti(const double* cur_allele_freqs, uint32_t cur_allele_ct);
444
GetMajIdx(const double * cur_allele_freqs,uint32_t cur_allele_ct)445 HEADER_INLINE uint32_t GetMajIdx(const double* cur_allele_freqs, uint32_t cur_allele_ct) {
446 if (cur_allele_freqs[0] >= 0.5) {
447 return 0;
448 }
449 if (cur_allele_ct == 2) {
450 return 1;
451 }
452 return GetMajIdxMulti(cur_allele_freqs, cur_allele_ct);
453 }
454
GetNonmajFreq(const double * cur_allele_freqs,uint32_t allele_ct)455 HEADER_INLINE double GetNonmajFreq(const double* cur_allele_freqs, uint32_t allele_ct) {
456 double tot_nonlast_freq = cur_allele_freqs[0];
457 double max_freq = tot_nonlast_freq;
458 const uint32_t allele_ct_m1 = allele_ct - 1;
459 for (uint32_t allele_idx = 1; allele_idx != allele_ct_m1; ++allele_idx) {
460 const double cur_alt_freq = cur_allele_freqs[allele_idx];
461 tot_nonlast_freq += cur_alt_freq;
462 if (cur_alt_freq > max_freq) {
463 max_freq = cur_alt_freq;
464 }
465 }
466 const double nonmajor_freq = 1.0 - max_freq;
467 return MINV(nonmajor_freq, tot_nonlast_freq);
468 }
469
GetAlleleFreq(const double * cur_allele_freqs,uint32_t allele_idx,uint32_t allele_ct)470 HEADER_INLINE double GetAlleleFreq(const double* cur_allele_freqs, uint32_t allele_idx, uint32_t allele_ct) {
471 const uint32_t allele_ct_m1 = allele_ct - 1;
472 if (allele_idx < allele_ct_m1) {
473 return cur_allele_freqs[allele_idx];
474 }
475 double last_freq = 1.0 - cur_allele_freqs[0];
476 for (uint32_t tmp_allele_idx = 1; tmp_allele_idx != allele_ct_m1; ++tmp_allele_idx) {
477 last_freq -= cur_allele_freqs[tmp_allele_idx];
478 }
479 return MAXV(last_freq, 0.0);
480 }
481
482
483 // Functions with Xid in their name deal with both FID/IID (with a single tab
484 // separator) and FID/IID/SID (two tabs) sample IDs. (Missing FID is
485 // represented as FID='0'.)
486
487 // With no header line, --keep/--remove and the like should interpret a single
488 // token as IID (treating FID as 0), and two tokens as FID/IID. --no-id-header
489 // does not support IID/SID output, so we don't need to worry about supporting
490 // FidIidSidOrIidSid.
491 // With a header line, all four {FID present/absent, SID present/absent}
492 // combinations are allowed.
493 FLAGSET_DEF_START()
494 kfXidMode0,
495
496 kfXidModeFlagOneTokenOk = (1 << 0),
497 kfXidModeFlagNeverFid = (1 << 1),
498 kfXidModeFlagSid = (1 << 2),
499
500 kfXidModeFidIid = 0,
501 kfXidModeFidIidOrIid = kfXidModeFlagOneTokenOk,
502 kfXidModeIid = (kfXidModeFlagOneTokenOk | kfXidModeFlagNeverFid),
503 kfXidModeFidIidSid = kfXidModeFlagSid,
504 kfXidModeIidSid = (kfXidModeFlagNeverFid | kfXidModeFlagSid)
505 FLAGSET_DEF_END(XidMode);
506
507 // Assumes fixed-width.
GetXidColCt(XidMode xid_mode)508 HEADER_INLINE uint32_t GetXidColCt(XidMode xid_mode) {
509 return 2 + (xid_mode == kfXidModeFidIidSid) - (xid_mode == kfXidModeIid);
510 }
511
512 // sample_xid_map allocated on bottom, to play well with --indiv-sort
513 PglErr SortedXidboxInitAlloc(const uintptr_t* sample_include, const SampleIdInfo* siip, uint32_t sample_ct, uint32_t allow_dups, XidMode xid_mode, uint32_t use_nsort, char** sorted_xidbox_ptr, uint32_t** xid_map_ptr, uintptr_t* max_xid_blen_ptr);
514
515 // returns slen for ID, or 0 on guaranteed mismatch (longer than max_xid_blen)
516 // or parse failure (*readpp set to nullptr in latter case).
517 uint32_t XidRead(uintptr_t max_xid_blen, uint32_t comma_delim, XidMode xid_mode, const char** read_pp, char* __restrict idbuf);
518
519 // returns 1 on missing token *or* if the sample ID is not present. cases can
520 // be distinguished by checking whether *read_pp_new == nullptr: if it is, a
521 // missing-tokens error should probably be reported.
522 // sample_id_map == nullptr is permitted
523 // *read_pp is now set to point to the end of IID/SID instead of the beginning
524 // of the next token; this is a change from plink 1.9.
SortedXidboxReadFind(const char * __restrict sorted_xidbox,const uint32_t * __restrict xid_map,uintptr_t max_xid_blen,uintptr_t xid_ct,uint32_t comma_delim,XidMode xid_mode,const char ** read_pp,uint32_t * sample_uidx_ptr,char * __restrict idbuf)525 HEADER_INLINE BoolErr SortedXidboxReadFind(const char* __restrict sorted_xidbox, const uint32_t* __restrict xid_map, uintptr_t max_xid_blen, uintptr_t xid_ct, uint32_t comma_delim, XidMode xid_mode, const char** read_pp, uint32_t* sample_uidx_ptr, char* __restrict idbuf) {
526 const uint32_t slen_final = XidRead(max_xid_blen, comma_delim, xid_mode, read_pp, idbuf);
527 if (!slen_final) {
528 return 1;
529 }
530 return SortedIdboxFind(idbuf, sorted_xidbox, xid_map, slen_final, max_xid_blen, xid_ct, sample_uidx_ptr);
531 }
532
533 // Matches a sample ID *prefix*. Thus, if FID/IID/SID is loaded, but the input
534 // file contains just FID/IID, and there are some FID/IID pairs which
535 // correspond to multiple samples, this lets you iterate over all of them.
536 // (Caller is responsible for looking up xid_map[] to perform xid_idx ->
537 // sample_uidx conversions.)
538 BoolErr SortedXidboxReadMultifind(const char* __restrict sorted_xidbox, uintptr_t max_xid_blen, uintptr_t xid_ct, uint32_t comma_delim, XidMode xid_mode, const char** read_pp, uint32_t* __restrict xid_idx_start_ptr, uint32_t* __restrict xid_idx_end_ptr, char* __restrict idbuf);
539
540 FLAGSET_DEF_START()
541 kfXidHeader0,
542
543 kfXidHeaderFixedWidth = (1 << 0),
544 kfXidHeaderIgnoreSid = (1 << 1),
545 kfXidHeaderFixedWidthIgnoreSid = (kfXidHeaderFixedWidth | kfXidHeaderIgnoreSid)
546 FLAGSET_DEF_END(XidHeaderFlags);
547
548 // - May return kPglRetEof, or other TextStream errors.
549 // - line_startp can be nullptr. If it isn't, it's set to the (lstripped)
550 // beginning of the current line.
551 // - line_iterp can be nullptr. If it isn't, it's set to the beginning of the
552 // first unparsed (i.e. not FID/IID/SID) token in the header line, or
553 // line_start if there's no header line.
554 // - line_idx must be zero unless initial lines were skipped.
555 // - If no header line is present, xid_mode will be set to kfXidModeFidIid if
556 // kfXidHeaderFixedWidth is set, and kfXidModeFidIidOrIid (which tolerates a
557 // mix of single-token and multitoken lines) otherwise.
558 // - TSTREAM_FAIL errors are not printed, but other errors are.
559 PglErr LoadXidHeader(const char* flag_name, XidHeaderFlags xid_header_flags, uintptr_t* line_idx_ptr, TextStream* txsp, XidMode* xid_mode_ptr, char** line_startp, char** line_iterp);
560
561 PglErr OpenAndLoadXidHeader(const char* fname, const char* flag_name, XidHeaderFlags xid_header_flags, uint32_t max_line_blen, TextStream* txsp, XidMode* xid_mode_ptr, uintptr_t* line_idx_ptr, char** line_startp, char** line_iterp);
562
563 // header line expected to start with FID1, ID1, or IID1
564 PglErr LoadXidHeaderPair(const char* flag_name, uint32_t sid_over_fid, uintptr_t* line_idx_ptr, TextStream* txsp, XidMode* xid_mode_ptr, char** line_startp, char** line_iterp);
565
566
567 // note that this is no longer divisible by 64
568 CONSTI32(kMaxContigs, 65274);
569 CONSTI32(kMaxChrCodeDigits, 5);
570
571 // change ChrIdx to uint32_t if (kMaxContigs + kChrOffsetCt) > 65536
572 typedef uint16_t ChrIdx;
573
574 // get_htable_min_size(kChrRawEnd) (use constexpr once sufficient
575 // compiler support is available)
576 // (not get_htable_fast_size since, an overwhelming majority of the time, we'll
577 // have far fewer than 2^16 codes)
578 CONSTI32(kChrHtableSize, 130579);
579
580 // (note that n+1, n+2, n+3, and n+4 are reserved for X/Y/XY/MT)
581 CONSTI32(kMaxChrTextnum, 95);
582
583 // get_chr_code_raw() needs to be modified if this changes
584 CONSTI32(kMaxChrTextnumSlen, 2);
585
586 ENUM_U31_DEF_START()
587 kChrOffsetX,
588 kChrOffsetY,
589
590 // old way of representing pseudo-autosomal regions. clumsy since this
591 // required changing chromosome order
592 kChrOffsetXY,
593
594 kChrOffsetMT,
595
596 // plink 2.x pseudo-autosomal regions.
597 kChrOffsetPAR1,
598 kChrOffsetPAR2,
599 kChrOffsetCt
600 ENUM_U31_DEF_END(XymtOffset);
601
602 CONSTI32(kChrRawX, kMaxContigs + kChrOffsetX);
603 CONSTI32(kChrRawY, kMaxContigs + kChrOffsetY);
604 CONSTI32(kChrRawXY, kMaxContigs + kChrOffsetXY);
605 CONSTI32(kChrRawMT, kMaxContigs + kChrOffsetMT);
606 CONSTI32(kChrRawPAR1, kMaxContigs + kChrOffsetPAR1);
607 CONSTI32(kChrRawPAR2, kMaxContigs + kChrOffsetPAR2);
608 CONSTI32(kChrRawEnd, kMaxContigs + kChrOffsetCt);
609
610 static_assert((!(kChrRawEnd % kBitsPerWord)), "kChrRawEnd expression must be updated.");
611 CONSTI32(kChrMaskWords, kChrRawEnd / kBitsPerWord);
612
613 #ifdef __LP64__
614 CONSTI32(kChrExcludeWords, 2);
615 #else
616 CONSTI32(kChrExcludeWords, 4);
617 #endif
618 static_assert(kChrExcludeWords * kBitsPerWord >= kMaxChrTextnum + 2 * kChrOffsetCt + 1, "kChrExcludeWords must be updated.");
619
620 ENUM_U31_DEF_START()
621 kChrsetSourceDefault,
622 kChrsetSourceCmdline,
623 kChrsetSourceFile
624 ENUM_U31_DEF_END(ChrsetSource);
625
626 FLAGSET_DEF_START()
627 kfChrOutput0,
628 kfChrOutputPrefix = (1 << 0),
629 kfChrOutputM = (1 << 1),
630 kfChrOutputMT = (1 << 2),
631 kfChrOutput0M = (1 << 3)
632 FLAGSET_DEF_END(ChrOutput);
633
634 typedef struct ChrInfoStruct {
635 // Main dynamic block intended to be allocated as a single aligned block of
636 // memory on the heap freeable with vecaligned_free(), with chr_mask at the
637 // base.
638 NONCOPYABLE(ChrInfoStruct);
639
640 uintptr_t* chr_mask; // which chromosomes aren't known to be absent?
641
642 // This includes chrX. As of alpha 2, it also includes MT again (like plink
643 // 1.07, and unlike 1.9 and 2.0a1), now that enough dosage functionality is
644 // in place.
645 uintptr_t* haploid_mask;
646
647 // order of chromosomes in input files
648 // currently tolerates out-of-order chromosomes, as long as all variants for
649 // any given chromosome are together
650 uint32_t* chr_file_order;
651
652 // if the second chromosome in the dataset is chr5, chr_file_order[1] == 5,
653 // the raw variant indexes for chr5 are in [chr_fo_vidx_start[1],
654 // chr_fo_vidx_start[2]). and chr_idx_to_foidx[5] == 1.
655 uint32_t* chr_fo_vidx_start;
656 uint32_t* chr_idx_to_foidx;
657
658 // --allow-extra-chr support
659 const char** nonstd_names;
660 uint32_t* nonstd_id_htable;
661 // end main dynamic block
662
663 uint32_t chr_ct; // number of distinct chromosomes/contigs
664 ChrsetSource chrset_source;
665
666 uintptr_t chr_exclude[kChrExcludeWords];
667
668 // X, Y, XY...; UINT32_MAXM1 = not in chromosome set
669 STD_ARRAY_DECL(uint32_t, kChrOffsetCt, xymt_codes);
670
671 uint32_t max_numeric_code;
672 uint32_t max_code; // no longer identical to max_numeric_code, with PARs
673
674 uint32_t autosome_ct;
675
676 // yet more --allow-extra-chr support
677 uint32_t zero_extra_chrs;
678 uint32_t name_ct;
679 LlStr* incl_excl_name_stack;
680 uint32_t is_include_stack;
681 ChrOutput output_encoding;
682 } ChrInfo;
683
684 extern const char g_xymt_log_names[kChrOffsetCt][5];
685
686 PglErr InitChrInfo(ChrInfo* cip);
687
688 void FinalizeChrset(MiscFlags misc_flags, ChrInfo* cip);
689
InitChrInfoHuman(ChrInfo * cip)690 HEADER_INLINE PglErr InitChrInfoHuman(ChrInfo* cip) {
691 // convenience wrapper
692 if (unlikely(InitChrInfo(cip))) {
693 return kPglRetNomem;
694 }
695 FinalizeChrset(kfMisc0, cip);
696 return kPglRetSuccess;
697 }
698
699 void ForgetExtraChrNames(uint32_t reinitialize, ChrInfo* cip);
700
701 // in the usual case where the number of chromosomes/contigs is much less than
702 // kMaxContigs, this reduces chr_info's memory consumption and improves
703 // locality.
704 PglErr FinalizeChrInfo(ChrInfo* cip);
705
706 void CleanupChrInfo(ChrInfo* cip);
707
708 // assumes chr_idx is valid
709 // note that chr_idx == 0 is always rendered as '0', never 'chr0'
710 char* chrtoa(const ChrInfo* cip, uint32_t chr_idx, char* buf);
711
712 uint32_t GetMaxChrSlen(const ChrInfo* cip);
713
714 uint32_t IsHaploidChrPresent(const ChrInfo* cip);
715
716 uint32_t IsAutosomalDiploidChrPresent(const ChrInfo* cip);
717
718 // any character <= ' ' is considered a terminator
719 // maps chrX -> kChrRawX, etc.
720 uint32_t GetChrCodeRaw(const char* str_iter);
721
722 // requires chr_name to be null-terminated
723 // maps chrX -> xymt_codes[kChrOffsetX], etc.
724 // error codes:
725 // UINT32_MAX = --allow-extra-chr ok
726 // UINT32_MAXM1 = total fail
727 uint32_t GetChrCode(const char* chr_name, const ChrInfo* cip, uint32_t name_slen);
728
729 // when the chromosome name isn't null-terminated
730 // requires chr_name[name_slen] to be mutable
731 uint32_t GetChrCodeCounted(const ChrInfo* cip, uint32_t name_slen, char* chr_name);
732
GetVariantChrFoIdx(const ChrInfo * cip,uintptr_t variant_uidx)733 HEADER_INLINE uint32_t GetVariantChrFoIdx(const ChrInfo* cip, uintptr_t variant_uidx) {
734 return CountSortedSmallerU32(&(cip->chr_fo_vidx_start[1]), cip->chr_ct, variant_uidx + 1);
735 }
736
GetVariantChr(const ChrInfo * cip,uintptr_t variant_uidx)737 HEADER_INLINE uint32_t GetVariantChr(const ChrInfo* cip, uintptr_t variant_uidx) {
738 return cip->chr_file_order[GetVariantChrFoIdx(cip, variant_uidx)];
739 }
740
XymtExists(const ChrInfo * cip,uint32_t xymt_offset,uint32_t * xymt_code_ptr)741 HEADER_INLINE uint32_t XymtExists(const ChrInfo* cip, uint32_t xymt_offset, uint32_t* xymt_code_ptr) {
742 // too easy to forget IsSet(chr_mask) check if we don't use this
743 const uint32_t xymt_code = cip->xymt_codes[xymt_offset];
744 *xymt_code_ptr = xymt_code;
745 return (!IsI32Neg(xymt_code)) && IsSet(cip->chr_mask, xymt_code);
746 }
747
GetXymtStartAndEnd(const ChrInfo * cip,uint32_t xymt_offset,uint32_t * xymt_start_ptr,uint32_t * xymt_end_ptr)748 HEADER_INLINE void GetXymtStartAndEnd(const ChrInfo* cip, uint32_t xymt_offset, uint32_t* xymt_start_ptr, uint32_t* xymt_end_ptr) {
749 uint32_t xymt_code;
750 if (!XymtExists(cip, xymt_offset, &xymt_code)) {
751 *xymt_start_ptr = 0;
752 *xymt_end_ptr = 0;
753 return;
754 }
755 const uint32_t chr_fo_idx = cip->chr_idx_to_foidx[xymt_code];
756 *xymt_start_ptr = cip->chr_fo_vidx_start[chr_fo_idx];
757 *xymt_end_ptr = cip->chr_fo_vidx_start[chr_fo_idx + 1];
758 }
759
GetXymtCodeStartAndEndUnsafe(const ChrInfo * cip,uint32_t xymt_offset,uint32_t * xymt_code_ptr,uint32_t * xymt_start_ptr,uint32_t * xymt_end_ptr)760 HEADER_INLINE void GetXymtCodeStartAndEndUnsafe(const ChrInfo* cip, uint32_t xymt_offset, uint32_t* xymt_code_ptr, uint32_t* xymt_start_ptr, uint32_t* xymt_end_ptr) {
761 // assumes xymt_exists was previously called, and is true
762 const uint32_t xymt_code = cip->xymt_codes[xymt_offset];
763 *xymt_code_ptr = xymt_code;
764 const uint32_t chr_fo_idx = cip->chr_idx_to_foidx[xymt_code];
765 *xymt_start_ptr = cip->chr_fo_vidx_start[chr_fo_idx];
766 *xymt_end_ptr = cip->chr_fo_vidx_start[chr_fo_idx + 1];
767 }
768
769 // now assumes chr_name is null-terminated
770 PglErr TryToAddChrName(const char* chr_name, const char* file_descrip, uintptr_t line_idx, uint32_t name_slen, uint32_t allow_extra_chrs, uint32_t* chr_idx_ptr, ChrInfo* cip);
771
GetOrAddChrCode(const char * chr_name,const char * file_descrip,uintptr_t line_idx,uint32_t name_slen,uint32_t allow_extra_chrs,ChrInfo * cip,uint32_t * chr_idx_ptr)772 HEADER_INLINE PglErr GetOrAddChrCode(const char* chr_name, const char* file_descrip, uintptr_t line_idx, uint32_t name_slen, uint32_t allow_extra_chrs, ChrInfo* cip, uint32_t* chr_idx_ptr) {
773 *chr_idx_ptr = GetChrCode(chr_name, cip, name_slen);
774 if (!IsI32Neg(*chr_idx_ptr)) {
775 return kPglRetSuccess;
776 }
777 return TryToAddChrName(chr_name, file_descrip, line_idx, name_slen, allow_extra_chrs, chr_idx_ptr, cip);
778 }
779
GetOrAddChrCodeDestructive(const char * file_descrip,uintptr_t line_idx,uint32_t allow_extra_chrs,char * chr_name,char * chr_name_end,ChrInfo * cip,uint32_t * chr_idx_ptr)780 HEADER_INLINE PglErr GetOrAddChrCodeDestructive(const char* file_descrip, uintptr_t line_idx, uint32_t allow_extra_chrs, char* chr_name, char* chr_name_end, ChrInfo* cip, uint32_t* chr_idx_ptr) {
781 *chr_name_end = '\0';
782 return GetOrAddChrCode(chr_name, file_descrip, line_idx, chr_name_end - chr_name, allow_extra_chrs, cip, chr_idx_ptr);
783 }
784
785 // Assumes sample_ct positive. Does not require trailing bits to be clear.
786 uint32_t AllGenoEqual(const uintptr_t* genoarr, uint32_t sample_ct);
787
788 // zeroes out samples not in the mask
789 void InterleavedMaskZero(const uintptr_t* __restrict interleaved_mask, uintptr_t vec_ct, uintptr_t* __restrict genovec);
790
791 // sets samples outside the mask to missing (0b11)
792 void InterleavedMaskMissing(const uintptr_t* __restrict interleaved_set, uintptr_t vec_ct, uintptr_t* __restrict genovec);
793
794 // sets samples in the mask to missing (0b11)
795 void InterleavedSetMissing(const uintptr_t* __restrict interleaved_set, uintptr_t vec_ct, uintptr_t* __restrict genovec);
796
797 void InterleavedSetMissingCleardosage(const uintptr_t* __restrict orig_set, const uintptr_t* __restrict interleaved_set, uintptr_t vec_ct, uintptr_t* __restrict genovec, uint32_t* __restrict write_dosage_ct_ptr, uintptr_t* __restrict dosagepresent, Dosage* dosage_main);
798
799 void SetMaleHetMissing(const uintptr_t* __restrict sex_male_interleaved, uint32_t vec_ct, uintptr_t* __restrict genovec);
800
801 void EraseMaleDphases(const uintptr_t* __restrict sex_male, uint32_t* __restrict write_dphase_ct_ptr, uintptr_t* __restrict dphasepresent, SDosage* dphase_delta);
802
803 void SetMaleHetMissingCleardosage(const uintptr_t* __restrict sex_male, const uintptr_t* __restrict sex_male_interleaved, uint32_t vec_ct, uintptr_t* __restrict genovec, uint32_t* __restrict write_dosage_ct_ptr, uintptr_t* __restrict dosagepresent, Dosage* dosage_main);
804
805 void SetMaleHetMissingKeepdosage(const uintptr_t* __restrict sex_male, const uintptr_t* __restrict sex_male_interleaved, uint32_t word_ct, uintptr_t* __restrict genovec, uint32_t* __restrict write_dosage_ct_ptr, uintptr_t* __restrict dosagepresent, Dosage* dosage_main);
806
807 // Clears each bit in bitarr which doesn't correspond to a genovec het.
808 // Assumes that either trailing bits of bitarr are already zero, or trailing
809 // bits of genovec are zero.
810 void MaskGenoarrHetsUnsafe(const uintptr_t* __restrict genoarr, uint32_t raw_sample_ctl2, uintptr_t* __restrict bitarr);
811
812 void MaskGenoarrHetsMultiallelicUnsafe(const uintptr_t* __restrict genoarr, const uintptr_t* __restrict patch_10_set, const AlleleCode* __restrict patch_10_vals, uint32_t raw_sample_ctl2, uintptr_t* __restrict bitarr);
813
814 // vertical popcount support
815 // VcountScramble1() and friends in plink2_cmdline
816 #ifdef __LP64__
817 # ifdef USE_AVX2
818 // 2->4: 0 2 ... 126 1 3 ... 127
819 // 4->8: 0 4 ... 124 2 6 ... 126 1 5 ... 125 3 7 ... 127
820 // 8->32: 0 16 ... 112 4 20 ... 116 ... 124 2 18 ... 114 6 22 ... 118 ... 126 1 17 ...
VcountScramble2(uint32_t orig_idx)821 HEADER_INLINE uint32_t VcountScramble2(uint32_t orig_idx) {
822 return (orig_idx & (~127)) + ((orig_idx & 1) * 64) + ((orig_idx & 2) * 16) + ((orig_idx & 12) * 2) + ((orig_idx & 112) / 16);
823 }
824 # else
VcountScramble2(uint32_t orig_idx)825 HEADER_INLINE uint32_t VcountScramble2(uint32_t orig_idx) {
826 return (orig_idx & (~63)) + ((orig_idx & 1) * 32) + ((orig_idx & 2) * 8) + (orig_idx & 12) + ((orig_idx & 48) / 16);
827 }
828 # endif
829 #else
830 // 2->4: 0 2 4 6 8 10 12 14 1 3 5 ...
831 // 4->8: 0 4 8 12 2 6 10 14 1 5 9 ...
832 // 8->32: 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15
VcountScramble2(uint32_t orig_idx)833 HEADER_INLINE uint32_t VcountScramble2(uint32_t orig_idx) {
834 return (orig_idx & (~15)) + ((orig_idx & 1) * 8) + ((orig_idx & 2) * 2) + ((orig_idx & 12) / 4);
835 }
836 #endif
837
VcountIncr2To4(const VecW * acc2_iter,uint32_t acc2_vec_ct,VecW * acc4_iter)838 HEADER_INLINE void VcountIncr2To4(const VecW* acc2_iter, uint32_t acc2_vec_ct, VecW* acc4_iter) {
839 const VecW m2 = VCONST_W(kMask3333);
840 for (uint32_t vidx = 0; vidx != acc2_vec_ct; ++vidx) {
841 VecW loader = *acc2_iter++;
842 *acc4_iter += loader & m2;
843 ++acc4_iter;
844 *acc4_iter += vecw_srli(loader, 2) & m2;
845 ++acc4_iter;
846 }
847 }
848
Vcount0Incr2To4(uint32_t acc2_vec_ct,VecW * acc2_iter,VecW * acc4_iter)849 HEADER_INLINE void Vcount0Incr2To4(uint32_t acc2_vec_ct, VecW* acc2_iter, VecW* acc4_iter) {
850 const VecW m2 = VCONST_W(kMask3333);
851 for (uint32_t vidx = 0; vidx != acc2_vec_ct; ++vidx) {
852 VecW loader = *acc2_iter;
853 *acc2_iter++ = vecw_setzero();
854 *acc4_iter += loader & m2;
855 ++acc4_iter;
856 *acc4_iter += vecw_srli(loader, 2) & m2;
857 ++acc4_iter;
858 }
859 }
860
861
862 // uint32_t chr_window_max(const uintptr_t* variant_include, const ChrInfo* cip, const uint32_t* variant_pos, uint32_t chr_fo_idx, uint32_t ct_max, uint32_t bp_max, uint32_t cur_window_max);
863
CountChrVariantsUnsafe(const uintptr_t * variant_include,const ChrInfo * cip,uint32_t chr_idx)864 HEADER_INLINE uint32_t CountChrVariantsUnsafe(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t chr_idx) {
865 assert(IsSet(cip->chr_mask, chr_idx));
866 const uint32_t chr_fo_idx = cip->chr_idx_to_foidx[chr_idx];
867 const uint32_t min_idx = cip->chr_fo_vidx_start[chr_fo_idx];
868 const uint32_t max_idx = cip->chr_fo_vidx_start[chr_fo_idx + 1];
869 return PopcountBitRange(variant_include, min_idx, max_idx);
870 }
871
ChrIsNonempty(const uintptr_t * variant_include,const ChrInfo * cip,uint32_t chr_idx)872 HEADER_INLINE uint32_t ChrIsNonempty(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t chr_idx) {
873 if (!IsSet(cip->chr_mask, chr_idx)) {
874 return 0;
875 }
876 const uint32_t chr_fo_idx = cip->chr_idx_to_foidx[chr_idx];
877 const uint32_t min_idx = cip->chr_fo_vidx_start[chr_fo_idx];
878 const uint32_t max_idx = cip->chr_fo_vidx_start[chr_fo_idx + 1];
879 return !AllBitsAreZero(variant_include, min_idx, max_idx);
880 }
881
XymtIsNonempty(const uintptr_t * variant_include,const ChrInfo * cip,uint32_t xymt_offset)882 HEADER_INLINE uint32_t XymtIsNonempty(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t xymt_offset) {
883 const uint32_t xymt_code = cip->xymt_codes[xymt_offset];
884 if (IsI32Neg(xymt_code) || (!IsSet(cip->chr_mask, xymt_code))) {
885 return 0;
886 }
887 return ChrIsNonempty(variant_include, cip, xymt_code);
888 }
889
890 // assumes there's at least one variant on specified chromosome
891 uint32_t NotOnlyXymt(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t raw_variant_ct, uint32_t xymt_offset);
892
893 uint32_t CountNonAutosomalVariants(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t count_x, uint32_t count_mt);
894
895 void ExcludeNonAutosomalVariants(const ChrInfo* cip, uintptr_t* variant_include);
896
897 PglErr ConditionalAllocateNonAutosomalVariants(const ChrInfo* cip, const char* calc_descrip, uint32_t raw_variant_ct, const uintptr_t** variant_include_ptr, uint32_t* variant_ct_ptr);
898
899 void FillSubsetChrFoVidxStart(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t* subset_chr_fo_vidx_start);
900
AllocAndFillSubsetChrFoVidxStart(const uintptr_t * variant_include,const ChrInfo * cip,uint32_t ** subset_chr_fo_vidx_start_ptr)901 HEADER_INLINE BoolErr AllocAndFillSubsetChrFoVidxStart(const uintptr_t* variant_include, const ChrInfo* cip, uint32_t** subset_chr_fo_vidx_start_ptr) {
902 const uint32_t chr_ct = cip->chr_ct;
903 if (unlikely(bigstack_alloc_u32(chr_ct + 1, subset_chr_fo_vidx_start_ptr))) {
904 return 1;
905 }
906 FillSubsetChrFoVidxStart(variant_include, cip, *subset_chr_fo_vidx_start_ptr);
907 return 0;
908 }
909
910 /*
911 // newval does not need to be null-terminated
912 // assumes *allele_ptr is not initialized
913 // (stop using these in main plink2 binary?)
914 BoolErr allele_set(const char* newval, uint32_t allele_slen, char** allele_ptr);
915
916 // *allele_ptr must be initialized; frees *allele_ptr if necessary
917 BoolErr allele_reset(const char* newval, uint32_t allele_slen, char** allele_ptr);
918
919 void cleanup_allele_storage(uint32_t max_allele_slen, uintptr_t allele_storage_entry_ct, const char** allele_storage);
920 */
921
922 CONSTI32(kMaxMissingPhenostrBlen, 32);
923 // might want g_input_missing_catname and/or g_output_missing_catname later,
924 // but let's start with the simplest implementation
925 extern char g_missing_catname[]; // default "NONE", not changeable for now
926
927 extern char g_output_missing_pheno[]; // default "NA"
928 extern char g_legacy_output_missing_pheno[]; // default "-9"
929
930 // don't care about kfUnsortedVarChrom
931 FLAGSET_DEF_START()
932 kfUnsortedVar0,
933 kfUnsortedVarBp = (1 << 0),
934 kfUnsortedVarCm = (1 << 1),
935 kfUnsortedVarSplitChr = (1 << 2)
936 FLAGSET_DEF_END(UnsortedVar);
937
938 FLAGSET_DEF_START()
939 kfFamCol0,
940 kfFamCol1 = (1 << 0),
941 kfFamCol34 = (1 << 1),
942 kfFamCol5 = (1 << 2),
943 kfFamCol6 = (1 << 3),
944 kfFamCol13456 = (kfFamCol1 | kfFamCol34 | kfFamCol5 | kfFamCol6)
945 FLAGSET_DEF_END(FamCol);
946
Sexchar(const uintptr_t * sex_nm,const uintptr_t * sex_male,uintptr_t sample_uidx)947 HEADER_INLINE char Sexchar(const uintptr_t* sex_nm, const uintptr_t* sex_male, uintptr_t sample_uidx) {
948 if (IsSet(sex_nm, sample_uidx)) {
949 return '2' - IsSet(sex_male, sample_uidx);
950 }
951 return '0';
952 }
953
954 // kPhenoDtypeCc and kPhenoDtypeQt currently can't change
955 // kPhenoDtypeOther currently used for --glm local covariates
956 ENUM_U31_DEF_START()
957 kPhenoDtypeCc,
958 kPhenoDtypeQt,
959 kPhenoDtypeCat,
960 kPhenoDtypeOther
961 ENUM_U31_DEF_END(PhenoDtype);
962
963 typedef union {
964 uintptr_t* cc; // bitvector
965 double* qt;
966 uint32_t* cat; // always 0 for missing, nonmiss[] check unnecessary
967 } PhenoData;
968
969 typedef struct PhenoColStruct {
970 // * If categorical phenotype, [0] points to g_missing_catname, while [1],
971 // [2], etc. point to category names. These are part of the same
972 // allocation as nonmiss, so no separate free is needed.
973 // Otherwise, this is nullptr.
974 // * When .sample categorical variables are imported, 'P' is added in front
975 // of the integers.
976 const char** category_names;
977
978 uintptr_t* nonmiss; // bitvector
979
980 // essentially a tagged union; part of the same allocation as nonmiss
981 PhenoData data;
982 PhenoDtype type_code;
983
984 uint32_t nonnull_category_ct;
985 } PhenoCol;
986
987 void InitPheno();
988
989
990 uint32_t IsCategoricalPhenostr(const char* phenostr_iter);
991
992 uint32_t IsCategoricalPhenostrNocsv(const char* phenostr_iter);
993
994 // returns 0xffffffffU if none exists
995 uint32_t FirstCcOrQtPhenoIdx(const PhenoCol* pheno_cols, uint32_t pheno_ct);
996
997 // "_covar" since this doesn't handle case/control
998 uint32_t IsConstCovar(const PhenoCol* covar_col, const uintptr_t* sample_include, uint32_t sample_ct);
999
1000 // returns number of nonempty categories, null included
1001 uint32_t IdentifyRemainingCats(const uintptr_t* sample_include, const PhenoCol* covar_col, uint32_t sample_ct, uintptr_t* observed_cat_bitarr);
1002
1003 // returns index of most common category
1004 uint32_t IdentifyRemainingCatsAndMostCommon(const uintptr_t* sample_include, const PhenoCol* covar_col, uint32_t sample_ct, uintptr_t* observed_cat_bitarr, uint32_t* cat_obs_buf);
1005
1006 uint32_t GetCatSamples(const uintptr_t* sample_include_base, const PhenoCol* cat_pheno_col, uint32_t raw_sample_ctl, uint32_t sample_ct, uint32_t cat_uidx, uintptr_t* cur_cat_samples);
1007
1008 // pheno_names is also allocated on the heap, but it can be handled with a
1009 // simple free_cond().
1010 void CleanupPhenoCols(uint32_t pheno_ct, PhenoCol* pheno_cols);
1011
1012 PglErr ParseChrRanges(const char* const* argvk, const char* flagname_p, const char* errstr_append, uint32_t param_ct, uint32_t allow_extra_chrs, uint32_t xymt_subtract, char range_delim, ChrInfo* cip, uintptr_t* chr_mask);
1013
1014
1015 // uint32_t MultiallelicVariantPresent(const uintptr_t* variant_include, const uintptr_t* allele_idx_offsets, uint32_t variant_ct);
1016
1017 uint32_t CountBiallelicVariants(const uintptr_t* variant_include, const uintptr_t* allele_idx_offsets, uint32_t variant_ct);
1018
1019 uintptr_t CountAlleles(const uintptr_t* variant_include, const uintptr_t* allele_idx_offsets, uint32_t raw_variant_ct, uint32_t variant_ct);
1020
1021 // This can be an overestimate (it won't return a value smaller than
1022 // read_block_size).
1023 uintptr_t GetMaxAltAlleleBlockSize(const uintptr_t* variant_include, const uintptr_t* allele_idx_offsets, uint32_t raw_variant_ct, uint32_t variant_ct, uint32_t read_block_size);
1024
1025 // Mhc stands for "multiallelic hardcall" here, not major histocompatibility
1026 // complex. Though multiallelic hardcalls are of course especially relevant
1027 // for analysis of major histocompatibility complex genomic data...
1028 uintptr_t GetMhcWordCt(uintptr_t sample_ct);
1029
1030 // sample_ct not relevant if genovecs_ptr == nullptr
1031 // only possible error is kPglRetNomem for now
1032 PglErr PgenMtLoadInit(const uintptr_t* variant_include, uint32_t sample_ct, uint32_t variant_ct, uintptr_t bytes_avail, uintptr_t pgr_alloc_cacheline_ct, uintptr_t thread_xalloc_cacheline_ct, uintptr_t per_variant_xalloc_byte_ct, uintptr_t per_alt_allele_xalloc_byte_ct, PgenFileInfo* pgfip, uint32_t* calc_thread_ct_ptr, uintptr_t*** genovecs_ptr, uintptr_t*** mhc_ptr, uintptr_t*** phasepresent_ptr, uintptr_t*** phaseinfo_ptr, uintptr_t*** dosage_present_ptr, Dosage*** dosage_mains_ptr, uintptr_t*** dphase_present_ptr, SDosage*** dphase_delta_ptr, uint32_t* read_block_size_ptr, uintptr_t* max_alt_allele_block_size_ptr, STD_ARRAY_REF(unsigned char*, 2) main_loadbufs, PgenReader*** pgr_pps, uint32_t** read_variant_uidx_starts_ptr);
1033
1034 // Returns number of variants in current block. Increases read_block_idx as
1035 // necessary (to get to a nonempty block).
1036 uint32_t MultireadNonempty(const uintptr_t* variant_include, const ThreadGroup* tgp, uint32_t raw_variant_ct, uint32_t read_block_size, PgenFileInfo* pgfip, uint32_t* read_block_idxp, PglErr* reterrp);
1037
1038 // Assumes mhc != nullptr, and is vector-aligned.
1039 void ExpandMhc(uint32_t sample_ct, uintptr_t* mhc, uintptr_t** patch_01_set_ptr, AlleleCode** patch_01_vals_ptr, uintptr_t** patch_10_set_ptr, AlleleCode** patch_10_vals_ptr);
1040
SetPgvMhc(uint32_t sample_ct,uintptr_t * mhc,PgenVariant * pgvp)1041 HEADER_INLINE void SetPgvMhc(uint32_t sample_ct, uintptr_t* mhc, PgenVariant* pgvp) {
1042 ExpandMhc(sample_ct, mhc, &(pgvp->patch_01_set), &(pgvp->patch_01_vals), &(pgvp->patch_10_set), &(pgvp->patch_10_vals));
1043 }
1044
SetPgvThreadMhcNull(uint32_t sample_ct,uint32_t tidx,uintptr_t ** thread_mhc,PgenVariant * pgvp)1045 HEADER_INLINE void SetPgvThreadMhcNull(uint32_t sample_ct, uint32_t tidx, uintptr_t** thread_mhc, PgenVariant* pgvp) {
1046 if (!thread_mhc) {
1047 pgvp->patch_01_set = nullptr;
1048 pgvp->patch_01_vals = nullptr;
1049 pgvp->patch_10_set = nullptr;
1050 pgvp->patch_10_vals = nullptr;
1051 } else {
1052 ExpandMhc(sample_ct, thread_mhc[tidx], &(pgvp->patch_01_set), &(pgvp->patch_01_vals), &(pgvp->patch_10_set), &(pgvp->patch_10_vals));
1053 }
1054 }
1055
1056
StoreStringAtBase(unsigned char * arena_top,const char * src,uintptr_t slen,unsigned char ** arena_bottom_ptr,char ** dst)1057 HEADER_INLINE BoolErr StoreStringAtBase(unsigned char* arena_top, const char* src, uintptr_t slen, unsigned char** arena_bottom_ptr, char** dst) {
1058 if (unlikely(slen >= S_CAST(uintptr_t, arena_top - (*arena_bottom_ptr)))) {
1059 return 1;
1060 }
1061 memcpyx(*arena_bottom_ptr, src, slen, '\0');
1062 *dst = R_CAST(char*, *arena_bottom_ptr);
1063 *arena_bottom_ptr += slen + 1;
1064 return 0;
1065 }
1066
StoreStringAtEnd(unsigned char * arena_bottom,const char * src,uintptr_t slen,unsigned char ** arena_top_ptr,char ** dst)1067 HEADER_INLINE BoolErr StoreStringAtEnd(unsigned char* arena_bottom, const char* src, uintptr_t slen, unsigned char** arena_top_ptr, char** dst) {
1068 // minor todo: verify that the tiny amount of additional safety provided by
1069 // PtrWSubCk has no/negligible performance cost
1070 if (PtrWSubCk(arena_bottom, slen + 1, arena_top_ptr)) {
1071 return 1;
1072 }
1073 memcpyx(*arena_top_ptr, src, slen, '\0');
1074 *dst = R_CAST(char*, *arena_top_ptr);
1075 return 0;
1076 }
1077
StoreStringAtEndK(unsigned char * arena_bottom,const char * src,uintptr_t slen,unsigned char ** arena_top_ptr,const char ** dst)1078 HEADER_INLINE BoolErr StoreStringAtEndK(unsigned char* arena_bottom, const char* src, uintptr_t slen, unsigned char** arena_top_ptr, const char** dst) {
1079 if (PtrWSubCk(arena_bottom, slen + 1, arena_top_ptr)) {
1080 return 1;
1081 }
1082 memcpyx(*arena_top_ptr, src, slen, '\0');
1083 *dst = R_CAST(char*, *arena_top_ptr);
1084 return 0;
1085 }
1086
StoreStringAndPrecharAtEnd(unsigned char * arena_bottom,const char * src,unsigned char prechar,uintptr_t slen,unsigned char ** arena_top_ptr,char ** dst)1087 HEADER_INLINE BoolErr StoreStringAndPrecharAtEnd(unsigned char* arena_bottom, const char* src, unsigned char prechar, uintptr_t slen, unsigned char** arena_top_ptr, char** dst) {
1088 if (PtrWSubCk(arena_bottom, slen + 2, arena_top_ptr)) {
1089 return 1;
1090 }
1091 **arena_top_ptr = prechar;
1092 char* dst_write = 1 + R_CAST(char*, *arena_top_ptr);
1093 memcpyx(dst_write, src, slen, '\0');
1094 *dst = dst_write;
1095 return 0;
1096 }
1097
1098 // These use g_textbuf.
1099 PglErr WriteSampleIdsOverride(const uintptr_t* sample_include, const SampleIdInfo* siip, const char* outname, uint32_t sample_ct, SampleIdFlags override_flags);
WriteSampleIds(const uintptr_t * sample_include,const SampleIdInfo * siip,const char * outname,uint32_t sample_ct)1100 HEADER_INLINE PglErr WriteSampleIds(const uintptr_t* sample_include, const SampleIdInfo* siip, const char* outname, uint32_t sample_ct) {
1101 return WriteSampleIdsOverride(sample_include, siip, outname, sample_ct, siip->flags);
1102 }
1103
1104 // read_realpath must be a buffer of size >= kPglFnamesize bytes
1105 uint32_t RealpathIdentical(const char* outname, const char* read_realpath, char* write_realpath_buf);
1106
1107 char* PrintHaploidNonintDosage(uint32_t rawval, char* start);
1108
1109 char* PrintMultiallelicHcAsDs(uint32_t hc1, uint32_t hc2, uint32_t allele_ct, char* start);
1110
1111 char* PrintMultiallelicHcAsHaploidDs(uint32_t hc1, uint32_t hc2, uint32_t allele_ct, char* start);
1112
1113
OutnameZstSet(const char * ext,uint32_t output_zst,char * outname_end)1114 HEADER_INLINE void OutnameZstSet(const char* ext, uint32_t output_zst, char* outname_end) {
1115 const uint32_t ext_slen = strlen(ext);
1116 assert(ext_slen < kMaxOutfnameExtBlen - 4);
1117 memcpy(outname_end, ext, ext_slen + 1);
1118 if (output_zst) {
1119 strcpy_k(&(outname_end[ext_slen]), ".zst");
1120 }
1121 }
1122
1123 ENUM_U31_DEF_START()
1124 kVfilterExtract,
1125 kVfilterExtractIntersect,
1126 kVfilterExclude
1127 ENUM_U31_DEF_END(VfilterType);
1128
1129 extern const char g_vft_names[3][18];
1130
CleanupPgfi2(const char * file_descrip,PgenFileInfo * pgfip,PglErr * reterrp)1131 HEADER_INLINE BoolErr CleanupPgfi2(const char* file_descrip, PgenFileInfo* pgfip, PglErr* reterrp) {
1132 if (CleanupPgfi(pgfip, reterrp)) {
1133 logerrprintfww(kErrprintfFread, file_descrip, strerror(errno));
1134 return 1;
1135 }
1136 return 0;
1137 }
1138
CleanupPgr2(const char * file_descrip,PgenReader * pgrp,PglErr * reterrp)1139 HEADER_INLINE BoolErr CleanupPgr2(const char* file_descrip, PgenReader* pgrp, PglErr* reterrp) {
1140 if (CleanupPgr(pgrp, reterrp)) {
1141 logerrprintfww(kErrprintfFread, file_descrip, strerror(errno));
1142 return 1;
1143 }
1144 return 0;
1145 }
1146
PgenErrPrintNEx(const char * file_descrip,PglErr reterr)1147 HEADER_INLINE void PgenErrPrintNEx(const char* file_descrip, PglErr reterr) {
1148 if (reterr == kPglRetReadFail) {
1149 logputs("\n");
1150 logerrprintfww(kErrprintfFread, file_descrip, strerror(errno));
1151 } else if (reterr == kPglRetMalformedInput) {
1152 logputs("\n");
1153 logerrprintfww("Error: Malformed %s.\n", file_descrip);
1154 }
1155 }
1156
PgenErrPrintN(PglErr reterr)1157 HEADER_INLINE void PgenErrPrintN(PglErr reterr) {
1158 PgenErrPrintNEx(".pgen file", reterr);
1159 }
1160
PgenErrPrintEx(const char * file_descrip,PglErr reterr)1161 HEADER_INLINE void PgenErrPrintEx(const char* file_descrip, PglErr reterr) {
1162 if (reterr == kPglRetReadFail) {
1163 logerrprintfww(kErrprintfFread, file_descrip, strerror(errno));
1164 } else if (reterr == kPglRetMalformedInput) {
1165 logerrprintfww("Error: Malformed %s.\n", file_descrip);
1166 }
1167 }
1168
PgenErrPrint(PglErr reterr)1169 HEADER_INLINE void PgenErrPrint(PglErr reterr) {
1170 PgenErrPrintEx(".pgen file", reterr);
1171 }
1172
1173 #ifdef __cplusplus
1174 } // namespace plink2
1175 #endif
1176
1177 #endif // __PLINK2_COMMON_H__
1178