libs/dtgf/detrend_code_C99.h

// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/dtgf/detrend_code_C99.h,v 1.1 2013/04/28 07:46:58 tplarson Exp $

#include "ctypes_def.h"

#if TYPE == FLOAT 
   #define DETREND sdetrend
   #define DETREND_DISCONTIG sdetrend_discontig
   #define LEGENDRE_POLY slegendre_poly
   #define ORTHO_POLY sortho_poly
#elif TYPE == DOUBLE 
   #define DETREND ddetrend
   #define DETREND_DISCONTIG ddetrend_discontig
   #define LEGENDRE_POLY dlegendre_poly
   #define ORTHO_POLY dortho_poly
#endif

//#define DEBUGFILE

static void LEGENDRE_POLY(int degree, int length, RTYPE **poly);
static void ORTHO_POLY(int degree, int length, RTYPE **poly, int *isgood);


void DETREND_DISCONTIG( int n, RTYPE *data, int *isgood, int degree, 
                         int length, int skip, int m, int *sect, int detrend_first)
{
  int i,first,last;
  // Seperately detrend m sections of data separated by 
  // discontinuities.
  // The first section begins on data[0], 
  // the m-1 first sections end on data[sect_last[i]], i=0,1,...,m-2,
  // and the last section ends on data[n-1].

  if (m==0)
    DETREND(n,data,isgood,degree,length,skip,detrend_first);
  else
  {
    for (i=0; i<m; i++)
    {
      first = sect[2*i];
      last = sect[2*i+1];
#ifdef DEBUGOUT
      printf("Detrending [%d:%d]\n",first,last);
#endif
      DETREND(last-first+1,&data[first],&isgood[first],degree,length,skip,detrend_first);
    }
  }
}
                

void DETREND( int n, RTYPE *data, int *isgood, int degree, 
              int length, int skip, int detrend_first)
{
  int first,i,d, overlap, start,end;
  RTYPE **legendre_poly, **gap_poly, *fit, gamma, *apodizer;
#ifdef DEBUGFILE
  char name[20];
  static int count=0;
  FILE *fh;
#endif

#ifdef DEBUGOUT
  printf("%s: n=%d, length=%d, skip=%d, degree=%d\n",
           __func__,n,length,skip,degree);
#endif

  fit = (RTYPE *)calloc(n,sizeof(RTYPE));
  // Set up Legendre polynomials
  legendre_poly = (RTYPE **)malloc((degree+1)*sizeof(RTYPE *));
  gap_poly = (RTYPE **)malloc((degree+1)*sizeof(RTYPE *));
  for (i=0; i<=degree; i++)
  {
    legendre_poly[i] = (RTYPE *)malloc(2*length*sizeof(RTYPE));
    gap_poly[i] = (RTYPE *)malloc(2*length*sizeof(RTYPE));
  }
  LEGENDRE_POLY(degree,length, legendre_poly);

  // Set up apodization to stitch together individual fits.
  overlap = length-skip;
  apodizer = (RTYPE *)malloc(overlap*sizeof(RTYPE));
  for (i=0;i<overlap;i++)
  {
    apodizer[i] = cos(M_PI/2*((RTYPE) i+1)/(overlap+1));
    apodizer[i] *= apodizer[i];
  }

//  for (first=0; first<=n-length; first+=skip)
  first=detrend_first;
  do
  {
/*
    // Possibly grow the last interval to match n.
    if (first+2*length>n )
    {
      length = n-first;
      LEGENDRE_POLY(degree, length, legendre_poly);
    }
*/
// don't grow it quite so much:
    if (first+3*length/2>n)
    {
      length = n-first;
      LEGENDRE_POLY(degree, length, legendre_poly);
    }

    start = first;
    while (start<first+length && !isgood[start])
      start++;
    if (start==first+length)
    {
      /* An empty interval. fill with zeros and skip. */
      for (i=first; i<first+length; i++)
        fit[i] = 0;
//      continue;
// continue is not correct here because it skips the increment
      goto emptyinterval;
    }

    end = first+length-1;
    while (end>start && !isgood[end])
      end--;
    end++;

    //    printf("first = %d, start = %d, end = %d.\n",first,start,end);
    if ((end-start) < length/2)
    {
#ifdef DEBUGOUT
      printf("special treatment of [%d:%d], length = %d\n",start,end,length);
#endif
      LEGENDRE_POLY(degree, (end-start), legendre_poly);
      // Modify polynomials to be orthonormal wrt. 
      // the local gap structure.
      for (i=0; i<=degree; i++)
        memcpy(gap_poly[i],legendre_poly[i],(end-start)*sizeof(RTYPE));
      ORTHO_POLY(degree,(end-start), gap_poly, &isgood[start]);
      for (i=start; i<end; i++)
        fit[i] = 0;
      for (d=0;d<=degree;d++)
      {
        gamma = 0;
        for (i=start; i<end; i++)
          if (isgood[i])
            gamma += gap_poly[d][i-start]*data[i];
        for (i=start; i<end; i++)
          fit[i] += gamma*gap_poly[d][i-start];
      }
      LEGENDRE_POLY(degree, length, legendre_poly);      
    }
    else
    {
#ifdef DEBUGOUT
      printf("normal treatment [%d:%d]\n",first,first+length-1);
#endif
      // Modify polynomials to be orthonormal wrt. 
      // the local gap structure.
      for (i=0; i<=degree; i++)
        memcpy(gap_poly[i],legendre_poly[i],length*sizeof(RTYPE));
      ORTHO_POLY(degree,length, gap_poly, &isgood[first]);

      // Subtract polynomial components one by one.
      if (first!=detrend_first)
      {
        for (i=first; i<first+overlap; i++)
        {
          RTYPE alpha = apodizer[i-first];        
          fit[i] = alpha*fit[i];
        }
      }
      for (d=0;d<=degree;d++)
      {
        gamma = 0;
        for (i=first; i<first+length; i++)
          if (isgood[i])
            gamma += gap_poly[d][i-first]*data[i];
        if (first!=detrend_first)
          for (i=first; i<first+overlap; i++)
          {
            RTYPE alpha = apodizer[i-first];      
            fit[i] += (1-alpha)*gamma*gap_poly[d][i-first];
          }
        else
          i=first;
        for (; i<first+length; i++)
          fit[i] += gamma*gap_poly[d][i-first];
      } 
    }
  emptyinterval:
  first += skip;
//  } while (first<=n-length);
  } while (first<=n-length/2);

  for (i=0; i<n; i++)
    if (!isgood[i])      
      fit[i] = 0;
  
  for (i=0; i<n; i++)
    if (isgood[i])      
      data[i] -= fit[i];
    else
      data[i] = 0;

#ifdef DEBUGFILE
  printf("writing debug file #%d\n",count);
  sprintf(name,"detrend_debug_%d.out",count++);
  fh = fopen(name,"w");
  assert(fh);
  for (i=0; i<n; i++)
    fprintf(fh,"%e %e\n",fit[i],data[i]);
  fclose(fh);
#endif

  for (i=0; i<=degree; i++)
  {
    free(legendre_poly[i]);
    free(gap_poly[i]);
  }
  free(legendre_poly);
  free(gap_poly);
  free(fit);  
  free(apodizer);
}


static void LEGENDRE_POLY(int degree, int length, RTYPE **poly)
{
  int d,i;
 
  // Set up Legendre polynomials using recurrence.
  for (i=0;i<(length+1)/2; i++)
  {
    poly[0][i] = 1.0;
    poly[1][i] = -1.0+(2.0*(i+1))/(length+1);
  }  
  for (d=1; d<degree; d++)
    for (i=0; i<(length+1)/2; i++)
      poly[d+1][i] = ((2*d+1)*poly[1][i]*poly[d][i]-d*poly[d-1][i])/(d+1);

  // Set up symmetric / anti-symmetric part.
  for (d=1; d<=degree; d+=2)
    for (i=0; i<length/2; i++)
      poly[d][length-i-1] = -poly[d][i];
  for (d=0; d<=degree; d+=2)
    for (i=0; i<length/2; i++)
      poly[d][length-i-1] = poly[d][i];
}


static void ORTHO_POLY(int degree, int length, RTYPE **poly, int *isgood)
{
  int d,i,j;
  RTYPE norm;
 
  // Normalize wrt. inner product on good part.
  for (d=0; d<=degree; d++)
  {
    norm = 0.0;
    for (i=0; i<length; i++)
      if (isgood[i])
        norm +=  poly[d][i]*poly[d][i];
    norm = sqrt(norm);
    for (i=0; i<length; i++)
      poly[d][i] /= norm;
  }

  // Orthogonalize  wrt. inner product on good part.
  for (d=0; d<=degree; d++)    
    for (j=d-1; j>=0; j--)
    {
      norm = 0;
      for (i=0; i<length; i++)
        if (isgood[i])
          norm +=  poly[j][i]*poly[d][i];
      for (i=0; i<length; i++)
        poly[d][i] -=  norm*poly[j][i];
    }

  // Normalize again wrt. inner product on good part.
  for (d=0; d<=degree; d++)
  {
    norm = 0.0;
    for (i=0; i<length; i++)
      if (isgood[i])
        norm +=  poly[d][i]*poly[d][i];
    norm = sqrt(norm);
    for (i=0; i<length; i++)
      poly[d][i] /= norm;
  }

}


#undef DETREND
#undef DETREND_DISCONTIG
#undef LEGENDRE_POLY
#undef ORTHO_POLY
#include "ctypes_undef.h"
Revision:	1.2
Committed:	Thu Jan 16 08:14:34 2014 UTC (9 years, 8 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	globalhs_version_5, Ver_8-7, Ver_8-5, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_8-3, globalhs_version_8, globalhs_version_9, globalhs_version_2, globalhs_version_3, globalhs_version_4, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-10, Ver_8-6, Ver_9-1, Ver_8-4, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.1:	+15 -38 lines
Log Message:	adjust behavior at end of detrending interval
#	Content
1	// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/dtgf/detrend_code_C99.h,v 1.1 2013/04/28 07:46:58 tplarson Exp $
2
3	#include "ctypes_def.h"
4
5	#if TYPE == FLOAT
6	#define DETREND sdetrend
7	#define DETREND_DISCONTIG sdetrend_discontig
8	#define LEGENDRE_POLY slegendre_poly
9	#define ORTHO_POLY sortho_poly
10	#elif TYPE == DOUBLE
11	#define DETREND ddetrend
12	#define DETREND_DISCONTIG ddetrend_discontig
13	#define LEGENDRE_POLY dlegendre_poly
14	#define ORTHO_POLY dortho_poly
15	#endif
16
17	//#define DEBUGFILE
18
19	static void LEGENDRE_POLY(int degree, int length, RTYPE **poly);
20	static void ORTHO_POLY(int degree, int length, RTYPE *poly, int isgood);
21
22
23	void DETREND_DISCONTIG( int n, RTYPE data, int isgood, int degree,
24	int length, int skip, int m, int *sect, int detrend_first)
25	{
26	int i,first,last;
27	// Seperately detrend m sections of data separated by
28	// discontinuities.
29	// The first section begins on data[0],
30	// the m-1 first sections end on data[sect_last[i]], i=0,1,...,m-2,
31	// and the last section ends on data[n-1].
32
33	if (m==0)
34	DETREND(n,data,isgood,degree,length,skip,detrend_first);
35	else
36	{
37	for (i=0; i<m; i++)
38	{
39	first = sect[2*i];
40	last = sect[2*i+1];
41	#ifdef DEBUGOUT
42	printf("Detrending [%d:%d]\n",first,last);
43	#endif
44	DETREND(last-first+1,&data[first],&isgood[first],degree,length,skip,detrend_first);
45	}
46	}
47	}
48
49
50
51	void DETREND( int n, RTYPE data, int isgood, int degree,
52	int length, int skip, int detrend_first)
53	{
54	int first,i,d, overlap, start,end;
55	RTYPE legendre_poly, gap_poly, fit, gamma, apodizer;
56	#ifdef DEBUGFILE
57	char name[20];
58	static int count=0;
59	FILE *fh;
60	#endif
61
62	#ifdef DEBUGOUT
63	printf("%s: n=%d, length=%d, skip=%d, degree=%d\n",
64	__func__,n,length,skip,degree);
65	#endif
66
67	fit = (RTYPE *)calloc(n,sizeof(RTYPE));
68	// Set up Legendre polynomials
69	legendre_poly = (RTYPE *)malloc((degree+1)sizeof(RTYPE *));
70	gap_poly = (RTYPE *)malloc((degree+1)sizeof(RTYPE *));
71	for (i=0; i<=degree; i++)
72	{
73	legendre_poly[i] = (RTYPE )malloc(2length*sizeof(RTYPE));
74	gap_poly[i] = (RTYPE )malloc(2length*sizeof(RTYPE));
75	}
76	LEGENDRE_POLY(degree,length, legendre_poly);
77
78	// Set up apodization to stitch together individual fits.
79	overlap = length-skip;
80	apodizer = (RTYPE )malloc(overlapsizeof(RTYPE));
81	for (i=0;i<overlap;i++)
82	{
83	apodizer[i] = cos(M_PI/2*((RTYPE) i+1)/(overlap+1));
84	apodizer[i] *= apodizer[i];
85	}
86
87	// for (first=0; first<=n-length; first+=skip)
88	first=detrend_first;
89	do
90	{
91	/*
92	// Possibly grow the last interval to match n.
93	if (first+2*length>n )
94	{
95	length = n-first;
96	LEGENDRE_POLY(degree, length, legendre_poly);
97	}
98	*/
99	// don't grow it quite so much:
100	if (first+3*length/2>n)
101	{
102	length = n-first;
103	LEGENDRE_POLY(degree, length, legendre_poly);
104	}
105
106	start = first;
107	while (start<first+length && !isgood[start])
108	start++;
109	if (start==first+length)
110	{
111	/* An empty interval. fill with zeros and skip. */
112	for (i=first; i<first+length; i++)
113	fit[i] = 0;
114	// continue;
115	// continue is not correct here because it skips the increment
116	goto emptyinterval;
117	}
118
119	end = first+length-1;
120	while (end>start && !isgood[end])
121	end--;
122	end++;
123
124	// printf("first = %d, start = %d, end = %d.\n",first,start,end);
125	if ((end-start) < length/2)
126	{
127	#ifdef DEBUGOUT
128	printf("special treatment of [%d:%d], length = %d\n",start,end,length);
129	#endif
130	LEGENDRE_POLY(degree, (end-start), legendre_poly);
131	// Modify polynomials to be orthonormal wrt.
132	// the local gap structure.
133	for (i=0; i<=degree; i++)
134	memcpy(gap_poly[i],legendre_poly[i],(end-start)*sizeof(RTYPE));
135	ORTHO_POLY(degree,(end-start), gap_poly, &isgood[start]);
136	for (i=start; i<end; i++)
137	fit[i] = 0;
138	for (d=0;d<=degree;d++)
139	{
140	gamma = 0;
141	for (i=start; i<end; i++)
142	if (isgood[i])
143	gamma += gap_poly[d][i-start]*data[i];
144	for (i=start; i<end; i++)
145	fit[i] += gamma*gap_poly[d][i-start];
146	}
147	LEGENDRE_POLY(degree, length, legendre_poly);
148	}
149	else
150	{
151	#ifdef DEBUGOUT
152	printf("normal treatment [%d:%d]\n",first,first+length-1);
153	#endif
154	// Modify polynomials to be orthonormal wrt.
155	// the local gap structure.
156	for (i=0; i<=degree; i++)
157	memcpy(gap_poly[i],legendre_poly[i],length*sizeof(RTYPE));
158	ORTHO_POLY(degree,length, gap_poly, &isgood[first]);
159
160	// Subtract polynomial components one by one.
161	if (first!=detrend_first)
162	{
163	for (i=first; i<first+overlap; i++)
164	{
165	RTYPE alpha = apodizer[i-first];
166	fit[i] = alpha*fit[i];
167	}
168	}
169	for (d=0;d<=degree;d++)
170	{
171	gamma = 0;
172	for (i=first; i<first+length; i++)
173	if (isgood[i])
174	gamma += gap_poly[d][i-first]*data[i];
175	if (first!=detrend_first)
176	for (i=first; i<first+overlap; i++)
177	{
178	RTYPE alpha = apodizer[i-first];
179	fit[i] += (1-alpha)gammagap_poly[d][i-first];
180	}
181	else
182	i=first;
183	for (; i<first+length; i++)
184	fit[i] += gamma*gap_poly[d][i-first];
185	}
186	}
187	emptyinterval:
188	first += skip;
189	// } while (first<=n-length);
190	} while (first<=n-length/2);
191
192	for (i=0; i<n; i++)
193	if (!isgood[i])
194	fit[i] = 0;
195
196	for (i=0; i<n; i++)
197	if (isgood[i])
198	data[i] -= fit[i];
199	else
200	data[i] = 0;
201
202	#ifdef DEBUGFILE
203	printf("writing debug file #%d\n",count);
204	sprintf(name,"detrend_debug_%d.out",count++);
205	fh = fopen(name,"w");
206	assert(fh);
207	for (i=0; i<n; i++)
208	fprintf(fh,"%e %e\n",fit[i],data[i]);
209	fclose(fh);
210	#endif
211
212	for (i=0; i<=degree; i++)
213	{
214	free(legendre_poly[i]);
215	free(gap_poly[i]);
216	}
217	free(legendre_poly);
218	free(gap_poly);
219	free(fit);
220	free(apodizer);
221	}
222
223
224	static void LEGENDRE_POLY(int degree, int length, RTYPE **poly)
225	{
226	int d,i;
227
228	// Set up Legendre polynomials using recurrence.
229	for (i=0;i<(length+1)/2; i++)
230	{
231	poly[0][i] = 1.0;
232	poly[1][i] = -1.0+(2.0*(i+1))/(length+1);
233	}
234	for (d=1; d<degree; d++)
235	for (i=0; i<(length+1)/2; i++)
236	poly[d+1][i] = ((2d+1)poly[1][i]poly[d][i]-dpoly[d-1][i])/(d+1);
237
238	// Set up symmetric / anti-symmetric part.
239	for (d=1; d<=degree; d+=2)
240	for (i=0; i<length/2; i++)
241	poly[d][length-i-1] = -poly[d][i];
242	for (d=0; d<=degree; d+=2)
243	for (i=0; i<length/2; i++)
244	poly[d][length-i-1] = poly[d][i];
245	}
246
247
248	static void ORTHO_POLY(int degree, int length, RTYPE *poly, int isgood)
249	{
250	int d,i,j;
251	RTYPE norm;
252
253	// Normalize wrt. inner product on good part.
254	for (d=0; d<=degree; d++)
255	{
256	norm = 0.0;
257	for (i=0; i<length; i++)
258	if (isgood[i])
259	norm += poly[d][i]*poly[d][i];
260	norm = sqrt(norm);
261	for (i=0; i<length; i++)
262	poly[d][i] /= norm;
263	}
264
265	// Orthogonalize wrt. inner product on good part.
266	for (d=0; d<=degree; d++)
267	for (j=d-1; j>=0; j--)
268	{
269	norm = 0;
270	for (i=0; i<length; i++)
271	if (isgood[i])
272	norm += poly[j][i]*poly[d][i];
273	for (i=0; i<length; i++)
274	poly[d][i] -= norm*poly[j][i];
275	}
276
277	// Normalize again wrt. inner product on good part.
278	for (d=0; d<=degree; d++)
279	{
280	norm = 0.0;
281	for (i=0; i<length; i++)
282	if (isgood[i])
283	norm += poly[d][i]*poly[d][i];
284	norm = sqrt(norm);
285	for (i=0; i<length; i++)
286	poly[d][i] /= norm;
287	}
288
289	}
290
291
292	#undef DETREND
293	#undef DETREND_DISCONTIG
294	#undef LEGENDRE_POLY
295	#undef ORTHO_POLY
296	#include "ctypes_undef.h"