1 /* $Id: thrdatd.c 461924 2015-03-13 16:44:08Z vasilche $
2 *===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * File Name: thrdatd.c
27 *
28 * Author: Stephen Bryant
29 *
30 * Initial Version Creation Date: 08/16/2000
31 *
32 * $Revision: 461924 $
33 *
34 * File Description: threader
35 */
36
37 #ifdef _MSC_VER
38 #pragma warning(disable:4244) // disable double->float warning in MSVC
39 #pragma warning(disable:4018) // disable signed/unsigned mismatch warning in MSVC
40 #endif
41
42 /*----------------------------------------------------------------------------*/
43 /* atd - adaptive threading of protein sequence through folding motif */
44 /* Steve Bryant, 2/94 */
45 /*----------------------------------------------------------------------------*/
46 /* Uses the Gibb's sampling algorithm to find optimal alignments of a */
47 /* query sequence and folding motif, and to locate the optimal core */
48 /* structure, given a set of constraints on the locations of each core */
49 /* segment. */
50 /*----------------------------------------------------------------------------*/
51 /* Test version of 8/96: Sets alignment-only flag to ttbi.c. */
52 /* Also passes trajectory data back to calling routine, satd.c. */
53 /* Also has test code to stop iterations in apparent local minima */
54 /*----------------------------------------------------------------------------*/
55
56 #include <stdlib.h>
57 #include <math.h>
58 #include <stdio.h>
59 #include <string.h>
60
61 #include <algo/structure/threader/thrdatd.h>
62 #include <algo/structure/threader/thrddecl.h>
63
64 typedef char Boolean;
65 #define TRUE 1
66 #define FALSE 0
67
atd(Fld_Mtf * mtf,Cor_Def * cdf,Qry_Seq * qsq,Rcx_Ptl * pmf,Gib_Scd * gsp,Thd_Tbl * ttb,Seq_Mtf * psm,float * trg,int zscs,double ScalingFactor,float PSSM_Weight)68 int atd(Fld_Mtf* mtf, Cor_Def* cdf, Qry_Seq* qsq, Rcx_Ptl* pmf,
69 Gib_Scd* gsp, Thd_Tbl* ttb, Seq_Mtf* psm, float* trg, int zscs,
70 double ScalingFactor, float PSSM_Weight) {
71 /*------------------------------------------------------------------*/
72 /* mtf: Contact matrices defining the folding motif */
73 /* cdf: Core segment locations and loop length limits */
74 /* qsq: Sequence to thread with alignment contraints */
75 /* pmf: Potential of mean force as a 3-d lookup table */
76 /* gsp: Various control parameters for Gibb's sampling */
77 /* ttb: Tables to hold Results of Gibbs sampled threading */
78 /* psm: Motif matching parameters */
79 /* trg: Will store trajectory energy values */
80 /* zscs: which z scores calculated? 0=none, 1=lowest energy, 2=all */
81 /*------------------------------------------------------------------*/
82
83 Cur_Loc *sli; /* Current locations of core segments in the motif */
84 Cur_Aln *sai; /* Current alignment of query sequence with core */
85 Seg_Nsm *spn; /* Partial sums of contact counts by segment pair */
86 Seg_Cmp *spc; /* Partial sums of residue counts by segment pair */
87 Thd_Cxe *cxe; /* Expected energy of a contact in current thread */
88 Seg_Gsm *spe; /* Partial sums of contact energies by segment pair */
89 Thd_Gsm *tdg; /* Total energies for the current threaded sequence */
90 Cxl_Los **cpr; /* Contacts by segment, largest possible set */
91 Cxl_Los **cpl; /* Contacts by segment, given current alignment */
92 Cxl_Als **cpa; /* Contacts by segment, given current location */
93 Seg_Ord *sgo; /* Segment samping order and related quantities */
94 Rnd_Smp *pvl; /* Alignment-location parameter value probabilities */
95 Thd_Tst *tts; /* Parameters for local min and convergence tests */
96
97 /*----------------------------------------------------------------------------*/
98 /* Parameters from argument data structures */
99 /*----------------------------------------------------------------------------*/
100 int nsc; /* Number of threaded segments in core definition */
101 int nmt; /* Number of residue positions in the folding motif */
102 int nlp; /* Number of loops, including n- and c-terminal tails */
103 /*int nrr;*/ /* Number of residue-residue contacts in motif */
104 /*int nrp;*/ /* Number of residue-peptide contacts in motif */
105 int nrt; /* Number of residue types in contact potential */
106 int ndi; /* Number of distance intervals in potential */
107
108 /*----------------------------------------------------------------------------*/
109 /* Countdown timers for Gibbs annealing schedule */
110 /*----------------------------------------------------------------------------*/
111 int nrs; /* Number of random starts for Gibb's sampling */
112 int nts; /* Number of temperature steps */
113 int nti; /* Number of iterations per tempeature step */
114 int nac; /* Number of alignment cycles per iteration */
115 int nlc; /* Number of extent/location cycles per iteration */
116 int ngs; /* Number of trajectory points recorded */
117
118 /*----------------------------------------------------------------------------*/
119 /* Other local variables */
120 /*----------------------------------------------------------------------------*/
121 int i,j,k,n; /* Counters */
122 Cxl_Als *ca; /* Pointer to alignment sampling pair lists */
123 Cxl_Los *cl; /* Pointer to location sampling lists */
124 Cxl_Los *cr; /* Pointer to contact pair lists */
125 int spcd; /* Flags a change in threaded sequence composition */
126 int mx,mn; /* Range of segment alignment or location */
127 int cs; /* Current segment in alignment or location sampling */
128 int ct=0; /* Current terminus for segment location sampling */
129 int al; /* Current alignment of a core segment */
130 int of; /* Current terminus offset from reference position */
131 /*int rf;*/ /* Reference position for a core element */
132 int tmp; /* temp holder for ttb->mx */
133 int dist, dist2; /* for loop distances */
134 /*float hh;*/
135
136 /*----------------------------------------------------------------------------*/
137 /* Function declarations for routines returning non-integer values */
138 /*----------------------------------------------------------------------------*/
139 /* float dgri(); */
140
141
142 /*----------------------------------------------------------------------------*/
143 /* if 2 aligned regions on the query sequence are constrained to match */
144 /* 2 adjacent core regions, then make sure the loop constraints are */
145 /* sufficiently big. */
146 /* this over-rides the automatically determined loop constraints. */
147 /*----------------------------------------------------------------------------*/
148 /* for each loop between core blocks */
149 for (i=0; i<(cdf->sll.n-1); i++) {
150 /* if 2 aligned blocks on query sequence are constrained to adjacent core blocks */
151 if ((qsq->sac.mn[i] > 0) && (qsq->sac.mn[i+1] > 0)) {
152 /* calculate loop distance of query sequence */
153 dist = (qsq->sac.mn[i+1] - cdf->sll.nomn[i+1]) - (qsq->sac.mn[i] + cdf->sll.comn[i]) - 1;
154 /* if 1.5 times this dist is greater than the loop constraint, then relax the loop constraint */
155 dist = (int) ((double)dist * 1.5);
156 if (dist > cdf->lll.llmx[i+1]) {
157 cdf->lll.llmx[i+1] = dist;
158 }
159 }
160 }
161
162
163 /*----------------------------------------------------------------------------*/
164 /* if 2 adjacent aligned regions on the query sequence are constrained to */
165 /* match 2 core regions that are separated by one intervening core region, */
166 /* then make sure the loop constraints are sufficiently big. */
167 /* this over-rides the automatically determined loop constraints. */
168 /*----------------------------------------------------------------------------*/
169 /* for each pair of core blocks separated by a core block */
170 for (i=0; i<(cdf->sll.n-2); i++) {
171 /* if 2 adjacent aligned blocks on the query sequence are constrained */
172 /* to 2 core blocks separated by an intervening core block */
173 if ((qsq->sac.mn[i] > 0) && (qsq->sac.mn[i+1] <= 0) && (qsq->sac.mn[i+2] > 0)) {
174 /* calculate loop distance of query sequence */
175 dist = (qsq->sac.mn[i+2] - cdf->sll.nomn[i+2]) - (qsq->sac.mn[i] + cdf->sll.comn[i]) - 1;
176 /* calculate length of intervening aligned block */
177 dist2 = cdf->sll.nomn[i+1] + cdf->sll.comn[i+1] + 1;
178 /* relax loop constraints between 2 core blocks and intervening core block */
179 if ((dist - dist2) > cdf->lll.llmx[i+1]) {
180 cdf->lll.llmx[i+1] = dist - dist2;
181 }
182 if ((dist - dist2) > cdf->lll.llmx[i+2]) {
183 cdf->lll.llmx[i+2] = dist - dist2;
184 }
185 }
186 }
187
188
189 /*----------------------------------------------------------------------------*/
190 /* Retrieve parameters */
191 /*----------------------------------------------------------------------------*/
192 nsc=cdf->sll.n;
193 nmt=mtf->n;
194 nlp=cdf->lll.n;
195 /*nrr=mtf->rrc.n;*/
196 /*nrp=mtf->rpc.n;*/
197 nrt=pmf->nrt;
198 ndi=pmf->ndi;
199
200 #ifdef ATD_DEBUG
201 printf("nsc %d\n",nsc);
202 printf("nmt %d\n",nmt);
203 printf("nlp %d\n",nlp);
204 /*printf("nrr %d\n",nrr);*/
205 /*printf("nrp %d\n",nrp);*/
206 printf("nrt %d\n",nrt);
207 printf("ndi %d\n",ndi);
208 #endif
209
210
211 /*----------------------------------------------------------------------------*/
212 /* Allocate storage for threading data structures */
213 /*----------------------------------------------------------------------------*/
214 sli=calloc(1,sizeof(Cur_Loc));
215 sai=calloc(1,sizeof(Cur_Aln));
216 spn=calloc(1,sizeof(Seg_Nsm));
217 spc=calloc(1,sizeof(Seg_Cmp));
218 cxe=calloc(1,sizeof(Thd_Cxe));
219 spe=calloc(1,sizeof(Seg_Gsm));
220 tdg=calloc(1,sizeof(Thd_Gsm));
221 cpr=calloc(nsc,sizeof(Cxl_Los *));
222 for(i=0;i<nsc;i++) cpr[i]=calloc(1,sizeof(Cxl_Los));
223 cpa=calloc(nsc,sizeof(Cxl_Als *));
224 for(i=0;i<nsc;i++) cpa[i]=calloc(1,sizeof(Cxl_Als));
225 cpl=calloc(nsc,sizeof(Cxl_Los *));
226 for(i=0;i<nsc;i++) cpl[i]=calloc(1,sizeof(Cxl_Los));
227 sgo=calloc(1,sizeof(Seg_Ord));
228 pvl=calloc(1,sizeof(Rnd_Smp));
229 tts=calloc(1,sizeof(Thd_Tst));
230
231 #ifdef ATD_DEBUG
232 printf("allocated threading data structures\n");
233 #endif
234
235
236 /*----------------------------------------------------------------------------*/
237 /* Allocate storage for current segment locations, type Cur_Loc */
238 /*----------------------------------------------------------------------------*/
239 sli->nsc=nsc;
240 sli->nlp=nlp;
241 sli->nmt=nmt;
242 sli->no=(int *)calloc(nsc,sizeof(int));
243 sli->co=(int *)calloc(nsc,sizeof(int));
244 sli->lp=(int *)calloc(nsc+1,sizeof(int));
245 sli->cr=(int *)calloc(nmt,sizeof(int));
246
247 #ifdef ATD_DEBUG
248 printf("sli->nsc:%d nlp:%d nmt%d\n",sli->nsc,sli->nlp,sli->nmt);
249 for(i=0;i<nsc;i++) printf("%d ",sli->no[i]); printf("sli->no\n");
250 for(i=0;i<nsc;i++) printf("%d ",sli->co[i]); printf("sli->co\n");
251 for(i=0;i<nlp;i++) printf("%d ",sli->lp[i]); printf("sli->lp\n");
252 for(i=0;i<nmt;i++) printf("%d ",sli->cr[i]); printf("sli->cr\n");
253 #endif
254
255 /* Allocate storage for current segment alignment and sequence, type Cur_Aln */
256 sai->nsc=nsc;
257 sai->nmt=nmt;
258 sai->al=(int *)calloc(nsc,sizeof(int));
259 sai->sq=(int *)calloc(nmt,sizeof(int));
260 sai->cf=(int *)calloc(nmt,sizeof(int));
261
262 #ifdef ATD_DEBUG
263 printf("sai->nsc:%d nmt%d\n",sai->nsc,sai->nmt);
264 for(i=0;i<nsc;i++) printf("%d ",sai->al[i]); printf("sai->al\n");
265 for(i=0;i<nmt;i++) printf("%d ",sai->sq[i]); printf("sai->sq\n");
266 #endif
267
268
269 /* Allocate storage for current segment pair energies, type Seg_Gsm */
270 spe->nsc=nsc;
271 spe->gss=calloc(nsc,sizeof(int *));
272 for(i=0;i<nsc;i++) spe->gss[i]=calloc(nsc,sizeof(int));
273 spe->gs=calloc(nsc,sizeof(int));
274 spe->ms=calloc(nsc,sizeof(int));
275 spe->cs=calloc(nsc,sizeof(int));
276 spe->ls=calloc(nsc,sizeof(int));
277 spe->s0=calloc(nsc,sizeof(int));
278
279 #ifdef ATD_DEBUG
280 for(i=0; i<nsc; i++) printf("%d ",spe->gs[i]); printf("spe->gs\n");
281 for(j=0; j<nsc; j++) {
282 for(i=0; i<nsc; i++) printf("%d ",spe->gss[j][i]);
283 printf("spe->gss[%d]\n",j); }
284 #endif
285
286 /* Allocate storage for current segment pair contact counts, type Seg_Nsm */
287 spn->ndi=ndi;
288 spn->nsc=nsc;
289 spn->nrt=nrt;
290 spn->nrr=calloc(ndi,sizeof(int **));
291 for(i=0; i<ndi; i++) spn->nrr[i]=calloc(nsc,sizeof(int *));
292 for(i=0; i<ndi; i++) for (j=0; j<nsc; j++)
293 spn->nrr[i][j]=calloc(nsc,sizeof(int));
294 spn->srr=calloc(ndi,sizeof(int));
295 spn->nrp=calloc(ndi,sizeof(int **));
296 for(i=0; i<ndi; i++) spn->nrp[i]=calloc(nsc,sizeof(int *));
297 for(i=0; i<ndi; i++) for (j=0; j<nsc; j++)
298 spn->nrp[i][j]=calloc(nsc,sizeof(int));
299 spn->srp=calloc(ndi,sizeof(int));
300 spn->nrf=calloc(ndi,sizeof(int **));
301 for(i=0; i<ndi; i++) spn->nrf[i]=calloc(nsc,sizeof(int *));
302 for(i=0; i<ndi; i++) for (j=0; j<nsc; j++)
303 spn->nrf[i][j]=calloc(nrt,sizeof(int));
304 spn->frf=calloc(ndi,sizeof(int *));
305 for(i=0;i<ndi;i++) spn->frf[i]=calloc(nrt,sizeof(int));
306 spn->srf=calloc(ndi,sizeof(int));
307
308 #ifdef ATD_DEBUG
309 printf("spn->ndi:%d nsc:%d nrt%d\n",spn->ndi,spn->nsc,spn->nrt);
310 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
311 for(i=0; i<nsc; i++) printf("%d ",spn->nrr[k][j][i]);
312 printf("spn->nrr[%d][%d]\n",k,j); }
313 for(i=0; i<ndi; i++) printf("%d ",spn->srr[i]); printf("spn->srr\n");
314 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
315 for(i=0; i<nsc; i++) printf("%d ",spn->nrp[k][j][i]);
316 printf("spn->nrp[%d][%d]\n",k,j); }
317 for(i=0; i<ndi; i++) printf("%d ",spn->srp[i]); printf("spn->srp\n");
318 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
319 for(i=0; i<nrt; i++) printf("%d ",spn->nrf[k][j][i]);
320 printf("spn->nrf[%d][%d]\n",k,j); }
321 for(j=0; j<ndi; j++) {
322 for(i=0;i<nrt;i++) printf("%d ",spn->frf[j][i]);
323 printf("spn->frf[%d]\n",j); }
324 for(i=0; i<ndi; i++) printf("%d ",spn->srf[i]); printf("spn->srf\n");
325 #endif
326
327 /* Allocate storage for threaded segment compostion, type Seg_Cmp */
328 spc->nsc=nsc;
329 spc->nrt=nrt;
330 spc->srt=calloc(nsc,sizeof(int *));
331 for(i=0; i<nsc; i++) spc->srt[i]=calloc(nrt,sizeof(int));
332 spc->nlp=nlp;
333 spc->lrt=calloc(nlp,sizeof(int *));
334 for(i=0; i<nlp; i++) spc->lrt[i]=calloc(nrt,sizeof(int));
335 spc->rt=calloc(nrt,sizeof(int));
336 spc->rto=calloc(nrt,sizeof(int));
337
338 #ifdef ATD_DEBUG
339 printf("spc->nsc:%d nrt:%d nlp:%d\n",spc->nsc,spc->nrt,spc->nlp);
340 for(j=0; j<nsc; j++) {
341 for(i=0;i<nrt;i++) printf("%d ",spc->srt[j][i]);
342 printf("spc->srt[%d]\n",j); }
343 for(j=0; j<nlp; j++) {
344 for(i=0;i<nrt;i++) printf("%d ",spc->lrt[j][i]);
345 printf("spc->lrt[%d]\n",j); }
346 for(i=0; i<nrt; i++) printf("%d ",spc->rt[i]); printf("spc->rt\n");
347 for(i=0; i<nrt; i++) printf("%d ",spc->rto[i]); printf("spc->rto\n");
348 #endif
349
350 /* Allocate storage for expected contact energies, type Thd_Cxe */
351 cxe->nrt=nrt;
352 cxe->ndi=ndi;
353 cxe->rp=calloc(nrt,sizeof(float));
354 cxe->rrp=calloc(nrt,sizeof(float *));
355 for(i=0; i<nrt; i++) cxe->rrp[i]=calloc(i+1,sizeof(float));
356 cxe->rre=calloc(ndi,sizeof(float));
357 cxe->rpe=calloc(ndi,sizeof(float));
358 cxe->rfe=calloc(ndi,sizeof(float));
359
360 #ifdef ATD_DEBUG
361 printf("cxe->nrt:%d ndi:%d\n",cxe->nrt,cxe->ndi);
362 for(i=0; i<nrt; i++) printf("%.4f ",cxe->rp[i]); printf("cxe->rp\n");
363 for(j=0; j<nrt; j++) {
364 for(i=0;i<nrt;i++) printf("%.4f ",cxe->rrp[j][i]);
365 printf("cxe->rrp[%d]\n",j); }
366 for(i=0; i<ndi; i++) printf("%.4f ",cxe->rre[i]); printf("cxe->rre\n");
367 for(i=0; i<ndi; i++) printf("%.4f ",cxe->rpe[i]); printf("cxe->rpe\n");
368 for(i=0; i<ndi; i++) printf("%.4f ",cxe->rfe[i]); printf("cxe->rfe\n");
369 #endif
370
371 /* Allocate storage for segment sampling order, type Seg_Ord */
372 sgo->nsc=nsc;
373 sgo->si=(int *)calloc(nsc,sizeof(int));
374 sgo->so=(int *)calloc(nsc,sizeof(int));
375 sgo->to=(int *)calloc(nsc,sizeof(int));
376
377 #ifdef ATD_DEBUG
378 printf("sgo->nsc:%d\n",sgo->nsc);
379 for(i=0; i<nsc; i++) printf("%d ",sgo->si[i]); printf("sgo->si\n");
380 for(i=0; i<nsc; i++) printf("%d ",sgo->so[i]); printf("sgo->so\n");
381 for(i=0; i<nsc; i++) printf("%d ",sgo->to[i]); printf("sgo->to\n");
382 #endif
383
384
385 /* Allocate storage for random sampling, type Rnd_Smp */
386 /* Need at most one value per residue in the query sequence */
387 pvl->p=(float *)calloc(qsq->n,sizeof(float));
388 pvl->e=(float *)calloc(qsq->n,sizeof(float));
389 pvl->sq=(int *)calloc(qsq->n,sizeof(int));
390 pvl->aqi=(int *)calloc(qsq->n,sizeof(int));
391 pvl->r=(int *)calloc(qsq->n,sizeof(int));
392 pvl->o=(int *)calloc(qsq->n,sizeof(int));
393 pvl->lsg=(int *)calloc(nsc+1,sizeof(int));
394
395
396 /* Allocate storage for local minima and convergence tests */
397 tts->nw=ttb->n;
398 tts->bw=(float *)calloc(tts->nw,sizeof(float)); /* For Boltzmann weights */
399 tts->nb=0; for(i=0; i<gsp->nts; i++) {
400 if(tts->nb<gsp->nti[i]) tts->nb=gsp->nti[i];
401 #ifdef ATD_DEBUG
402 printf("nts=%d nti=%d tts->nb=%d\n",i,gsp->nti[i],tts->nb);
403 #endif
404 }
405 tts->ib=(float *)calloc(tts->nb,sizeof(float)); /* Best energy per iteration */
406
407
408 /* For each core segment identify all possible contact pairs */
409 cprl(mtf,cdf,pmf,cpr,0);
410
411 #ifdef ATD_DEBUG
412 for(i=0; i<nsc; i++) {
413 cr=cpr[i];
414 printf("rr:%d rp:%d rf:%d cpr[%d]\n",cr->rr.n,cr->rp.n,cr->rf.n,i); }
415 #endif
416
417
418 /* Allocate storage for master contact lists, type Cxl_Los */
419 for(i=0; i<nsc; i++) {
420 cr=cpr[i];
421 cr->rr.r1=calloc(cr->rr.n,sizeof(int));
422 cr->rr.r2=calloc(cr->rr.n,sizeof(int));
423 cr->rr.d=calloc(cr->rr.n,sizeof(int));
424 cr->rp.r1=calloc(cr->rp.n,sizeof(int));
425 cr->rp.p2=calloc(cr->rp.n,sizeof(int));
426 cr->rp.d=calloc(cr->rp.n,sizeof(int));
427 if(cr->rf.n>0) {cr->rf.r1=calloc(cr->rf.n,sizeof(int));
428 cr->rf.t2=calloc(cr->rf.n,sizeof(int));
429 cr->rf.d=calloc(cr->rf.n,sizeof(int)); } }
430
431 /* Construct master contact lists */
432 cprl(mtf,cdf,pmf,cpr,1);
433
434 #ifdef ATD_DEBUG
435 for(i=0; i<nsc; i++) {
436 cr=cpr[i];
437 printf("rr:%d rp:%d rf:%d cpr[%d]\n",cr->rr.n,cr->rp.n,cr->rf.n,i); }
438 #endif
439
440 /* Allocate storage for extent-sampling contact lists, type Cxl_Los */
441 for(i=0; i<nsc; i++) {
442 cl=cpl[i]; cr=cpr[i];
443 cl->rr.r1=calloc(cr->rr.n,sizeof(int));
444 cl->rr.r2=calloc(cr->rr.n,sizeof(int));
445 cl->rr.d=calloc(cr->rr.n,sizeof(int));
446 cl->rr.e=calloc(cr->rr.n,sizeof(int));
447 cl->rp.r1=calloc(cr->rp.n,sizeof(int));
448 cl->rp.p2=calloc(cr->rp.n,sizeof(int));
449 cl->rp.d=calloc(cr->rp.n,sizeof(int));
450 cl->rp.e=calloc(cr->rp.n,sizeof(int));
451 cl->rf.r1=calloc(cr->rf.n,sizeof(int));
452 cl->rf.t2=calloc(cr->rf.n,sizeof(int));
453 cl->rf.d=calloc(cr->rf.n,sizeof(int));
454 cl->rf.e=calloc(cr->rf.n,sizeof(int)); }
455
456
457 /* Allocate storage for alignment-sampling contact lists, type Cxl_Als */
458 for(i=0; i<nsc; i++) {
459 ca=cpa[i];
460 n=cpr[i]->rr.n;
461 ca->rr.r1=calloc(n,sizeof(int));
462 ca->rr.r2=calloc(n,sizeof(int));
463 ca->rr.e=calloc(n,sizeof(int **));
464 n=cdf->sll.nomx[i]+cdf->sll.comx[i]+1;
465 ca->r.r=calloc(n,sizeof(int));
466 ca->r.e=calloc(n,sizeof(int *));
467 for(j=0;j<n;j++) ca->r.e[j]=calloc(nrt,sizeof(int));}
468
469
470 /* Initialize pseudo-random number generation. */
471 /* srand48((long)gsp->rsd); */
472 /* RandomSeed((long)gsp->rsd); */
473 Rand01(&gsp->rsd);
474
475 /* Initialize linked list for storage of top threads */
476 ttb0(ttb);
477
478
479 /* Initialize trajectory record counter */
480 ngs=0;
481
482 #ifdef ATD_DEBUG
483 printf("ttb->n: %d mx: %d mn: %d\n",ttb->n,ttb->mx,ttb->mn);
484 for(i=0; i<ttb->n; i++) printf("i:%d tg:%.4f tf:%d pr:%d nx:%d\n",
485 i,ttb->tg[i],ttb->tf[i],ttb->pr[i],ttb->nx[i]);
486 #endif
487
488 /* ----------------------------------------------------------------- */
489 /* Begin sampling. Loop over the number of random starts specified. */
490 /* ----------------------------------------------------------------- */
491
492 for(nrs=0;nrs<gsp->nrs;nrs++) {
493 #ifdef ATD_DEBUG
494 printf("nrs:%d\n",nrs);
495 #endif
496
497 /* Choose the initial location and extent of core elements */
498
499 /* Flag all core segments as missing extent assignment */
500 for(i=0; i<nsc; i++) { sli->no[i]=-1; sli->co[i]=-1; }
501
502 /* Flag all core residue positions as not assigned to any core segment */
503 for(i=0; i<nmt; i++) sli->cr[i]=-1;
504
505 #ifdef ATD_DEBUG
506 for(i=0;i<nsc;i++) printf("%d ",sli->no[i]); printf("sli->no\n");
507 for(i=0;i<nsc;i++) printf("%d ",sli->co[i]); printf("sli->co\n");
508 for(i=0;i<nlp;i++) printf("%d ",sli->lp[i]); printf("sli->lp\n");
509 for(i=0;i<nmt;i++) printf("%d ",sli->cr[i]); printf("sli->cr\n");
510 #endif
511
512 /* Determine the order of segment location. */
513 sgoi(gsp->iso,gsp->ito,pvl,sgo);
514 /* sgoi(1,2,pvl,sgo); */
515
516 #ifdef ATD_DEBUG
517 printf("gsp->iso: %d gsp->ito: %d\n",gsp->iso,gsp->ito);
518 for(i=0;i<nsc;i++) printf("%d ",sgo->si[i]); printf("sgo->si\n");
519 for(i=0;i<nsc;i++) printf("%d ",sgo->so[i]); printf("sgo->so\n");
520 for(i=0;i<nsc;i++) printf("%d ",sgo->to[i]); printf("sgo->to\n");
521 #endif
522
523 /* Loop over core segments and terminii */
524 for(i=0; i<nsc; i++) {
525 cs=sgo->si[i];
526 /*rf=cdf->sll.rfpt[cs];*/
527 #ifdef ATD_DEBUG
528 /*printf("cs:%d rf:%d\n",cs,rf);*/
529 #endif
530 for(j=0;j<=1;j++) {
531
532 switch(j) { case 0: { ct=sgo->to[cs];
533 break; }
534 case 1: { ct=(sgo->to[cs]==0) ? 1:0;
535 break; }
536 }
537
538 /* Find allowed extent range for current terminus */
539 /* printf("ct:%d\n",ct); */
540 if(!slo0(mtf,cdf,qsq,sli,cs,ct,&mn,&mx)) {
541 printf("slo0 failed at nrs:%d cs:%d ct:%d\n",nrs,cs,ct);
542 return(0); }
543 /* printf("cs:%d ct:%d mn:%d mx:%d\n",cs,ct,mn,mx); */
544
545 /* Assign extent, using the method specified */
546 switch(gsp->isl) {
547
548 case 1: { /* Assign minimum extent */
549 switch(ct) {
550 case 0: {
551 sli->no[cs]=cdf->sll.nomn[cs];
552 break; }
553 case 1: {
554 sli->co[cs]=cdf->sll.comn[cs];
555 break;}
556 }
557 break;
558 }
559
560 default: { /* Choose extent at random */
561 switch(ct) {
562 case 0: {
563 pvl->n=mx-mn+1;
564 for(k=0;k<pvl->n;k++)
565 pvl->p[k]=1./pvl->n;
566 sli->no[cs]=mn+rsmp(pvl);
567 break; }
568 case 1: {
569 pvl->n=mx-mn+1;
570 for(k=0;k<pvl->n;k++)
571 pvl->p[k]=1./pvl->n;
572 sli->co[cs]=mn+rsmp(pvl);
573 break; }
574 }
575 break;
576 }
577
578 case 2: { /* Choose maximum extent */
579 switch(ct) {
580 case 0: {
581
582 /* disable this test since the assignment algorithm is based on alignment of
583 minimum extent core elements */ /* if(mx<cdf->sll.nomx[cs]) {
584 printf( "over maximum extent nrs:%d cs:%d ct:0\n",nrs,cs);
585 return(0); } */
586
587 sli->no[cs]=cdf->sll.nomx[cs];
588 break; }
589 case 1: {
590
591 /* if(mx<cdf->sll.comx[cs]) {
592 printf("over maximum extent nrs:%d cs:%d ct:1\n",nrs,cs);
593 return(0); } */
594
595 sli->co[cs]=cdf->sll.comx[cs];
596 break; }
597 }
598 break;
599 }
600 }
601
602 }
603 }
604
605 /* for(k=0;k<nsc;k++) printf("%d ",sli->no[k]);printf("sli->no\n");
606 for(k=0;k<nsc;k++) printf("%d ",sli->co[k]);printf("sli->co\n");*/
607
608 /* Store the initial location of the core segments */
609 for(i=0;i<nsc;i++) {
610 mn=cdf->sll.rfpt[i]-sli->no[i];
611 mx=cdf->sll.rfpt[i]+sli->co[i];
612 for(j=mn; j<=mx; j++) sli->cr[j]=i;
613 }
614
615 /* for(i=0;i<nmt;i++) printf("%d ",sli->cr[i]); printf("sli->cr\n"); */
616
617 /* Store the initial motif-derived loop lengths */
618 for(i=1;i<nsc;i++) {
619 mn=cdf->sll.rfpt[i-1]+sli->co[i-1];
620 mx=cdf->sll.rfpt[i]-sli->no[i];
621 sli->lp[i]=mtf->mll[mn][mx]; }
622 /* for(i=0;i<nlp;i++) printf("%d ",sli->lp[i]); printf("sli->lp\n"); */
623
624 /* Choose the initial alignment of the query sequence and core motif */
625
626 /* Flag all core segments as unaligned */
627 for(i=0; i<nsc; i++) sai->al[i]=-1;
628
629 /* Determine the order of segment alignment. */
630 sgoi(gsp->iso,gsp->ito,pvl,sgo);
631
632 /* Randomly align core segments in the order determined */
633 for(i=0; i<nsc; i++) { cs=sgo->si[i];
634
635 /* Find alignment constraints arising from previously aligned */
636 /* core segment, explicitly constrained core segments, and/or */
637 /* consideration of the lengths of the query sequence and core */
638 /* segments and length ranges for loops. */
639 if(!sal0(cdf,qsq,sli,sai,cs,&mn,&mx)) {
640 printf("failed sal0 nrs:%d cs:%d\n",nrs,cs);
641 return(0);}
642 /* printf("cs:%d mn:%d mx:%d\n",cs,mn,mx); */
643
644 /* Choose an alignment at random, given the min/max constraints */
645 pvl->n=mx-mn+1;
646 for(j=0;j<pvl->n;j++) pvl->p[j]=1./pvl->n;
647 sai->al[cs]=mn+rsmp(pvl);
648 }
649 /* for(k=0;k<nsc;k++) printf("%d ",sai->al[k]);printf("sai->al\n"); */
650
651 /* Assign the threaded sequence, i.e. the residue types aligned with each */
652 /* residue position in the current core model. These are stored in sai->sq. */
653
654 for(i=0; i<nsc; i++) {
655 mn=sai->al[i]-sli->no[i];
656 mx=sai->al[i]+sli->co[i];
657 j=cdf->sll.rfpt[i]-sli->no[i];
658 for(k=mn; k<=mx; k++) {
659 sai->sq[j]=qsq->sq[k];
660 j++;} }
661
662
663 /* Compute threaded sequence compostion. */
664 for(i=0; i<nsc; i++) spci(cdf,qsq,sli,sai,i,spc);
665 /* for(i=0;i<nrt;i++) printf("%d ",spc->rt[i]); printf("spc->rt\n"); */
666
667 /*
668 for(i=0;i<nsc;i++) {
669 for(j=0;j<nrt;j++) printf("%d ",spc->srt[i][j]);
670 printf("spc->srt[%d]\n",i); }
671 for(i=0;i<nlp;i++) {
672 for(j=0;j<nrt;j++) printf("%d ",spc->lrt[i][j]);
673 printf("spc->lrt[%d]\n",i); }
674 for(i=0;i<nrt;i++) printf("%d ",spc->rt[i]); printf("spc->rt\n");
675 for(i=0;i<nrt;i++) printf("%d ",spc->rto[i]); printf("spc->rto\n");
676 */
677
678 /* Compute contact counts. */
679 for(i=0; i<nsc; i++) spni(cpr,sli,i,spn);
680
681 /*
682 printf("spn->ndi:%d nsc:%d nrt%d\n",spn->ndi,spn->nsc,spn->nrt);
683 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
684 for(i=0; i<nsc; i++) printf("%d ",spn->nrr[k][j][i]);
685 printf("spn->nrr[%d][%d]\n",k,j); }
686 for(i=0; i<ndi; i++) printf("%d ",spn->srr[i]); printf("spn->srr\n");
687 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
688 for(i=0; i<nsc; i++) printf("%d ",spn->nrp[k][j][i]);
689 printf("spn->nrp[%d][%d]\n",k,j); }
690 for(i=0; i<ndi; i++) printf("%d ",spn->srp[i]); printf("spn->srp\n");
691 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
692 for(i=0; i<nrt; i++) printf("%d ",spn->nrf[k][j][i]);
693 printf("spn->nrf[%d][%d]\n",k,j); }
694 for(j=0; j<ndi; j++) {
695 for(i=0;i<nrt;i++) printf("%d ",spn->frf[j][i]);
696 printf("spn->frf[%d]\n",j); }
697 for(i=0; i<ndi; i++) printf("%d ",spn->srf[i]); printf("spn->srf\n");
698 */
699
700 /* Compute expected contact energy. */
701 cxei(spn,spc,pmf,sli,psm,cdf,cxe);
702
703
704 /* printf("cxe->nrt:%d ndi:%d\n",cxe->nrt,cxe->ndi);
705 for(i=0; i<nrt; i++) printf("%.4f ",cxe->rp[i]); printf("cxe->rp\n");
706 for(j=0; j<nrt; j++) {
707 for(i=0;i<nrt;i++) printf("%.4f ",cxe->rrp[j][i]);
708 printf("cxe->rrp[%d]\n",j); }
709 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rre[k]); printf("cxe->rre\n");
710 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rpe[k]); printf("cxe->rpe\n");
711 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rfe[i]); printf("cxe->rfe\n"); */
712
713 /* Construct contact pair lists for extent sampling, including energies */
714 cpll(cdf,pmf,qsq,cpr,sai,cpl);
715
716 /* Compute partial contact energy sums using the extent sampling lists. */
717 for(i=0; i<nsc; i++) spel(cpl,sai,sli,i,spe,cdf,psm,spc);
718
719 /* for(k=0; k<nsc; k++) printf("%d ",spe->gs[k]); printf("spe->gs\n");
720 for(k=0; k<nsc; k++) {
721 for(l=0; l<nsc; l++) printf("%d ",spe->gss[k][l]);
722 printf("spe->gss[%d]\n",k); } */
723
724 /* Compute thread energy for the initial model */
725 dgri(spe,spn,cxe,tdg,psm,spc);
726 /*printf("g:%d g0:%.2f dg:%.2f\n",tdg->g, tdg->g0, tdg->dg); */
727
728 /* Push the initial model onto the stack */
729 ttbi(sai,sli,tdg,ttb,nrs,gsp->als);
730
731 /* If indicated, record this value in the trajectory */
732 if(gsp->trg!=0) { trg[ngs]=tdg->dg; ngs++; }
733
734 /*
735 for(k=0;k<nsc;k++) printf("%d ",sai->al[k]); printf("sai->al\n");
736 for(k=0;k<nmt;k++) printf("%d ",sai->sq[k]); printf("sai->sq\n");
737
738 printf("ttb->n: %d mx: %d mn: %d\n",ttb->n,ttb->mx,ttb->mn);
739 for(k=0; k<ttb->n; k++) {
740 printf("k:%d tg:%.4f tf:%d ",k,ttb->tg[k],ttb->tf[k]);
741 printf("al: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->al[l][k]);
742 printf("no: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->no[l][k]);
743 printf("co: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->co[l][k]);
744 printf("nx:%d pr:%d\n",ttb->nx[k],ttb->pr[k]); }
745 */
746
747 /* Loop over the specified number of Gibb's sampling iterations. For each, */
748 /* perform the specified number of alignment and location sampling cycles, */
749 /* and push the resulting threads onto the best threads stack. */
750
751 /* Iterations may follow an annealing schedule as specified by temperature */
752 /* parameters associated with that iteration. The loop over iterations is */
753 /* therefore nested according to the number of temperature changes, and the */
754 /* number of iterations to conduct at each temperature. The temperature */
755 /* parameters for alignment and location cycles are specified separately. */
756
757 for(nts=0;nts<gsp->nts;nts++) {
758 /* printf("nts:%d\n",nts); */
759
760 for(nti=0;nti<gsp->nti[nts];nti++) {
761
762 /* Initialize storage for lowest energy found in this iteration */
763 tts->ib[nti]=BIGNEG;
764
765 if(gsp->nac[nts]>0) {
766
767 /* Construct contact pair lists for alignment sampling */
768 cpal(pmf,cpr,sli,cpa);
769
770 /* Initialize partial energy sums for alignment sampling */
771 for(i=0; i<nsc; i++) spea(cdf,cpa,sai,sli,i,1,spe,psm,spc);
772
773 /* Perform the indicated number of alignment cycles. In each cycle */
774 /* a new alignment for each core segments is chosen by the Gibbs */
775 /* sampling algorithm, conditioned on the current extent of core */
776 /* elements in the folding motif. */
777
778 for(nac=0;nac<gsp->nac[nts];nac++) {
779 /* printf("nac:%d\n",nac); */
780
781 /* Determine order for segment alignment. */
782 sgoi(gsp->iso,gsp->ito,pvl,sgo);
783
784 /* Loop over core segments. */
785 for(i=0; i<nsc; i++) { cs=sgo->si[i];
786
787 /* Find allowed alignment range for current segment */
788 if(salr(cdf,qsq,sli,sai,cs,&mn,&mx)==0) {
789 printf("failed salr nrs:%d, nts:%d nti:%d cs:%d\n",
790 nrs,nts,nti,cs);
791
792 for(i=0;i<nsc;i++) printf("%d ",sai->al[i]); printf("sai->al\n");
793 for(i=0;i<sai->nmt;i++) printf("%d ",sai->sq[i]); printf("sai->sq\n");
794
795 for(i=0;i<nsc;i++) printf("%d ",sli->no[i]); printf("sli->no\n");
796 for(i=0;i<nsc;i++) printf("%d ",sli->co[i]); printf("sli->co\n");
797 for(i=0;i<sli->nlp;i++) printf("%d ",sli->lp[i]); printf("sli->lp\n");
798 for(i=0;i<sli->nmt;i++) printf("%d ",sli->cr[i]); printf("sli->cr\n");
799
800 return(0); }
801
802
803 /* Loop over values of the aligment variable. */
804 pvl->n=mx-mn+1;
805 /* printf("cs:%d mn:%d mx:%d\n",cs,mn,mx); */
806 for(j=0; j<=mx-mn; j++) { al=mn+j;
807
808 /* Assign the alignment. */
809 salu(cdf,qsq,sli,cs,al,sai);
810 /* printf("assigned cs:%d al:%d\n",cs,al); */
811
812 /* Update contact energies. */
813 spea(cdf,cpa,sai,sli,cs,1,spe,psm,spc);
814 /* printf("updated energies cs:%d\n",cs); */
815
816 /* Update composition of threaded sequence */
817 spcd=spci(cdf,qsq,sli,sai,cs,spc);
818 /* printf("spcd:%d\n",spcd); */
819
820 /* If necessary update expected energy. */
821 if(spcd!=0) cxei(spn,spc,pmf,sli,psm,cdf,cxe);
822
823 /* Assign thread energy for this aligment. */
824 pvl->e[j]=dgri(spe,spn,cxe,tdg,psm,spc);
825 }
826
827
828 /* Calculate probabilities and choose an alignment. */
829 al=mn+algs(pvl,gsp->tma[nts]);
830
831 /* Assign the chosen alignment for this segment. */
832 salu(cdf,qsq,sli,cs,al,sai);
833 /* printf("choose by sampling cs:%d al:%d\n",cs,al); */
834
835 /* Update contact energies. */
836 spea(cdf,cpa,sai,sli,cs,1,spe,psm,spc);
837
838 /* Update composition. */
839 spcd=spci(cdf,qsq,sli,sai,cs,spc);
840
841 /* If necessary update expected energies. */
842 if(spcd) cxei(spn,spc,pmf,sli,psm,cdf,cxe);
843
844 /* Assign thread energy */
845 dgri(spe,spn,cxe,tdg,psm,spc);
846 }
847
848
849 /* Push the sampled thread onto the stack. */
850 ttbi(sai,sli,tdg,ttb,nrs,gsp->als);
851
852 /* Record the lowest energy found in this iteration */
853 if(tdg->dg>tts->ib[nti]) tts->ib[nti]=tdg->dg;
854
855 /* If indicated, record this value in the trajectory */
856 if(gsp->trg!=0) { trg[ngs]=tts->ib[nti]; ngs++; }
857
858 } }
859
860 /* printf("After alignment cycle\n");
861
862 for(k=0;k<nsc;k++) printf("%d ",sai->al[k]); printf("sai->al\n");
863 for(k=0;k<nmt;k++) printf("%d ",sai->sq[k]); printf("sai->sq\n");
864
865 for(k=0; k<nsc; k++) printf("%d ",spe->gs[k]); printf("spe->gs\n");
866 for(k=0; k<nsc; k++) {
867 for(l=0; l<nsc; l++) printf("%d ",spe->gss[k][l]);
868 printf("spe->gss[%d]\n",k); }
869 printf("nts:%d nti:%d\n",nts,nti);
870 printf("ps:%.2f ms:%.2f dg:%.2f\n",tdg->ps, tdg->ms, tdg->dg);
871
872 printf("ttb->n: %d mx: %d mn: %d\n",ttb->n,ttb->mx,ttb->mn);
873 for(k=0; k<ttb->n; k++) {
874 printf("k:%d tg:%.4f tf:%d ",k,ttb->tg[k],ttb->tf[k]);
875 printf("al: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->al[l][k]);
876 printf("no: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->no[l][k]);
877 printf("co: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->co[l][k]);
878 printf("nx:%d pr:%d\n",ttb->nx[k],ttb->pr[k]); }
879
880 printf("spn->ndi:%d nsc:%d nrt%d\n",spn->ndi,spn->nsc,spn->nrt);
881 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
882 for(i=0; i<nsc; i++) printf("%d ",spn->nrr[k][j][i]);
883 printf("spn->nrr[%d][%d]\n",k,j); }
884 for(i=0; i<ndi; i++) printf("%d ",spn->srr[i]); printf("spn->srr\n");
885 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
886 for(i=0; i<nsc; i++) printf("%d ",spn->nrp[k][j][i]);
887 printf("spn->nrp[%d][%d]\n",k,j); }
888 for(i=0; i<ndi; i++) printf("%d ",spn->srp[i]); printf("spn->srp\n");
889 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
890 for(i=0; i<nrt; i++) printf("%d ",spn->nrf[k][j][i]);
891 printf("spn->nrf[%d][%d]\n",k,j); }
892 for(j=0; j<ndi; j++) {
893 for(i=0;i<nrt;i++) printf("%d ",spn->frf[j][i]);
894 printf("spn->frf[%d]\n",j); }
895 for(i=0; i<ndi; i++) printf("%d ",spn->srf[i]); printf("spn->srf\n");
896
897 printf("spc->nsc:%d nrt:%d nlp:%d\n",spc->nsc,spc->nrt,spc->nlp);
898 for(j=0; j<nsc; j++) {
899 for(i=0;i<nrt;i++) printf("%d ",spc->srt[j][i]);
900 printf("spc->srt[%d]\n",j); }
901 for(j=0; j<nlp; j++) {
902 for(i=0;i<nrt;i++) printf("%d ",spc->lrt[j][i]);
903 printf("spc->lrt[%d]\n",j); }
904 for(i=0; i<nrt; i++) printf("%d ",spc->rt[i]); printf("spc->rt\n");
905 for(i=0; i<nrt; i++) printf("%d ",spc->rto[i]); printf("spc->rto\n");
906
907 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rre[k]); printf("cxe->rre\n");
908 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rpe[k]); printf("cxe->rpe\n");
909 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rfe[i]); printf("cxe->rfe\n");
910 */
911
912 if(gsp->nlc[nts]>0) {
913
914 /* Construct contact pair lists for segment extent sampling. */
915 cpll(cdf,pmf,qsq,cpr,sai,cpl);
916
917 /* Initialize partial energy sums for segment extent sampling. */
918 for(i=0; i<nsc; i++) spel(cpl,sai,sli,i,spe,cdf,psm,spc);
919
920 /* Perform the indicated number of segment extent sampling cycles. */
921 /* For each cycle, a new extent for the n- and c-termini of each */
922 /* core element is chosen by the Gibbs sampling algorithm, */
923 /* conditioned on the current aligment of the query sequence. */
924
925 for(nlc=0;nlc<gsp->nlc[nts];nlc++) {
926
927 /* Determine the order of segment location. */
928 sgoi(gsp->iso,gsp->ito,pvl,sgo);
929
930 /* Loop over core segments */
931 for(i=0; i<nsc; i++) { cs=sgo->si[i];
932
933 /* Loop over termini */
934 for(j=0;j<=1;j++) {
935
936 switch(j) { case 0: { ct=sgo->to[cs];
937 break; }
938 case 1: { ct=(sgo->to[cs]==0) ? 1:0;
939 break; }
940 }
941
942 /* Find allowed extent range for current terminus. */
943 if(slor(mtf,cdf,qsq,sli,sai,cs,ct,&mn,&mx)==0) {
944 printf("failed slor nrs:%d nts:%d nti:%d cs:%d ct:%d\n",
945 nrs,nts,nti,cs,ct);
946
947 for(k=0;k<nsc;k++) printf("%d ",sai->al[k]); printf("sai->al\n");
948 for(k=0;k<nmt;k++) printf("%d ",sai->sq[k]); printf("sai->sq\n");
949
950 for(i=0;i<nsc;i++) printf("%d ",sli->no[i]); printf("sli->no\n");
951 for(i=0;i<nsc;i++) printf("%d ",sli->co[i]); printf("sli->co\n");
952 for(i=0;i<nlp;i++) printf("%d ",sli->lp[i]); printf("sli->lp\n");
953 for(i=0;i<nmt;i++) printf("%d ",sli->cr[i]); printf("sli->cr\n");
954
955 return(0); }
956 /* printf("cs:%d ct:%d mn:%d mx:%d\n",cs,ct,mn,mx); */
957
958
959 /* Loop over possible extent values. */
960 pvl->n=mx-mn+1;
961 for(k=0; k<=mx-mn; k++) { of=mn+k;
962
963 /* Assign the extent. */
964 slou(mtf,cdf,cs,ct,of,sli,sai,qsq);
965 /* printf("assigned cs:%d ct:%d of:%d\n",
966 cs,ct,of); */
967
968 /* Update contact energies. */
969 spel(cpl,sai,sli,cs,spe,cdf,psm,spc);
970
971 /* Update composition. */
972 spci(cdf,qsq,sli,sai,cs,spc);
973
974 /* Update contact counts. */
975 if (PSSM_Weight < 0.9999999) {
976 spni(cpl,sli,cs,spn);
977 }
978
979 /* Update expected contact energy. */
980 cxei(spn,spc,pmf,sli,psm,cdf,cxe);
981
982 /* Assign thread energy for this extent. */
983 pvl->e[k]=dgri(spe,spn,cxe,tdg,psm,spc);
984 }
985
986 /* Calculate probabilities and choose an extent. */
987 of=mn+algs(pvl,gsp->tml[nts]);
988
989 /* Assign the chosen extent for this terminus. */
990 slou(mtf,cdf,cs,ct,of,sli,sai,qsq);
991 /* printf("chose by sampling cs:%d ct:%d of:%d\n",
992 cs,ct,of); */
993
994 /* Update contact energies. */
995 spel(cpl,sai,sli,cs,spe,cdf,psm,spc);
996
997 /* Update composition. */
998 spci(cdf,qsq,sli,sai,cs,spc);
999
1000 /* Update contact counts. */
1001 if (PSSM_Weight < 0.9999999) {
1002 spni(cpl,sli,cs,spn);
1003 }
1004
1005 /* Update expected contact energy. */
1006 cxei(spn,spc,pmf,sli,psm,cdf,cxe);
1007
1008 /* Update thread energy */
1009 dgri(spe,spn,cxe,tdg,psm,spc); } }
1010
1011 /* Push the sampled thread onto the stack. */
1012 ttbi(sai,sli,tdg,ttb,nrs,gsp->als);
1013
1014 /* Record the lowest energy found in this iteration */
1015 if(tdg->dg>tts->ib[nti]) {tts->ib[nti]=tdg->dg;
1016 /*hh=tdg->ms;*/}
1017
1018 /*printf("%d,%d,%d\n",nrs,nts,nti);
1019 printf("%f,%f\n",tts->ib[nti],hh);*/
1020
1021 /* printf("tts->ib[%d]=%.2f\n",nti,tts->ib[nti]); */
1022
1023 /* If indicated, record this value in the trajectory */
1024 if(gsp->trg!=0) { trg[ngs]=tts->ib[nti]; ngs++; }
1025
1026 } }
1027
1028 /*
1029 for(k=0;k<nsc;k++) printf("%d ",sai->al[k]); printf("sai->al\n");
1030 for(k=0;k<nmt;k++) printf("%d ",sai->sq[k]); printf("sai->sq\n");
1031
1032 for(k=0; k<nsc; k++) printf("%d ",spe->gs[k]); printf("spe->gs\n");
1033 for(k=0; k<nsc; k++) {
1034 for(l=0; l<nsc; l++) printf("%d ",spe->gss[k][l]);
1035 printf("spe->gss[%d]\n",k); }
1036 printf("nts:%d nti:%d\n",nts,nti);
1037 printf("ps:%.2f ms:%.2f dg:%.2f\n",tdg->ps, tdg->ms, tdg->dg);
1038
1039 printf("ttb->n: %d mx: %d mn: %d\n",ttb->n,ttb->mx,ttb->mn);
1040 for(k=0; k<ttb->n; k++) {
1041 printf("k:%d tg:%.4f tf:%d ",k,ttb->tg[k],ttb->tf[k]);
1042 printf("al: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->al[l][k]);
1043 printf("no: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->no[l][k]);
1044 printf("co: "); for(l=0; l<sli->nsc; l++) printf("%d ",ttb->co[l][k]);
1045 printf("nx:%d pr:%d\n",ttb->nx[k],ttb->pr[k]); }
1046
1047 printf("spn->ndi:%d nsc:%d nrt%d\n",spn->ndi,spn->nsc,spn->nrt);
1048 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
1049 for(i=0; i<nsc; i++) printf("%d ",spn->nrr[k][j][i]);
1050 printf("spn->nrr[%d][%d]\n",k,j); }
1051 for(i=0; i<ndi; i++) printf("%d ",spn->srr[i]); printf("spn->srr\n");
1052 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
1053 for(i=0; i<nsc; i++) printf("%d ",spn->nrp[k][j][i]);
1054 printf("spn->nrp[%d][%d]\n",k,j); }
1055 for(i=0; i<ndi; i++) printf("%d ",spn->srp[i]); printf("spn->srp\n");
1056 for(k=0; k<ndi; k++) for(j=0; j<nsc; j++) {
1057 for(i=0; i<nrt; i++) printf("%d ",spn->nrf[k][j][i]);
1058 printf("spn->nrf[%d][%d]\n",k,j); }
1059 for(j=0; j<ndi; j++) {
1060 for(i=0;i<nrt;i++) printf("%d ",spn->frf[j][i]);
1061 printf("spn->frf[%d]\n",j); }
1062 for(i=0; i<ndi; i++) printf("%d ",spn->srf[i]); printf("spn->srf\n");
1063
1064 printf("spc->nsc:%d nrt:%d nlp:%d\n",spc->nsc,spc->nrt,spc->nlp);
1065 for(j=0; j<nsc; j++) {
1066 for(i=0;i<nrt;i++) printf("%d ",spc->srt[j][i]);
1067 printf("spc->srt[%d]\n",j); }
1068 for(j=0; j<nlp; j++) {
1069 for(i=0;i<nrt;i++) printf("%d ",spc->lrt[j][i]);
1070 printf("spc->lrt[%d]\n",j); }
1071 for(i=0; i<nrt; i++) printf("%d ",spc->rt[i]); printf("spc->rt\n");
1072 for(i=0; i<nrt; i++) printf("%d ",spc->rto[i]); printf("spc->rto\n");
1073
1074 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rre[k]); printf("cxe->rre\n");
1075 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rpe[k]); printf("cxe->rpe\n");
1076 for(k=0; k<ndi; k++) printf("%.4f ",cxe->rfe[i]); printf("cxe->rfe\n");
1077 */
1078
1079 /* Test for local minimum within this iteration. Set flag if found. */
1080 tts->lm=0;
1081 /* printf("tts->ib[%d]=%.2f mx=%.2f\n",
1082 nti,tts->ib[nti],ttb->tg[ttb->mx]); */
1083 if(gsp->lms[nts]>0 && (nti+1)>=gsp->lms[nts]) {
1084 tts->gb=BIGNEG; for(i=nti-gsp->lmw[nts]+1;i<=nti;i++)
1085 if(tts->gb<tts->ib[i]) tts->gb=tts->ib[i];
1086 if(tts->gb<(((float)gsp->lmf[nts])/100.)*ttb->tg[ttb->mx]) {
1087 /* printf("gb:%.2f mx:%.2f\n", tts->gb, ttb->tg[ttb->mx]); */
1088 tts->lm=1; break; }}
1089
1090 /* End of loop over iterations per temperature step */
1091 }
1092
1093
1094 /* If local minimum flag is set, stop sampling for this random start */
1095 if(tts->lm==1) break;
1096
1097 /* End of loop over temperature steps */
1098 }
1099
1100
1101 /* Test for convergence based on recurrence of the same alignments */
1102
1103 /* bwfi(ttb,gsp,tts); printf("nrs:%d ts:%d tf:%d\n",nrs,tts->ts,tts->tf); */
1104 if(nrs>=gsp->crs) {
1105 bwfi(ttb,gsp,tts);
1106
1107 /* Stop if the number of starts and frequency crosses a threshold */
1108 if(tts->tf>=gsp->cfm && tts->ts>=gsp->csm) break; }
1109
1110
1111 /* End of loop over random starts */
1112 }
1113
1114 /* put results in order, from highest score to lowest */
1115 OrderThdTbl(ttb);
1116
1117 /* Calculate z_scores */
1118 switch(zscs) {
1119 case 0:
1120 /* don't calculate any z-scores */
1121 break;
1122 case 1:
1123 /* calculate z-score for top alignment */
1124 zsc(ttb,psm,qsq,cpr,cdf,pmf,spe,sai,pvl,ScalingFactor);
1125 break;
1126 default:
1127 /* calculate z-scores for all top alignments */
1128 tmp = ttb->mx;
1129 for (i=0; i<ttb->n; i++) {
1130 ttb->mx = i;
1131 zsc(ttb,psm,qsq,cpr,cdf,pmf,spe,sai,pvl,ScalingFactor);
1132 }
1133 ttb->mx = tmp;
1134 break;
1135 }
1136
1137 /* scale energies back down */
1138 ScaleThdTbl(ttb, 1.0/ScalingFactor);
1139
1140 /*
1141
1142 for(i=0;i<ttb->n;i++) {
1143 printf("i:%d\n",i);
1144 printf("ttb->tf:%d",ttb->tf[i]);
1145 printf("ttb->ts:%d",ttb->ts[i]);
1146 printf("ttb->ls:%d\n",ttb->ls[i]);
1147
1148 for(j=0;j<nsc;j++){
1149 printf("ttb->al:%d",ttb->al[j][i]);
1150 printf("ttb->no:%d",ttb->no[j][i]);
1151 printf("ttb->co:%dn",ttb->co[j][i]);}
1152
1153 }
1154
1155 */
1156
1157 /* Free all allocated storage */
1158
1159 free(sli->no); /* In sli, type Cur_Loc */
1160 free(sli->co);
1161 free(sli->lp);
1162 free(sli->cr);
1163 free(sli);
1164 /* printf("freed sli\n"); */
1165
1166 free(sai->al); /* In sai, type Cur_Aln */
1167 free(sai->sq);
1168 free(sai->cf);
1169 free(sai);
1170 /* printf("freed sai\n"); */
1171
1172 free(spe->gs); /* In spe, type Seg_Gsm */
1173 free(spe->ms);
1174 free(spe->cs);
1175 free(spe->ls);
1176 free(spe->s0);
1177 for(i=0;i<nsc;i++) free(spe->gss[i]); free(spe->gss);
1178 free(spe);
1179 /* printf("freed spe\n"); */
1180
1181 free(spn->srr); /* In spn, type Seg_Nsm */
1182 for(i=0; i<ndi; i++)
1183 for (j=0; j<nsc; j++) free(spn->nrr[i][j]);
1184 for(i=0; i<ndi; i++) free(spn->nrr[i]); free(spn->nrr);
1185 free(spn->srp);
1186 for(i=0; i<ndi; i++)
1187 for (j=0; j<nsc; j++) free(spn->nrp[i][j]);
1188 for(i=0; i<ndi; i++) free(spn->nrp[i]); free(spn->nrp);
1189 free(spn->srf);
1190 for(i=0; i<ndi; i++)
1191 for (j=0; j<nsc; j++) free(spn->nrf[i][j]);
1192 for(i=0; i<ndi; i++) free(spn->nrf[i]); free(spn->nrf);
1193 for(i=0; i<ndi; i++) free(spn->frf[i]); free(spn->frf);
1194 free(spn);
1195 /* printf("freed spn\n"); */
1196
1197
1198 for(i=0;i<nsc;i++) { /* In cpr, type Cxl_Los */
1199 cr=cpr[i];
1200 free(cr->rr.r1);
1201 free(cr->rr.r2);
1202 free(cr->rr.d);
1203 free(cr->rp.r1);
1204 free(cr->rp.p2);
1205 free(cr->rp.d);
1206 free(cr->rf.r1);
1207 free(cr->rf.t2);
1208 free(cr->rf.d); }
1209 for(i=0;i<nsc;i++) free(cpr[i]);
1210 free(cpr);
1211 /* printf("freed cpr\n"); */
1212
1213
1214 for(i=0;i<nsc;i++) { /* In cpl, type Cxl_Los */
1215 cl=cpl[i];
1216 free(cl->rr.r1);
1217 free(cl->rr.r2);
1218 free(cl->rr.d);
1219 free(cl->rr.e);
1220 free(cl->rp.r1);
1221 free(cl->rp.p2);
1222 free(cl->rp.d);
1223 free(cl->rp.e);
1224 free(cl->rf.r1);
1225 free(cl->rf.t2);
1226 free(cl->rf.d);
1227 free(cl->rf.e); }
1228 for(i=0;i<nsc;i++) free(cpl[i]);
1229 free(cpl);
1230 /* printf("freed cpl\n"); */
1231
1232
1233 for(i=0;i<nsc;i++) { /* In cpa, type Cxl_Als */
1234 ca=cpa[i];
1235 free(ca->rr.r1);
1236 free(ca->rr.r2);
1237 free(ca->rr.e);
1238 free(ca->r.r);
1239 n=cdf->sll.nomx[i]+cdf->sll.comx[i]+1;
1240 for(j=0;j<n;j++) free(ca->r.e[j]);
1241 free(ca->r.e);}
1242 for(i=0;i<nsc;i++) free(cpa[i]);
1243 free(cpa);
1244 /* printf("freed cpa\n"); */
1245
1246
1247 free(spc->rt); /* In spc, type Seg_Cmp */
1248 free(spc->rto);
1249 for(i=0; i<nsc; i++) free(spc->srt[i]);
1250 free(spc->srt);
1251 for(i=0; i<nlp; i++) free(spc->lrt[i]);
1252 free(spc->lrt);
1253 free(spc);
1254 /* printf("freed spc\n"); */
1255
1256 free(cxe->rp); /* In cxe, type Thd_Cxe */
1257 for(i=0; i<nrt; i++) free(cxe->rrp[i]);
1258 free(cxe->rrp);
1259 free(cxe->rre);
1260 free(cxe->rpe);
1261 free(cxe->rfe);
1262 free(cxe);
1263 /* printf("freed cxe\n"); */
1264
1265 free(sgo->si); /* In sgo; type Seg_Ord */
1266 free(sgo->so);
1267 free(sgo->to);
1268 free(sgo);
1269 /* printf("freed sgo\n"); */
1270
1271 free(pvl->p); /* In pvl; type Rnd_Smp */
1272 free(pvl->e);
1273 free(pvl->sq);
1274 free(pvl->aqi);
1275 free(pvl->r);
1276 free(pvl->o);
1277 free(pvl->lsg);
1278 free(pvl);
1279 /* printf("freed pvl\n"); */
1280
1281 free(tts->bw); /* In tts; type Thd_Tst */
1282 free(tts->ib);
1283 free(tts);
1284 /* printf("freed tts\n"); */
1285
1286 free(tdg);
1287
1288 return(1); /* successful completion */
1289 }
1290
1291
NewCorDef(int NumBlocks)1292 Cor_Def* NewCorDef(int NumBlocks) {
1293 /*----------------------------------------------------------*/
1294 /* allocate space for a new Cor_Def. */
1295 /*----------------------------------------------------------*/
1296 Cor_Def* cdfp;
1297
1298 cdfp = (Cor_Def*) calloc(1,sizeof(Cor_Def));
1299 cdfp->sll.rfpt = (int*) calloc(1,sizeof(int) * NumBlocks);
1300 cdfp->sll.nomn = (int*) calloc(1,sizeof(int) * NumBlocks);
1301 cdfp->sll.nomx = (int*) calloc(1,sizeof(int) * NumBlocks);
1302 cdfp->sll.comn = (int*) calloc(1,sizeof(int) * NumBlocks);
1303 cdfp->sll.comx = (int*) calloc(1,sizeof(int) * NumBlocks);
1304 cdfp->sll.n = NumBlocks;
1305 cdfp->lll.llmn = (int*) calloc(1,sizeof(int) * (NumBlocks+1));
1306 cdfp->lll.llmx = (int*) calloc(1,sizeof(int) * (NumBlocks+1));
1307 cdfp->lll.lrfs = (int*) calloc(1,sizeof(int) * (NumBlocks+1));
1308 cdfp->lll.n = NumBlocks+1;
1309 return(cdfp);
1310 }
1311
1312
FreeCorDef(Cor_Def * cdfp)1313 Cor_Def* FreeCorDef(Cor_Def* cdfp) {
1314 /*----------------------------------------------------------*/
1315 /* free Cor_Def. */
1316 /*----------------------------------------------------------*/
1317 free(cdfp->sll.rfpt);
1318 free(cdfp->sll.nomn);
1319 free(cdfp->sll.nomx);
1320 free(cdfp->sll.comn);
1321 free(cdfp->sll.comx);
1322 free(cdfp->lll.llmn);
1323 free(cdfp->lll.llmx);
1324 free(cdfp->lll.lrfs);
1325 free(cdfp);
1326 return NULL;
1327 }
1328
1329
PrintCorDef(Cor_Def * cdf,FILE * pFile)1330 void PrintCorDef(Cor_Def* cdf, FILE* pFile) {
1331 /*----------------------------------------------------------*/
1332 /* for debugging. Print CorDef. */
1333 /*----------------------------------------------------------*/
1334 int i;
1335 static Boolean FirstPass = TRUE;
1336
1337 if (FirstPass == FALSE) return;
1338 FirstPass = FALSE;
1339
1340 fprintf(pFile, "Core Definition:\n");
1341 fprintf(pFile, "number of core segments: %4d\n", cdf->sll.n);
1342 fprintf(pFile, " nomx nomn rfpt comn comx\n");
1343 for (i=0; i<cdf->sll.n; i++) {
1344 fprintf(pFile, " %6d %6d %6d %6d %6d\n",
1345 cdf->sll.nomx[i], cdf->sll.nomn[i], cdf->sll.rfpt[i], cdf->sll.comn[i], cdf->sll.comx[i]);
1346 }
1347 fprintf(pFile, "number of loops (one more than core segs): %4d\n", cdf->lll.n);
1348 fprintf(pFile, " llmn llmx lrfs\n");
1349 for (i=0; i<cdf->lll.n; i++) {
1350 fprintf(pFile, " %6d %6d %6d\n", cdf->lll.llmn[i], cdf->lll.llmx[i], cdf->lll.lrfs[i]);
1351 }
1352 fprintf(pFile, "number of fixed segments: %4d\n", cdf->fll.n);
1353 fprintf(pFile, " nt ct sq\n");
1354 for (i=0; i<cdf->fll.n; i++) {
1355 fprintf(pFile, " %6d %6d %6d\n", cdf->fll.nt[i], cdf->fll.ct[i], cdf->fll.sq[i]);
1356 }
1357 fprintf(pFile, "\n");
1358 }
1359
1360
NewSeqMtf(int NumResidues,int AlphabetSize)1361 Seq_Mtf* NewSeqMtf(int NumResidues, int AlphabetSize) {
1362 /*----------------------------------------------------------*/
1363 /* allocate space for a new pssm. */
1364 /*----------------------------------------------------------*/
1365 int i;
1366 Seq_Mtf* pssm;
1367
1368 pssm = (Seq_Mtf*) calloc(1,sizeof(Seq_Mtf));
1369 pssm->n = NumResidues;
1370 pssm->AlphabetSize = AlphabetSize;
1371 pssm->ww = (int**) calloc(1,NumResidues * sizeof(int*));
1372 pssm->freqs = (int**) calloc(1,NumResidues * sizeof(int*));
1373 for (i=0; i<NumResidues; i++) {
1374 pssm->ww[i] = (int*) calloc(1,AlphabetSize * sizeof(int));
1375 pssm->freqs[i] = (int*) calloc(1,AlphabetSize * sizeof(int));
1376 }
1377 return(pssm);
1378 }
1379
1380
FreeSeqMtf(Seq_Mtf * pssm)1381 Seq_Mtf* FreeSeqMtf(Seq_Mtf* pssm) {
1382 /*----------------------------------------------------------*/
1383 /* free pssm. */
1384 /*----------------------------------------------------------*/
1385 int i;
1386
1387 for (i=0; i<pssm->n; i++) {
1388 free(pssm->ww[i]);
1389 free(pssm->freqs[i]);
1390 }
1391 free(pssm->ww);
1392 free(pssm->freqs);
1393 free(pssm);
1394 return NULL;
1395 }
1396
1397
PrintSeqMtf(Seq_Mtf * psm,FILE * pFile)1398 void PrintSeqMtf(Seq_Mtf* psm, FILE* pFile) {
1399 /*----------------------------------------------------------*/
1400 /* print the pssm. */
1401 /*----------------------------------------------------------*/
1402 char OutputOrder[] = "ARNDCQEGHILKMFPSTWYV";
1403 int i, j;
1404
1405 fprintf(pFile, "PSSM:\n");
1406 fprintf(pFile, "Number of residues: %4d\n", psm->n);
1407 fprintf(pFile, "Weights:\n");
1408 for (i=0; i<strlen(OutputOrder); i++) {
1409 fprintf(pFile, "%8c ", OutputOrder[i]);
1410 }
1411 fprintf(pFile, "\n");
1412 for (i=0; i<psm->n; i++) {
1413 for (j=0; j<strlen(OutputOrder); j++) {
1414 fprintf(pFile, "%8d ", psm->ww[i][j]);
1415 }
1416 fprintf(pFile, "\n");
1417 }
1418 fprintf(pFile, "Frequencies:\n");
1419 for (i=0; i<strlen(OutputOrder); i++) {
1420 fprintf(pFile, "%8c ", OutputOrder[i]);
1421 }
1422 fprintf(pFile, "\n");
1423 for (i=0; i<psm->n; i++) {
1424 for (j=0; j<strlen(OutputOrder); j++) {
1425 fprintf(pFile, "%8d ", psm->freqs[i][j]);
1426 }
1427 fprintf(pFile, "\n");
1428 }
1429 fprintf(pFile, "\n");
1430 }
1431
1432
NewQrySeq(int NumResidues,int NumBlocks)1433 Qry_Seq* NewQrySeq(int NumResidues, int NumBlocks) {
1434 /*----------------------------------------------------------*/
1435 /* allocate space for a new query sequence. */
1436 /*----------------------------------------------------------*/
1437 Qry_Seq* qsqp;
1438 int i;
1439
1440 qsqp = (Qry_Seq*) calloc(1,sizeof(Qry_Seq));
1441 qsqp->n = NumResidues;
1442 qsqp->sq = (int*) calloc(1,sizeof(int) * NumResidues);
1443 qsqp->sac.n = NumBlocks;
1444 qsqp->sac.mn = (int*) calloc(1,sizeof(int) * NumBlocks);
1445 qsqp->sac.mx = (int*) calloc(1,sizeof(int) * NumBlocks);
1446 /* for no constraints, set mn and mx to -1 for each block */
1447 /* make no constraints the default */
1448 for (i=0; i<NumBlocks; i++) {
1449 qsqp->sac.mn[i] = -1;
1450 qsqp->sac.mx[i] = -1;
1451 }
1452 return(qsqp);
1453 }
1454
1455
FreeQrySeq(Qry_Seq * qsqp)1456 Qry_Seq* FreeQrySeq(Qry_Seq* qsqp) {
1457 /*----------------------------------------------------------*/
1458 /* free query sequence. */
1459 /*----------------------------------------------------------*/
1460 free(qsqp->sac.mn);
1461 free(qsqp->sac.mx);
1462 free(qsqp->sq);
1463 free(qsqp);
1464 return NULL;
1465 }
1466
1467
PrintQrySeq(Qry_Seq * qsqp,FILE * pFile)1468 void PrintQrySeq(Qry_Seq* qsqp, FILE* pFile) {
1469 /*----------------------------------------------------------*/
1470 /* for debugging. print query sequence. */
1471 /*----------------------------------------------------------*/
1472 int i, j, Count;
1473 char ResLetters[32];
1474
1475 /* lookup to give one-letter names of residues */
1476 strcpy(ResLetters, "ARNDCQEGHILKMFPSTWYV");
1477
1478 fprintf(pFile, "Query Sequence:\n");
1479 fprintf(pFile, "number of residues: %4d\n", qsqp->n);
1480 Count = 0;
1481 for (i=0; i<=qsqp->n/30; i++) {
1482 for (j=0; j<30; j++) {
1483 fprintf(pFile, "%c ", ResLetters[qsqp->sq[i*30 + j]]);
1484 Count++;
1485 if (Count == qsqp->n) {
1486 fprintf(pFile, "\n");
1487 goto Here;
1488 }
1489 }
1490 fprintf(pFile, "\n");
1491 }
1492
1493 Here:
1494 fprintf(pFile, "number of constraints: %4d\n", qsqp->sac.n);
1495 fprintf(pFile, " mn mx\n");
1496 for (i=0; i<qsqp->sac.n; i++) {
1497 fprintf(pFile, " %6d %6d\n", qsqp->sac.mn[i], qsqp->sac.mx[i]);
1498 }
1499 fprintf(pFile, "\n");
1500 }
1501
1502
NewRcxPtl(int NumResTypes,int NumDistances,int PeptideIndex)1503 Rcx_Ptl* NewRcxPtl(int NumResTypes, int NumDistances, int PeptideIndex) {
1504 /*----------------------------------------------------------*/
1505 /* allocate space for the contact potential */
1506 /*----------------------------------------------------------*/
1507 Rcx_Ptl* pmf;
1508 int i, j;
1509
1510 pmf = (Rcx_Ptl*) calloc(1,sizeof(Rcx_Ptl));
1511 pmf->rre = (int***) calloc(1,NumDistances * sizeof(int**));
1512 pmf->rrt = (int***) calloc(1,NumDistances * sizeof(int**));
1513 pmf->re = (int**) calloc(1,NumDistances * sizeof(int*));
1514 for (i=0; i<NumDistances; i++) {
1515 pmf->rre[i] = (int**) calloc(1,NumResTypes * sizeof(int*));
1516 pmf->rrt[i] = (int**) calloc(1,NumResTypes * sizeof(int*));
1517 for (j=0; j<NumResTypes; j++) {
1518 pmf->rre[i][j] = (int*) calloc(1,NumResTypes * sizeof(int));
1519 pmf->rrt[i][j] = (int*) calloc(1,NumResTypes * sizeof(int));
1520 }
1521 pmf->re[i] = (int*) calloc(1,NumResTypes * sizeof(int));
1522 }
1523 pmf->nrt = NumResTypes;
1524 pmf->ndi = NumDistances;
1525 pmf->ppi = PeptideIndex;
1526 return(pmf);
1527 }
1528
1529
FreeRcxPtl(Rcx_Ptl * pmf)1530 Rcx_Ptl* FreeRcxPtl(Rcx_Ptl* pmf) {
1531 /*----------------------------------------------------------*/
1532 /* free contact potential. */
1533 /*----------------------------------------------------------*/
1534 int i, j;
1535
1536 for (i=0; i<pmf->ndi; i++) {
1537 for (j=0; j<pmf->nrt; j++) {
1538 free(pmf->rre[i][j]);
1539 free(pmf->rrt[i][j]);
1540 }
1541 free(pmf->rre[i]);
1542 free(pmf->rrt[i]);
1543 free(pmf->re[i]);
1544 }
1545 free(pmf->rre);
1546 free(pmf->rrt);
1547 free(pmf->re);
1548 free(pmf);
1549 return NULL;
1550 }
1551
1552
NewGibScd(int NumTempSteps)1553 Gib_Scd* NewGibScd(int NumTempSteps) {
1554 /*----------------------------------------------------------*/
1555 /* allocate space for annealing schedule. */
1556 /*----------------------------------------------------------*/
1557 Gib_Scd* gsp;
1558
1559 gsp = (Gib_Scd*) calloc(1,sizeof(Gib_Scd));
1560 gsp->nts = NumTempSteps;
1561 gsp->nti = (int*) calloc(1,NumTempSteps * sizeof(int));
1562 gsp->nac = (int*) calloc(1,NumTempSteps * sizeof(int));
1563 gsp->nlc = (int*) calloc(1,NumTempSteps * sizeof(int));
1564 gsp->tma = (int*) calloc(1,NumTempSteps * sizeof(int));
1565 gsp->tml = (int*) calloc(1,NumTempSteps * sizeof(int));
1566 gsp->lms = (int*) calloc(1,NumTempSteps * sizeof(int));
1567 gsp->lmw = (int*) calloc(1,NumTempSteps * sizeof(int));
1568 gsp->lmf = (int*) calloc(1,NumTempSteps * sizeof(int));
1569 return(gsp);
1570 }
1571
1572
FreeGibScd(Gib_Scd * gsp)1573 Gib_Scd* FreeGibScd(Gib_Scd* gsp) {
1574 /*----------------------------------------------------------*/
1575 /* free the annealing schedule. */
1576 /*----------------------------------------------------------*/
1577 free(gsp->nti);
1578 free(gsp->nac);
1579 free(gsp->nlc);
1580 free(gsp->tma);
1581 free(gsp->tml);
1582 free(gsp->lms);
1583 free(gsp->lmw);
1584 free(gsp->lmf);
1585 free(gsp);
1586 return NULL;
1587 }
1588
1589
NewFldMtf(int NumResidues,int NumResResContacts,int NumResPepContacts)1590 Fld_Mtf* NewFldMtf(int NumResidues, int NumResResContacts, int NumResPepContacts) {
1591 /*-------------------------------------------------------------*/
1592 /* allocate space for contact lists and minimum loop lengths */
1593 /*-------------------------------------------------------------*/
1594 int i;
1595 Fld_Mtf* mtf;
1596
1597 mtf = (Fld_Mtf*) calloc(1,sizeof(Fld_Mtf));
1598 mtf->n = NumResidues;
1599 mtf->rrc.n = NumResResContacts;
1600 mtf->rrc.r1 = (int*) calloc(1,NumResResContacts * sizeof(int));
1601 mtf->rrc.r2 = (int*) calloc(1,NumResResContacts * sizeof(int));
1602 mtf->rrc.d = (int*) calloc(1,NumResResContacts * sizeof(int));
1603 mtf->rpc.n = NumResPepContacts;
1604 mtf->rpc.r1 = (int*) calloc(1,NumResPepContacts * sizeof(int));
1605 mtf->rpc.p2 = (int*) calloc(1,NumResPepContacts * sizeof(int));
1606 mtf->rpc.d = (int*) calloc(1,NumResPepContacts * sizeof(int));
1607 mtf->mll = (int**) calloc(1,NumResidues * sizeof(int*));
1608 for (i=0; i<NumResidues; i++) {
1609 mtf->mll[i] = (int*) calloc(1,NumResidues * sizeof(int));
1610 }
1611 return(mtf);
1612 }
1613
1614
FreeFldMtf(Fld_Mtf * mtf)1615 Fld_Mtf* FreeFldMtf(Fld_Mtf* mtf) {
1616 /*----------------------------------------------------------*/
1617 /* free contact lists and loop lengths */
1618 /*----------------------------------------------------------*/
1619 int i;
1620
1621 free(mtf->rrc.r1);
1622 free(mtf->rrc.r2);
1623 free(mtf->rrc.d);
1624 free(mtf->rpc.r1);
1625 free(mtf->rpc.p2);
1626 free(mtf->rpc.d);
1627 for (i=0; i<mtf->n; i++) {
1628 free(mtf->mll[i]);
1629 }
1630 free(mtf->mll);
1631 free(mtf);
1632 return NULL;
1633 }
1634
1635
PrintFldMtf(Fld_Mtf * mtf,FILE * pFile)1636 void PrintFldMtf(Fld_Mtf* mtf, FILE* pFile) {
1637 /*----------------------------------------------------------*/
1638 /* print contact lists and loop lengths. */
1639 /*----------------------------------------------------------*/
1640 int i, j;
1641
1642 fprintf(pFile, "Contact Lists, Loop Lengths:\n");
1643 fprintf(pFile, "Number of residues: %4d\n", mtf->n);
1644 fprintf(pFile, "Number of residue-residue contacts: %6d\n", mtf->rrc.n);
1645 fprintf(pFile, "res index 1, res index 2, Distance interval\n");
1646 for (i=0; i<mtf->rrc.n; i++) {
1647 fprintf(pFile, " %9d %12d %18d\n", mtf->rrc.r1[i], mtf->rrc.r2[i], mtf->rrc.d[i]);
1648 }
1649 fprintf(pFile, "Number of residue-peptide contacts: %6d\n", mtf->rpc.n);
1650 fprintf(pFile, "res index, peptide index, Distance interval\n");
1651 for (i=0; i<mtf->rpc.n; i++) {
1652 fprintf(pFile, " %7d %14d %18d\n", mtf->rpc.r1[i], mtf->rpc.p2[i], mtf->rpc.d[i]);
1653 }
1654 fprintf(pFile, "Minimum loop lengths:\n");
1655 for (i=0; i<mtf->n; i++) {
1656 for (j=0; j<mtf->n; j++) {
1657 fprintf(pFile, " %4d", mtf->mll[i][j]);
1658 }
1659 fprintf(pFile, "\n");
1660 }
1661 fprintf(pFile, "\n");
1662 }
1663
1664
NewThdTbl(int NumResults,int NumCoreElements)1665 Thd_Tbl* NewThdTbl(int NumResults, int NumCoreElements) {
1666 /*----------------------------------------------------------*/
1667 /* allocate space for results. */
1668 /*----------------------------------------------------------*/
1669 int i;
1670 Thd_Tbl* ttb;
1671
1672 ttb = (Thd_Tbl*) calloc(1,sizeof(Thd_Tbl));
1673 ttb->n = NumResults;
1674 ttb->nsc = NumCoreElements;
1675 ttb->tg = (float*) calloc(1,NumResults * sizeof(float));
1676 ttb->ps = (float*) calloc(1,NumResults * sizeof(float));
1677 ttb->ms = (float*) calloc(1,NumResults * sizeof(float));
1678 ttb->cs = (float*) calloc(1,NumResults * sizeof(float));
1679 ttb->lps = (float*) calloc(1,NumResults * sizeof(float));
1680 ttb->zsc = (float*) calloc(1,NumResults * sizeof(float));
1681 ttb->g0 = (float*) calloc(1,NumResults * sizeof(float));
1682 ttb->m0 = (float*) calloc(1,NumResults * sizeof(float));
1683 ttb->errm = (float*) calloc(1,NumResults * sizeof(float));
1684 ttb->errp = (float*) calloc(1,NumResults * sizeof(float));
1685 ttb->tf = (int*) calloc(1,NumResults * sizeof(int));
1686 ttb->ts = (int*) calloc(1,NumResults * sizeof(int));
1687 ttb->ls = (int*) calloc(1,NumResults * sizeof(int));
1688 ttb->pr = (int*) calloc(1,NumResults * sizeof(int));
1689 ttb->nx = (int*) calloc(1,NumResults * sizeof(int));
1690 ttb->al = (int**) calloc(1,NumCoreElements * sizeof(int*));
1691 ttb->no = (int**) calloc(1,NumCoreElements * sizeof(int*));
1692 ttb->co = (int**) calloc(1,NumCoreElements * sizeof(int*));
1693 for (i=0; i<NumCoreElements; i++) {
1694 ttb->al[i] = (int*) calloc(1,NumResults *sizeof(int));
1695 ttb->no[i] = (int*) calloc(1,NumResults *sizeof(int));
1696 ttb->co[i] = (int*) calloc(1,NumResults *sizeof(int));
1697 }
1698 return(ttb);
1699 }
1700
1701
FreeThdTbl(Thd_Tbl * ttb)1702 Thd_Tbl* FreeThdTbl(Thd_Tbl* ttb) {
1703 /*----------------------------------------------------------*/
1704 /* free results. */
1705 /*----------------------------------------------------------*/
1706 int i;
1707
1708 free(ttb->tg);
1709 free(ttb->ps);
1710 free(ttb->ms);
1711 free(ttb->cs);
1712 free(ttb->lps);
1713 free(ttb->zsc);
1714 free(ttb->g0);
1715 free(ttb->m0);
1716 free(ttb->errm);
1717 free(ttb->errp);
1718 free(ttb->tf);
1719 free(ttb->ts);
1720 free(ttb->ls);
1721 free(ttb->pr);
1722 free(ttb->nx);
1723 for (i=0; i<ttb->nsc; i++) {
1724 free(ttb->al[i]);
1725 free(ttb->no[i]);
1726 free(ttb->co[i]);
1727 }
1728 free(ttb->al);
1729 free(ttb->no);
1730 free(ttb->co);
1731 free(ttb);
1732 return NULL;
1733 }
1734
1735
PrintThdTbl(Thd_Tbl * ttb,FILE * pFile)1736 void PrintThdTbl(Thd_Tbl* ttb, FILE* pFile) {
1737 /*------------------------------------------------------------------------------*/
1738 /* print the results. */
1739 /*------------------------------------------------------------------------------*/
1740 int i, j;
1741
1742 fprintf(pFile, "Threading Results:\n");
1743 fprintf(pFile, "number of threads: %6d\n", ttb->n);
1744 fprintf(pFile, "number of core segments: %6d\n", ttb->nsc);
1745 fprintf(pFile, "index of lowest energy thread: %6d\n", ttb->mn);
1746 fprintf(pFile, "index of highest energy thread: %6d\n", ttb->mx);
1747 fprintf(pFile, "for each thread:\n");
1748 fprintf(pFile, " tg ps ms cs lps zsc\n");
1749 for (i=0; i<ttb->n; i++) {
1750 fprintf(pFile, " %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
1751 ttb->tg[i], ttb->ps[i], ttb->ms[i], ttb->cs[i], ttb->lps[i], ttb->zsc[i]);
1752 }
1753 fprintf(pFile, " g0 m0 errm errp\n");
1754 for (i=0; i<ttb->n; i++) {
1755 fprintf(pFile, " %12.5e %12.5e %12.5e %12.5e\n",
1756 ttb->g0[i], ttb->m0[i], ttb->errm[i], ttb->errp[i]);
1757 }
1758 fprintf(pFile, " tf ts ls pr nx\n");
1759 for (i=0; i<ttb->n; i++) {
1760 fprintf(pFile, " %6d %6d %6d %6d %6d\n",
1761 ttb->tf[i], ttb->ts[i], ttb->ls[i], ttb->pr[i], ttb->nx[i]);
1762 }
1763 fprintf(pFile, "threading alignments: %6d\n", i+1);
1764 fprintf(pFile, "Centers:\n");
1765 for (i=0; i<ttb->n; i++) {
1766 for (j=0; j<ttb->nsc; j++) {
1767 fprintf(pFile, " %4d", ttb->al[j][i] + 1);
1768 }
1769 fprintf(pFile, "\n");
1770 }
1771 fprintf(pFile, "N offsets:\n");
1772 for (i=0; i<ttb->n; i++) {
1773 for (j=0; j<ttb->nsc; j++) {
1774 fprintf(pFile, " %4d", ttb->no[j][i]);
1775 }
1776 fprintf(pFile, "\n");
1777 }
1778 fprintf(pFile, "C offsets:\n");
1779 for (i=0; i<ttb->n; i++) {
1780 for (j=0; j<ttb->nsc; j++) {
1781 fprintf(pFile, " %4d", ttb->co[j][i]);
1782 }
1783 fprintf(pFile, "\n");
1784 }
1785 fprintf(pFile, "\n");
1786 }
1787
1788
CopyResult(Thd_Tbl * pFromResults,Thd_Tbl * pToResults,int from,int to)1789 int CopyResult(Thd_Tbl* pFromResults, Thd_Tbl* pToResults, int from, int to) {
1790 /*------------------------------------------------------------------------------*/
1791 /* copy the from-th result of pFromResults to the to-th result of pToResults. */
1792 /*------------------------------------------------------------------------------*/
1793 int i;
1794
1795 pToResults->tg[to] = pFromResults->tg[from];
1796 pToResults->ps[to] = pFromResults->ps[from];
1797 pToResults->ms[to] = pFromResults->ms[from];
1798 pToResults->cs[to] = pFromResults->cs[from];
1799 pToResults->lps[to] = pFromResults->lps[from];
1800 pToResults->zsc[to] = pFromResults->zsc[from];
1801 pToResults->g0[to] = pFromResults->g0[from];
1802 pToResults->m0[to] = pFromResults->m0[from];
1803 pToResults->errm[to] = pFromResults->errm[from];
1804 pToResults->errp[to] = pFromResults->errp[from];
1805 pToResults->tf[to] = pFromResults->tf[from];
1806 pToResults->ts[to] = pFromResults->ts[from];
1807 pToResults->ls[to] = pFromResults->ls[from];
1808 pToResults->pr[to] = pFromResults->pr[from];
1809 pToResults->nx[to] = pFromResults->nx[from];
1810 for (i=0; i<pFromResults->nsc; i++) {
1811 pToResults->al[i][to] = pFromResults->al[i][from];
1812 pToResults->no[i][to] = pFromResults->no[i][from];
1813 pToResults->co[i][to] = pFromResults->co[i][from];
1814 }
1815 return(1);
1816 }
1817
1818
OrderThdTbl(Thd_Tbl * pResults)1819 void OrderThdTbl(Thd_Tbl* pResults) {
1820 /*----------------------------------------------------------*/
1821 /* put the results in order of highest score to lowest. */
1822 /* it's O(n**2), but that's ok because n is very small. */
1823 /* return the new ordered results, free original results. */
1824 /*----------------------------------------------------------*/
1825 int* Order;
1826 Boolean* CheckList;
1827 int i, j, SaveIndex, Index;
1828 double HighestScore;
1829 Thd_Tbl* pOrderedResults;
1830
1831 /* mem allocation for array where order of results is stored */
1832 Order = (int*) calloc(1,sizeof(int) * pResults->n);
1833 /* mem allocation for checklist to tell which results have been ordered */
1834 CheckList = (Boolean*) calloc(1,sizeof(Boolean) * pResults->n);
1835
1836 /* look through the list n times */
1837 for (i=0; i<pResults->n; i++) {
1838 HighestScore = -9999999999.0;
1839 SaveIndex = -1;
1840 /* each time, find the highest remaining score */
1841 for (j=0; j<pResults->n; j++) {
1842 if (CheckList[j] == FALSE) {
1843 if (pResults->tg[j] > HighestScore) {
1844 HighestScore = pResults->tg[j];
1845 SaveIndex = j;
1846 }
1847 }
1848 }
1849 /* save the index, and mark that the score has been ordered */
1850 Order[i] = SaveIndex;
1851 CheckList[SaveIndex] = TRUE;
1852 }
1853
1854 /* now that the order is determined, do the re-ordering */
1855 pOrderedResults = NewThdTbl(pResults->n, pResults->nsc);
1856 for (i=0; i<pResults->n; i++) {
1857 CopyResult(pResults, pOrderedResults, Order[i], i);
1858 }
1859 pOrderedResults->mx = 0;
1860 pOrderedResults->mn = pResults->n - 1;
1861
1862 /* now copy results back to original structure */
1863 for (i=0; i<pResults->n; i++) {
1864 CopyResult(pOrderedResults, pResults, i, i);
1865 }
1866
1867 /* update next and previous arrays to match new order */
1868 for (i=0; i<pResults->n; i++) {
1869 /* point to next-lower energy thread */
1870 pResults->nx[i] = i+1 < pResults->n ? i+1 : 0;
1871 /* point to next-higher energy thread */
1872 Index = pResults->n - i - 1;
1873 pResults->pr[Index] = Index-1 < 0 ? 0 : Index-1;
1874 }
1875
1876 /* update min and max indices to match new order */
1877 pResults->mx = pOrderedResults->mx;
1878 pResults->mn = pOrderedResults->mn;
1879
1880 /* free memory allocated for this routine */
1881 free(Order);
1882 free(CheckList);
1883 FreeThdTbl(pOrderedResults);
1884
1885 return;
1886 }
1887
1888
ScaleThdTbl(Thd_Tbl * ttb,double ScalingFactor)1889 void ScaleThdTbl(Thd_Tbl* ttb, double ScalingFactor) {
1890 /*----------------------------------------------------------*/
1891 /* scale the energies in ttb by ScalingFactor */
1892 /*----------------------------------------------------------*/
1893 int i;
1894 float SF;
1895
1896 SF = (float) ScalingFactor;
1897
1898 for (i=0; i<ttb->n; i++) {
1899 ttb->tg[i] *= SF;
1900 ttb->ps[i] *= SF;
1901 ttb->ms[i] *= SF;
1902 ttb->cs[i] *= SF;
1903 ttb->lps[i] *= SF;
1904 ttb->zsc[i] *= SF;
1905 ttb->g0[i] *= SF;
1906 ttb->m0[i] *= SF;
1907 ttb->errm[i] *= SF;
1908 ttb->errp[i] *= SF;
1909 }
1910 }
1911
1912
ThrdRound(double Num)1913 int ThrdRound(double Num) {
1914 /*----------------------------------------------------------*/
1915 /* round fp to nearest int */
1916 /*----------------------------------------------------------*/
1917 if (Num > 0) {
1918 return((int)(Num + 0.5));
1919 }
1920 else {
1921 return((int)(Num - 0.5));
1922 }
1923 }
1924