home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e002 / 1.ddi / FIT.CN_ / FIT.CN
Encoding:
Text File  |  1993-12-06  |  14.3 KB  |  569 lines
  1. /**********************************************
  2.  * 
  3.  * 
  4. we create a new menu, called 'Fit', that appears
  5. after 'data' in the menu bar
  6.  * 
  7.  * 
  8. **********************************************/
  9.  
  10. NumPeaks = 1;/* default to single peak for peak related functions */
  11.  
  12. menu -plot /* make sure plot window menu */
  13.  
  14. menu 7 F&it
  15.  
  16. /* this macro is used below as well as
  17.  * for the apparent linear fit menu
  18.  */
  19. def ShowLinearFit {
  20.     /* no need to have "" except if there is '=' involved */
  21.     type "Param\tValue\tsd";
  22.     type "A\t$(LR.A)\t$(LR.sdA)";
  23.     type "B\t$(LR.B)\t$(LR.sdB)";
  24.     type "R  = $(LR.r)";
  25.     type "SD = $(LR.sd), N = $(LR.N)";
  26.     df=LR.N-2;TINY=1e-20;
  27.     t = LR.r*sqrt(df/((1-LR.r+TINY)*(1+LR.r+TINY)));
  28.     P=incbeta(df/(df+t*t),0.5*df,0.5);
  29.     type "P = $(P)";    
  30. }
  31.  
  32. /* this macro is used in initializing various
  33.  * fitting functions
  34.  * The parameters in limit.??? is assumed
  35.  * set x0, y0 accordingly
  36.  */
  37.  
  38. def SetXYoffset {
  39.     minoffset = 0.05;
  40.         /* consider zero offset if 5% from zero */
  41.     x0 = limit.xmin;
  42.     if(abs(x0/(limit.xmax-x0)) < minoffset) {
  43.         x0 = 0
  44.     };
  45.     
  46.     y0 = limit.ymin;
  47.     if(abs(y0/(limit.ymax-y0)) < minoffset) {
  48.         y0 = 0
  49.     };
  50. }
  51.  
  52. def TypeBasicParam { /* type out Chisqr and x0,y0 */
  53.     ty -a "Chisqr = $(Chisqr)";
  54.     type "X offset (x0) = $(x0)";
  55.     type "Y offset (y0) = $(y0)";
  56.     Separator 3;
  57. }
  58.  
  59. def FitGuess {/* a rough fit to get init guess */
  60.     NoFitCurve = -1;/* do not generate fit curve */
  61.     nlsf -H %C;
  62.     nlsf -1;
  63.     nlsf -2;
  64.     nlsf -5;
  65.     nlsf -t;/* terminate fit */
  66. }
  67.     
  68. def KineticInit {
  69.     limit %C;
  70.             /* find max/min in both X and Y for the data */
  71.             
  72.     SetXYoffset;
  73.         /* set y0 and x0 based on limit.###
  74.          * will set y0 or x0 to zero if
  75.          * they are close to zero
  76.          */
  77.          
  78.     /* next we need to find the time constant(t1) and amplitude (A1) */
  79.     del dummy;/* incase it was not deleted from other parts */
  80.  
  81.     copy -s 20 %C dummy;
  82.         /* we interpolate the data with 20 points and put into
  83.          * a temporary data set called dummy
  84.          */
  85.     
  86.     set dummy -f 0;/* this remove any possible offset in x */
  87.     dummy-y0;    /* also remove y offset */
  88.     dummy = ln(dummy);
  89.     LR -B dummy;
  90.         /* linear regression on the beginning */
  91.     del dummy;/* we are interested only in the parameters */
  92.     
  93.     t1=-1/LR.B;t1=$(t1,*3);
  94.     A1=exp(LR.A);A1=$(A1,*3);/* limit them to 3 sig digits */
  95.     if(y0==0) (pv1=0);/* hold y0 if it is zero */
  96.     if(x0==0) (pv2=0);/* hold x0 if it is zero */
  97.     A2=0;A3=0;/* to disable the other two exponentials */
  98. }
  99.  
  100. menu (&Linear Regression) {
  101.     %W=%[%C,'_'];%B=%[%C,>'_'];
  102.     del %W.%B_LINE;/* resulting curve of LR is in wksName.colName_LINE */
  103.     LR %C;        /* linear Regression */
  104.     type -a;    /* open the script window is it is not */
  105.     type Linear Regression for %C:;
  106.     type "Y = A + B * X";
  107.     ShowLinearFit;
  108. #Linear fit to %C
  109. }
  110.  
  111.  
  112. def CheckPRName {
  113.     if(%[%C,>"Polyfit"] != "") {/* > for beyond */
  114.         ty -N "Are you sure you want to do
  115. polynomial regression on %C.
  116. which is already a polynomial fit?";
  117.     };
  118. }
  119.  
  120. /* The next menu is for polynomial fit
  121.  * we will initialize some parameters
  122.  * here, so they will not be init as
  123.  * zero
  124.  */
  125. menu (&Polynomial Regression..) {
  126.     CheckPRName;/* to make sure we are not doing a PR on a PR fit curve */
  127.     CheckVar order 2;/* see if order is defined, if not set it to 2 */
  128.     CheckVar PR.Fit 1; /* 1 = create fit curve, 0 = no fit curve */
  129.     CheckVar PR.Formula 0;
  130.  
  131.     limit %C;fitx1=limit.xmin;fitx2=limit.xmax;
  132.  
  133.     GetNumber [Order (1-9, 1 = linear)] order
  134.     [Show Formula on Plot?] PR.Formula:2
  135.         /* :2 for binery check box instead of edittext */
  136.     [Create Fit Curve?] PR.Fit:2s
  137.         /* s for show items below if yes */
  138.     [Fit curve # pts] fitnpts
  139.     [Fit curve Xmin] fitx1
  140.     [Fit curve Xmax] fitx2
  141.     [Polynomial Fit to %C];
  142.     /* user can cancel, otherwise we will continue */
  143.     
  144.     if(order<=0 || order>9 || fitnpts < 2 || fitx1==fitx2) {
  145.         type -b Parameters invalid; break 1;
  146.     };
  147.     
  148.     /* if no more cancel button, any
  149.      * error must be cause by
  150.      * the DLL routine
  151.      */
  152.     def DLLError {
  153.         type -b Polynomial fit failed! Please check parameters and try again.
  154.     };
  155.  
  156.     PR.A=data(0,0,30);/* dataset to hold parameters */
  157.     dll orgmath PR $(order) PR.A xof(%C) %C;
  158.     /* Polynomial Regression is in the ORGMATH.DLL */
  159.  
  160.     layer -c %C;/* Position of %C in the layer, result in count */
  161.  
  162.     if(PR.Fit){
  163.         GetCuvName PR$(order);
  164.         /* %B is the name of the fit curve */;
  165.         
  166.         (%B)=data(0,0,fitnpts);
  167.         /* (%B) instead of %B, otherwise just the content of %B
  168.          * will be set instead of the name in %B be assigned 
  169.          */
  170.         SetDataRange %B;/* use fitx1, fitx2 and fitnpts */
  171.         /* the fitted curve is now ready */
  172.  
  173.         (%B)=poly(xof(%B),PR.A); /* get the y values for the fit */
  174.  
  175.         graph %B;/* show it in the current layer */
  176.         set %B -c 2;/* color red */
  177.     };
  178.  
  179.     Para2ScriptWindow PR.A order;
  180.     /* Output to the script window is always enabled */
  181.     if(PR.Formula){
  182.         %Z="$(PR.A[1])";/* string could be very long, use Z */
  183.         for(i=2;i<=order+1;i+=1) {
  184.             if(PR.A[i]<0)
  185.                 {%K="$(PR.A[i])"}
  186.             else
  187.                 {%K="+$(PR.A[i])"};
  188.  
  189.             if(i>2)
  190.                 {%Z="%Z%K x\+($(i-1))"}
  191.             else
  192.                 {%Z="%Z%K x"};
  193.         };
  194.         label -s -n PR.F.$(count) "y=%Z";
  195.         /* -s = string substitution, -n = name */
  196.     };
  197. #Polynomial fit to %C
  198. }
  199.  
  200. /* 
  201.  * %1 = coeff data set
  202.  * %2 = order
  203.  */
  204. def Para2ScriptWindow {
  205.     type -a "Polynomial Regression on %C";
  206.     type "y = A0 + A1 x + A2 x^2 + A3 x^3 +...";
  207.  
  208.     Separator 4;/* macro defined in ORGSYS.CNF */
  209.     type -o Tab 55;/* set tab to 10 chars */
  210.     ty Parameter\tValue\tsd;
  211.     
  212.     Separator 4;
  213.  
  214.     for(i=1;i<=%2+1;i+1) {
  215.         ty "A$(i-1)\t$(%1[i],*8)\t$(%1[i+10])";
  216.     };/* The $() notation with *8 means 8 significant digits */
  217.  
  218.     Separator 4;
  219.     
  220.     ty "R  =$(%1[21])";
  221.     ty "R^2=$(%1[23])";
  222.     ty "SD =$(%1[22])";
  223. }
  224.  
  225. Menu
  226.  
  227. Menu "Apparent Linear Fit" {
  228.     del %C.LINE;/* resulting curve of LR is in %C.LINE */
  229.     LR -S %C;        /* linear Regression on linearized scales */
  230.     type -a;    /* open the script window is it is not */
  231.     type Linear Fit for %C on linearized scales;
  232.     type "yscale(Y) = A + B * xscale(X)";
  233.     type where scale() is the current axis scale function;
  234.     ShowLinearFit;
  235. #Linear fit to %C on linearized X Y scales.
  236. }
  237.  
  238. Menu
  239.  
  240. Menu "1 &Exponential Decay" {
  241.     nlsf -f ExpDecay 4;
  242.             /* ExpDecay function can be triple exponential
  243.              * so we restrict it to single by
  244.              * using only the first 4 parameters.
  245.              * params:y0 x0 A1 t1. Must define function first
  246.              */
  247.  
  248.     KineticInit;
  249.     for(i=5;i<=8;i+=1) {/* turn off 2nd and 3rd exp */
  250.         p$(i)=0;/* the parameters */
  251.         pv$(i)=0;/* the vary flag */
  252.     };
  253.     lsqfit;/* the generic fitting macro */
  254.     queue {
  255.         ty -a Fit y0+Ae^(-(x-x0)/t) to %C;
  256.         TypeBasicParam;
  257.         ty "t=$(t1)\tA=$(A1)";
  258.     };#Fit y0+Ae^(-x/t) to %C
  259. }
  260.  
  261. Menu "2 &Exponential Decay" {
  262.     nlsf -f ExpDecay 4;
  263.             /* ExpDecay function can be triple exponential
  264.              * so we restrict it to single by
  265.              * using only the first 4 parameters.
  266.              * params:y0 x0 A1 t1 A2 t2. Must define function first
  267.              */
  268.  
  269.     KineticInit;
  270.     FitGuess;/* the generic fitting macro for init guess */
  271.     t1 = prec(t1,1);/* just the order of magnitude */
  272.     A2=0.5*A1;t2=10*t1;
  273.     GetNum 
  274.         [First Decay Time(t1)] t1
  275.         [Second Decay Time(t2)] t2
  276.         [Initial Time] x0
  277.         [Init Time vary?] pv2:2
  278.         [Y Offset] y0
  279.         [Y Offset Vary?] pv1:2
  280.         [Verify initial guesses];
  281.     pv5=1;pv6=1;/* turn on A2 and t2 variables */
  282.     lsqfit;/* the generic fitting macro */
  283.     queue {
  284.         ty -a Fit y0 + A1e^(-x/t1) + A2e^(-x/t2) to %C;
  285.         TypeBasicParam;
  286.         ty "t1=$(t1,*4)\tA1=$(A1,*4)";
  287.         ty "t2=$(t2,*4)\tA2=$(A2,*4)";
  288.     };#Fit y0 + A1e^(-x/t1) + A2e^(-x/t2) to %C
  289. }
  290.  
  291. Menu "3 &Exponential Decay" {
  292.     nlsf -f ExpDecay 4;
  293.             /* ExpDecay function can be triple exponential
  294.              * so we restrict it to single by
  295.              * using only the first 4 parameters.
  296.              * params:y0 x0 A1 t1 A2 t2. Must define function first
  297.              */
  298.  
  299.     KineticInit;
  300.     FitGuess;/* the generic fitting macro for init guess */
  301.     t1 = prec(t1,1);/* just the order of magnitude */
  302.     A2=0.5*A1;t2=10*t1;A3=0.5*A2;t3=10*t2;
  303.     GetNum 
  304.         [First Decay Time(t1)] t1
  305.         [Second Decay Time(t2)] t2
  306.         [Third Decay Time(t3)] t3
  307.         [Initial Time] x0
  308.         [Init Time vary?] pv2:2
  309.         [Y Offset] y0
  310.         [Y Offset Vary?] pv1:2
  311.         [Verify initial guesses];
  312.  
  313.     pv5=1;pv6=1;/* turn on A2 and t2 variables */
  314.     pv7=1;pv8=1;/* turn on A3 and t3 variables */
  315.     lsqfit;/* the generic fitting macro */
  316.         
  317.     queue {
  318.         ty -a Fit y0 + A1e^(-x/t1) + A2e^(-x/t2) + A3e^(-x/t3) to %C;
  319.         TypeBasicParam;
  320.         ty "t1=$(t1,*4)\tA1=$(A1,*4)";
  321.         ty "t2=$(t2,*4)\tA2=$(A2,*4)";
  322.         ty "t3=$(t3,*4)\tA3=$(A3,*4)";
  323.     };#Fit y0 + A1e^(-x/t1) + A2e^(-x/t2) + A3e^(-x/t3) to %C
  324. }
  325.  
  326. Menu
  327.  
  328. Menu "1 &Exponential Growth" {
  329.     nlsf -f ExpGrow 4;
  330.             /* ExpGrow function can be double exponential
  331.              * so we restrict it to single by
  332.              * using only the first 4 parameters.
  333.              * params:y0 x0 A1 t1. Must define function first
  334.              */
  335.  
  336.     KineticInit;
  337.     t1 = -t1;/* see above, it was for decay */
  338.     
  339.     lsqfit;/* the generic fitting macro */
  340.     queue {
  341.         ty -a Fit y0+Ae^(x/t) to %C;
  342.         ty "y0=$(y0), A=$(A1), t=$(t1)";
  343.     };#Fit y0+Ae^(-x/t) to %C
  344. }
  345.  
  346. /* The Fit1Peak macro is used for both Gaussian and
  347.  * Lorentzian fit. The macro does not need to be defined
  348.  * before the menu, but it is clearer to put them
  349.  * together.
  350.  */
  351.  
  352. def Fit1Peak {
  353.     NumPeaks = 1;/* Gauss,Lorentz functions are multi-peaks functions */
  354.     
  355.     integ %C;nlsf -func %1;
  356.         /* %1 is the argument passed to the macro */
  357.  
  358.     if(integ.dx < 0) (integ.dx *=-1;integ.area *= -1;);
  359.     /* if x is in descending order, we need to turn these around */
  360.  
  361.     limit %C;
  362.  
  363.     SetXYoffset;
  364.         /* see Kinetic above */
  365. /* paramters
  366.  * y0 = P1
  367.  * xc = P2
  368.  * w  = P3
  369.  * A  = P4
  370.  */
  371.     if(y0==0) (pv1=0);/* y0 is the 1st parameter, hold it if zero */
  372.     
  373.     P2=integ.x0;
  374.     P4=integ.Area;
  375.     P3=integ.dx;
  376.         /* parameter initialization from the
  377.          * integration result(integ.x0,integ.area,integ.dx).
  378.          */
  379.     lsqfit;
  380.         /* the generic macro for iterating 10 times */
  381.     queue {
  382.         x = P4/(P3*sqrt(pi/2));/* peak height */
  383.         ty -a %1 fit to %C;
  384.         ty "Area\tCenter\tWidth\tOffset\tHeight";
  385.         ty "$(P4,*5)\t$(P2,*5)\t$(P3,*5)\t$(P1,*5)\t$(x,*5)"
  386.         /* *5 if a formatting statement for 5 significant digits */
  387.     };
  388.     /* using a queue command will 
  389.      * ensure proper windows updating.
  390.      */
  391. }
  392.  
  393. Menu (1 &Gaussian) {
  394.     Fit1Peak Gaussian;#Gaussian fit to %C
  395. }
  396.  
  397. Menu (1 &Lorentzian) {
  398.     Fit1Peak Lorentz;#Lorentzian fit to %C
  399. }
  400.  
  401. Menu Sig&moidal {
  402.     axis -pget X Scale ScaleType;
  403.     if(ScaleType==1) (%S = Logistic) else (%S = Boltzman);
  404.     nlsf -f %S;
  405.     lsqfit;
  406.         /* the generic macro for iterating 12 times */
  407.     queue {
  408.         ty -a Sigmoidal(%S) fit to %C;
  409.         ty "   Chisqr = $(Chisqr)";
  410.         Separator 3;
  411.         ty " Init(A1)  = $(A1,*5)\t$(PE1,*3)";
  412.         ty "Final(A2)  = $(A2,*5)\t$(PE2,*3)";
  413.         ty "XatY50(x0) = $(x0,*5)\t$(PE3,*3)";
  414.         if(ScaleType==1) {/* log scale, must be logistic */
  415.             ty "Power(p)  = $(p,*5)\t$(PE4,*3)";
  416.             Separator 3;
  417.             ty "   XatY20 = $( x0*(0.2/0.8)^(1/p) )";
  418.             ty "   XatY80 = $( x0*(0.8/0.2)^(1/p) )";
  419.         }
  420.         else {
  421.             ty "Width(dx) = $(dx,*5)\t$(PE4,*3)";
  422.             Separator 3;
  423.             ty "   XatY20 = $( x0+dx*ln(0.2/0.8)  )";
  424.             ty "   XatY80 = $( x0+dx*ln(0.8/0.2)  )";
  425.         };
  426.         /* *5 if a formatting statement for 5 significant digits */
  427.     };
  428.     /* using a queue command will 
  429.      * ensure proper windows updating.
  430.      */
  431.     
  432. }
  433.  
  434. /* assume fit to %C and all fitting parameters are in place */
  435. def MakePeaks {
  436.     %B=xof(%C);
  437.     if(%B=="")
  438.         return 1;/* no x column, can not work */
  439.     nlsf.funcx$=%B;
  440.  
  441.     %W=%[%C,'_'];/* get wks name */
  442.     /* del fit and all peak columns */
  443.     del %W_fit;
  444.     for(i=1;i<=12;i+=1) (del %W_peak$(i)) /* max 12 peaks */
  445.  
  446.     /* put columns back to wks */
  447.     win -o %W {
  448.         work -v fit;
  449.         for(i=1;i<=NumPeaks;i+=1) (work -v peak$(i));
  450.     };
  451.     break -b Generating curves...;
  452.     break -r 0 (numPeaks+1);
  453.     nlsf.funcCol$=%W_Fit;
  454.     nlsf.make(func);
  455.  
  456.     break -p 1;
  457.     layer -i %W_Fit;
  458.  
  459.     set %W_Fit -c color(red);
  460.  
  461.     for(i=1;i<=NumPeaks;i+=1) {
  462.         break -p (i+1);
  463.         nlsf -p Save;/* save copy of parameters */
  464.         nlsf -y 0 3 3;/* value=0, init = 3(0,1,2), inc = 3 */
  465.         y = A$(i);
  466.         p$(1+3*i)= y;
  467.  
  468.         nlsf.funcCol$=%W_peak$(i);
  469.         nlsf.make(func);
  470.  
  471.         nlsf -p Restore;
  472.         layer -i %W_peak$(i);
  473.         set %W_peak$(i) -c 4;/* blue */
  474.         set %W_peak$(i) -d 2;/* dotted */
  475.     };
  476.     break -e;
  477. }
  478.  
  479. def FitNPeaks {
  480.     GetNumber -s [Number of Peaks] NumPeaks;
  481.     if(NumPeaks > 12) {ty -b Too many peaks;break 1};
  482.  
  483.     del _Peak_Y;
  484.     create _Peak_Y -X 20;
  485.     
  486.     integ %C;nlsf -func %1;
  487.         /* %1 is the argument passed to the macro */
  488.  
  489.     if(integ.dx < 0) (integ.dx *=-1;integ.area *= -1;);
  490.     /* if x is in descending order, we need to turn these around */
  491.  
  492.     def QuitFit [del _Peak_Y];
  493.  
  494.     limit %C;
  495.     SetXYoffset;
  496.     /* integration was carried out from zero, we need
  497.      * to remove the offset */
  498.     integ.area -= y0*(limit.xmax-limit.xmin);
  499.     integ.area/=NumPeaks;
  500.     integ.dx/=NumPeaks;
  501.     GetNumber -s [Initial half width estimate] integ.dx;
  502.     for(i=1;i<=NumPeaks;i+=1) {
  503.         A$(i)=integ.area;
  504.         w$(i)=integ.dx
  505.     };
  506.  
  507.     if(y0==0) (pv1=0);/* y0 is the 1st parameter, hold it if zero */
  508.     def EndToolbox {
  509.         doTool 0;/* reset toolbox */
  510.         for(i=1;i<=NumPeaks;i+=1) { XC$(i)=_Peak_A[i] };
  511.         nlsf -H %C;
  512.         noFitCurve=-1;/* we will generate the fit cureves by script */
  513.         nlsf -v 0 1 3;/* hold all xc's constant */
  514.         step*=2;
  515.         nlsf -2;step/=2;nlsf -2;/* carry out four iterations */
  516.         step/=2;nlsf -2;
  517.         
  518.         nlsf -v 1 1 3;/* release all holds on xc's */
  519.         step*=2;
  520.         nlsf -2;step/=2;
  521.         nlsf -5;
  522.         nlsf -t;/* terminate fit */
  523.         queue {
  524.             ty -a %1($(NumPeaks)) fit to %C;
  525.             ty "Peak\tArea\tCenter\tWidth\tHeight";
  526.             for(i=1;i<=NumPeaks;i+=1) {
  527.                 x = A$(i);y=Xc$(i);t=w$(i);
  528.                 z = x/(t*sqrt(pi/2));/* peak height */
  529.                 ty "$(i)\t$(x,*5)\t$(y,*5)\t$(t,*5)\t$(z,*5)";
  530.             };
  531.             ty "Yoffset = $(y0)";
  532.             MakePeaks;/* macro defined above */
  533.         };
  534.     };
  535.     GetPts -n _Peak_Y NumPeaks {Double-click(or Enter) to set Peak$(1+count)};
  536. }
  537.  
  538. Menu
  539. Menu (Multiple &Gaussian) {
  540.     FitNPeaks Gaussian;#Multiple Gaussian fit to %C
  541. }
  542.  
  543. Menu (Multiple &Lorentzian) {
  544.     FitNPeaks Lorentz;#Multiple Lorentzian fit to %C
  545. }
  546. menu
  547.  
  548. Menu "&Select Fitting Function.." {
  549.     nlsf -f;
  550. #Select one of the built-in \
  551. fitting functions, or enter your own.
  552. }
  553.  
  554. Menu "Start &Fitting Session.." {
  555.     nlsf %C;
  556.     def QuitFit (doTool -wc nlsf);
  557.     doTool -w nlsf;/* load the nlsf toolbar */
  558. #Start the user controlled \
  559. non-linear least square fitting session
  560. }
  561.  
  562. Menu "&Parameters/Simulate.." {
  563.     fitX1=X1;fitX2=X2;
  564.     nlsf -e;
  565. #Open the Fitting Function dialog box
  566. }
  567.  
  568. /* End of Fit menu */
  569.