src/Objects/DataManip/ctfInfo/src/org/powerspectraTangentLine.c.20150124

   1 /*
   2 # powerspectraTangentLine : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : powerspectraTangentLine
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 /* $Log$ */
  12 static char __sccs_id[] = "%Z%lctfDetermination ver%I%; Date:%D% %Z%";
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 #include "./powerspectraTangentLine.h"
  17 #include "./ctfInfoWrite2.h"
  18 #define DEBUG
  19 #include "genUtil.h"
  20 #include "Vector.h"
  21 #include "nr2.h"
  22 #include "Memory.h"
  23
  24 typedef struct stepInfo{
  25         float min;
  26         float max;
  27         float step;
  28         long  Number;
  29 }stepInfo;
  30
  31 typedef struct iterationInfo{
  32         long  degree;
  33         long  CutLowSize;
  34         long  CutHighSize;
  35         float CutLine;
  36         float Emin;
  37         float Emax;
  38 }iterationInfo;
  39
  40
  41 double __power(double x,int a)
  42 {
  43         int i;
  44         double y = 1.0;
  45
  46         for(i=0 ; i<a ; i++){
  47                 y *= x;
  48         }
  49         return(y);
  50 }
  51
  52 void __stepNumberCulc(stepInfo *sInfo, long degree)
  53 {
  54         long i;
  55         for(i=0 ; i<= degree ; i++){
  56                 sInfo[i].Number =  (long) ((sInfo[i].max - sInfo[i].min )/ sInfo[i].step);
  57                 if(sInfo[i].Number < 1){
  58                         sInfo[i].Number = 1;
  59                 }
  60         }
  61 }
  62
  63 /*mode = 1  belowLine
  64        = -1 aboveLine*/
  65 void
  66 __CoefficientDetermine(floatVector* spacing ,floatVector* Plog ,floatVector* tmpcoefficient ,stepInfo* sInfo,iterationInfo* iInfo, float mode)
  67 {
  68         long i,j,k,l;
  69         long AdvanceRatio = 1;
  70         long Number[iInfo->degree+1];
  71         unsigned long long int  belowN;
  72         unsigned long long int fewerN;
  73         unsigned long n;
  74         unsigned long long int  sum = 1;
  75         float coefficient[iInfo->degree+1];
  76         float E;
  77         float minTangent;
  78         floatVector* Tangent;
  79
  80 /*      if(mode == 1.0){
  81                 fprintf(stderr,"BelowLineEstimation ");
  82         }
  83         else if(mode == -1.0){
  84                 fprintf(stderr,"AboveLineEstimation ");
  85         }
  86         else{
  87                 fprintf(stderr,"CoefficientDetermine don't support this mode\nprogram stop");
  88                 exit(1);
  89         }
  90
  91         fprintf(stderr,"degree:%2ld\n",iInfo->degree);
  92 */
  93         Tangent = floatVectorInit(NULL,spacing->size);
  94
  95         iInfo-> Emax;
  96         for(j=0 ; j<=iInfo->degree ; j++){
  97                 sum *= sInfo[j].Number;
  98                 //fprintf(stderr,"StepCondition[%2ld] Min: %+8.3lf Max: %+8.3lf StepSize: %6.3f StepNumber: %4ld\n",j,sInfo[j].min,sInfo[j].max,sInfo[j].step,sInfo[j].Number);
  99         }
 100         //fprintf(stderr,"sum:%llu\n",sum);
 101
 102         for(n=0 ; n< sum; n++){
 103                 E = 0.0;
 104                 minTangent = 0.0;
 105                 for(j=0 ; j<=iInfo->degree ; j++){
 106                         fewerN = 1;
 107                         for(k=0 ; k<j ; k++){
 108                                 fewerN *= sInfo[k].Number;
 109                         }
 110                         belowN = fewerN * sInfo[j].Number;
 111                         coefficient[j] = (long)((n % belowN)/fewerN)* sInfo[j].step + sInfo[j].min;
 112                         minTangent += coefficient[j]*__power(spacing->data[iInfo->CutHighSize],j);
 113                 }
 114
 115                 if(minTangent < iInfo->CutLine){
 116                         E = iInfo->Emax;
 117                 }
 118
 119                 for(i=iInfo->CutLowSize ; i<=iInfo->CutHighSize ; i++){
 120                         Tangent->data[i] = coefficient[0];
 121                         for(j=1 ; j<=iInfo->degree ; j++){
 122                                 Tangent->data[i] += coefficient[j]*__power(spacing->data[i],j);
 123                         }
 124                         E += fabs(Tangent->data[i] - Plog->data[i]);
 125
 126                         if((Tangent->data[i] - Plog->data[i])*mode > 0 || Tangent->data[i] < iInfo->CutLine || E>iInfo->Emin){
 127                                 E = iInfo->Emax;
 128                                 break;
 129                         }
 130                 }
 131
 132
 133                 if(E <= iInfo->Emin){
 134                         for(j=0 ; j<=iInfo->degree ; j++){
 135                                 tmpcoefficient->data[j] = coefficient[j];
 136                                 //fprintf(stdout,"%2ld:%+10.3e ",j,tmpcoefficient->data[j]);
 137                         }
 138                         //fprintf(stdout,"E=%10.3e Emin=%10.3e \n",E,iInfo->Emin);
 139                         iInfo->Emin = E;
 140                 }
 141         }
 142 /*
 143         fprintf(stderr,"result: ");
 144         for(j=0 ; j<=iInfo->degree ; j++){
 145                 fprintf(stderr,"%2ld: %+10.3e ",j,tmpcoefficient->data[j]);
 146         }
 147         fprintf(stderr,"Emin= %10.3e\n",iInfo->Emin);
 148 */
 149 }
 150
 151 /*mode 0  above & below
 152           -1  below
 153            1  above*/
 154 void
 155 TangentLine(floatVector* spacing, floatVector* scatter, float** belowCoefficient, float** aboveCoefficient, ctfInfo* ini, lctfDetermineInfo* linfo, int mode)
 156 {
 157         long i,j;
 158         floatVector* tmpcoefficient;
 159         float E;
 160         float Pmax = 0.0;
 161         float tmp;
 162         float tmpEmin;
 163         long CutLowSize = 0;
 164         long CutHighSize = 0;
 165         floatVector* Plog;
 166         stepInfo sInfo[linfo->degree+1];
 167         iterationInfo iInfo;
 168
 169         Plog = floatVectorInit(NULL,spacing->size);
 170         iInfo.CutLine = 0.0;
 171         iInfo.degree = linfo->degree;
 172
 173         for(i=0;i<=linfo->degree;i++){
 174                 sInfo[i].step = linfo->step;
 175                 sInfo[i].min = linfo->min;
 176                 sInfo[i].max = linfo->max;
 177         }
 178         __stepNumberCulc(sInfo,linfo->degree);
 179
 180         iInfo.CutLowSize = 0;
 181         iInfo.CutHighSize =0;
 182         iInfo.Emax = 0.0;
 183
 184         tmpcoefficient = floatVectorInit(NULL,linfo->degree);
 185
 186         for(i=0;i<spacing->size;i++){
 187                 Plog->data[i] = log10(scatter->data[i]);
 188
 189                 if(Pmax < Plog->data[i]){
 190                         Pmax = Plog->data[i];
 191                 }
 192
 193                 if(spacing->data[i] < ini->CutLow){
 194                         iInfo.CutLowSize = i+1;
 195                 }
 196                 if(spacing->data[i] < ini->CutHigh){
 197                         iInfo.CutHighSize = i;
 198                 }
 199
 200                 if((spacing->data[i] > ini->CutLow) && (spacing->data[i] < ini->CutHigh)){
 201                         iInfo.Emax += Plog->data[i];
 202                 }
 203         }
 204         iInfo.Emin = iInfo.Emax;
 205         tmpEmin = iInfo.Emin;
 206
 207
 208 /*below*/
 209         if(0 == mode || -1 == mode){
 210                 fprintf(stdout,"below\n");
 211                 iInfo.degree = 0;
 212                 sInfo[0].step = 0.01;
 213                 __stepNumberCulc(sInfo,linfo->degree);
 214                 __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 215                 belowCoefficient[0][0] = tmpcoefficient->data[0];
 216                 sInfo[0].step = linfo->step;
 217                 iInfo.CutLine = tmpcoefficient->data[0] -0.1;
 218                 fprintf(stdout,"%2ld: %10.3e  \n",0,tmpcoefficient->data[0]);
 219
 220                 if(linfo->degree > 0){
 221                         iInfo.degree = 1;
 222                         __stepNumberCulc(sInfo,linfo->degree);
 223                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 224                         for(j=0 ; j<=iInfo.degree ; j++){
 225                                 belowCoefficient[1][j] = tmpcoefficient->data[j];
 226                                 fprintf(stdout,"%2ld: %10.3e  ",j,tmpcoefficient->data[j]);
 227                         }
 228                         fprintf(stdout,"\n");
 229                 }
 230                 for(i=2 ; i<linfo->degree+1 ; i++){
 231                         iInfo.degree = i;
 232                         sInfo[i-1].step = sInfo[i-1].step * 100;
 233                         sInfo[i].step = sInfo[i].step * 100;
 234                         __stepNumberCulc(sInfo,linfo->degree);
 235                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 236
 237                         sInfo[i-2].min = tmpcoefficient->data[i-2];
 238                         sInfo[i-2].max = tmpcoefficient->data[i-2];
 239                         sInfo[i-1].step = sInfo[i-1].step/100;
 240                         sInfo[i].step = sInfo[i].step/100;
 241
 242                         sInfo[i-1].min = tmpcoefficient->data[i-1] - 100*sInfo[i-1].step;
 243                         sInfo[i-1].max = tmpcoefficient->data[i-1] + 100*sInfo[i-1].step;
 244                         sInfo[i].min = tmpcoefficient->data[i] - 100*sInfo[i].step;
 245                         sInfo[i].max = tmpcoefficient->data[i] + 100*sInfo[i].step;
 246
 247                         __stepNumberCulc(sInfo,linfo->degree);
 248                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 249                         sInfo[i-1].min = linfo->min;
 250                         sInfo[i-1].max = linfo->max;
 251                         sInfo[i].min = linfo->min;
 252                         sInfo[i].max = linfo->max;
 253                         for(j=0 ; j<=i ; j++){
 254                                  belowCoefficient[i][j] = tmpcoefficient->data[j];
 255                                  fprintf(stdout,"%2ld: %10.3e  ",j,tmpcoefficient->data[j]);
 256                         }
 257                         fprintf(stdout,"\n");
 258                 }
 259                 iInfo.CutLine = 0.0;
 260         }
 261 /*below end*/
 262 /*above initial*/
 263         if(0 == mode || 1 == mode){
 264                 for(i=0 ; i<=linfo->degree ; i++){
 265                         sInfo[i].step =linfo->step * 10;
 266                         sInfo[i].min = linfo->min * 10;
 267                         sInfo[i].max = linfo->max * 10;
 268                 }
 269                 tmpcoefficient = floatVectorInit(NULL,linfo->degree+1);
 270                 iInfo.Emin = tmpEmin;
 271
 272         /*subtract noise*/
 273         /*      for(i=iInfo.CutLowSize ; i<=iInfo.CutHighSize ; i++){
 274                         fprintf(stdout,"%10.3f %10.3f",spacing->data[i],Plog->data[i]);
 275                         for(j=0 ; j<=linfo->degree ; j++){
 276                                 Plog->data[i] -= belowCoefficient[linfo->degree][j]*__power(spacing->data[i] ,j);
 277                         }
 278                         fprintf(stdout,"%10.3f\n",Plog->data[i]);
 279                 }
 280         */
 281         /*above start*/
 282                 /*above*/
 283                 fprintf(stdout,"above\n");
 284                 iInfo.degree = 0;
 285                 __stepNumberCulc(sInfo,linfo->degree);
 286                 __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,-1.0);
 287                 aboveCoefficient[0][0] = tmpcoefficient->data[0];
 288                 fprintf(stdout,"%2ld: %10.3e  \n",0,tmpcoefficient->data[0]);
 289
 290                 if(linfo->degree > 0){
 291                         iInfo.degree = 1;
 292                         __stepNumberCulc(sInfo,linfo->degree);
 293                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,-1.0);
 294                         aboveCoefficient[1][0] = tmpcoefficient->data[0];
 295                         aboveCoefficient[1][1] = tmpcoefficient->data[1];
 296                         for(i=0 ; i<iInfo.degree+1 ; i++){
 297                                 fprintf(stdout,"%2ld: %10.3e  ",i,tmpcoefficient->data[i]);
 298                         }
 299                         fprintf(stdout,"\n");
 300                 }
 301
 302                 for(j= 2 ; j<=linfo->degree ; j++){
 303                         iInfo.degree = j;
 304                         sInfo[j].step = sInfo[j].step * 100;
 305                         sInfo[j-1].step = sInfo[j-1].step * 100;
 306                         __stepNumberCulc(sInfo,linfo->degree);
 307                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,-1.0);
 308
 309                         sInfo[j-2].step = sInfo[j-2].step / 10;
 310                         sInfo[j-2].min = tmpcoefficient->data[j-2] - 1000*sInfo[j-2].step;
 311                         sInfo[j-2].max = tmpcoefficient->data[j-2] + 1000*sInfo[j-2].step;
 312
 313                         sInfo[j-1].step = sInfo[j-1].step /10;
 314                         sInfo[j-1].min = tmpcoefficient->data[j-1] - 10*sInfo[j-1].step;
 315                         sInfo[j-1].max = tmpcoefficient->data[j-1] + 10*sInfo[j-1].step;
 316
 317                         sInfo[j].step = sInfo[j].step /10;
 318                         sInfo[j].min = tmpcoefficient->data[j] - 10*sInfo[j].step;
 319                         sInfo[j].max = tmpcoefficient->data[j] + 10*sInfo[j].step;
 320
 321                         __stepNumberCulc(sInfo,linfo->degree);
 322                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,-1.0);
 323                         sInfo[j-2].min = tmpcoefficient->data[j-2];
 324                         sInfo[j-2].max = tmpcoefficient->data[j-2];
 325
 326                         sInfo[j-1].step = linfo->step;
 327                         sInfo[j-1].min = tmpcoefficient->data[j-1] - 100*sInfo[j-1].step;
 328                         sInfo[j-1].max = tmpcoefficient->data[j-1] + 100*sInfo[j-1].step;
 329
 330                         sInfo[j].step = linfo->step;
 331                         sInfo[j].min = tmpcoefficient->data[j] - 100*sInfo[j].step;
 332                         sInfo[j].max = tmpcoefficient->data[j] + 100*sInfo[j].step;
 333
 334                         __stepNumberCulc(sInfo,linfo->degree);
 335                         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,-1.0);
 336
 337                         sInfo[j-1].min = linfo->min;
 338                         sInfo[j-1].max = linfo->max;
 339                         sInfo[j-1].step = linfo->step;
 340                         sInfo[j].min = linfo->min;
 341                         sInfo[j].max = linfo->max;
 342                         sInfo[j].step = linfo->step;
 343
 344                         for(i=0 ; i<=j ; i++){
 345                                 aboveCoefficient[j][i] = tmpcoefficient->data[i];
 346                                 fprintf(stdout,"%2ld: %10.3e  ",i,tmpcoefficient->data[i]);
 347                         }
 348                         fprintf(stdout,"\n");
 349                 }
 350         }
 351 /*above end*/
 352 /*      for(i=0 ; i<=linfo->degree ; i++){
 353                 sInfo[i].step =linfo->step;
 354                 sInfo[i].min = linfo->min;
 355                 sInfo[i].max = linfo->max;
 356         }
 357         tmpcoefficient = floatVectorInit(NULL,linfo->degree+1);
 358         iInfo.Emin = tmpEmin;
 359 */
 360 /*below*/
 361 /*      fprintf(stdout,"below\n");
 362         iInfo.degree = 0;
 363         sInfo[0].step = 0.01;
 364         __stepNumberCulc(sInfo,linfo->degree);
 365         __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 366         belowCoefficient[0][0] = tmpcoefficient->data[0];
 367         sInfo[0].step = linfo->step;
 368         iInfo.CutLine = tmpcoefficient->data[0] -0.1;
 369         fprintf(stdout,"%2ld: %10.3e  \n",0,tmpcoefficient->data[0]);
 370
 371         if(linfo->degree > 0){
 372                 iInfo.degree = 1;
 373                 __stepNumberCulc(sInfo,linfo->degree);
 374                 __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 375                 for(j=0 ; j<=iInfo.degree ; j++){
 376                         belowCoefficient[1][j] = tmpcoefficient->data[j];
 377                         fprintf(stdout,"%2ld: %10.3e  ",j,tmpcoefficient->data[j]);
 378                 }
 379                 fprintf(stdout,"\n");
 380         }
 381         for(i=2 ; i<linfo->degree+1 ; i++){
 382                 iInfo.degree = i;
 383                 sInfo[i-1].step = sInfo[i-1].step * 100;
 384                 sInfo[i].step = sInfo[i].step * 100;
 385                 __stepNumberCulc(sInfo,linfo->degree);
 386                 __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 387
 388                 sInfo[i-2].min = tmpcoefficient->data[i-2];
 389                 sInfo[i-2].max = tmpcoefficient->data[i-2];
 390                 sInfo[i-1].step = sInfo[i-1].step/100;
 391                 sInfo[i].step = sInfo[i].step/100;
 392
 393                 sInfo[i-1].min = tmpcoefficient->data[i-1] - 100*sInfo[i-1].step;
 394                 sInfo[i-1].max = tmpcoefficient->data[i-1] + 100*sInfo[i-1].step;
 395                 sInfo[i].min = tmpcoefficient->data[i] - 100*sInfo[i].step;
 396                 sInfo[i].max = tmpcoefficient->data[i] + 100*sInfo[i].step;
 397
 398                 __stepNumberCulc(sInfo,linfo->degree);
 399                 __CoefficientDetermine(spacing,Plog,tmpcoefficient,sInfo,&iInfo,1.0);
 400                 sInfo[i-1].min = linfo->min;
 401                 sInfo[i-1].max = linfo->max;
 402                 sInfo[i].min = linfo->min;
 403                 sInfo[i].max = linfo->max;
 404                 for(j=0 ; j<=i ; j++){
 405                          belowCoefficient[i][j] = tmpcoefficient->data[j];
 406                          fprintf(stdout,"%2ld: %10.3e  ",j,tmpcoefficient->data[j]);
 407                 }
 408                 fprintf(stdout,"\n");
 409         }
 410 */
 411 /*
 412         fprintf(stdout,"degree  : ");
 413         for(i=0 ; i<=linfo->degree ; i++){
 414                 fprintf(stdout,"%8ld",i);
 415         }
 416         fprintf(stdout,"\n");
 417
 418         fprintf(stdout,"above\n");
 419         for(i=0 ; i<=linfo->degree ; i++){
 420                 fprintf(stdout,"%6ld   ",i);
 421                 for(j=0 ; j<=i ; j++){
 422                         fprintf(stdout,"%+7.3f ",aboveCoefficient[i][j]);
 423                 }
 424                 fprintf(stdout,"\n");
 425         }
 426
 427         fprintf(stdout,"below\n");
 428         for(i=0 ; i<=linfo->degree ; i++){
 429                 fprintf(stdout,"%6ld   ",i);
 430                 for(j=0 ; j<=i ; j++){
 431                         fprintf(stdout,"%+7.3f ",belowCoefficient[i][j]);
 432                 }
 433                 fprintf(stdout,"\n");
 434         }
 435         fprintf(stdout,"\n");
 436 */
 437 }