unit Qanal;
{ ======================================================
  WinWCP - Quantal Analysis Module (c) J. Dempster 1996
  18/2/98 ... Now has All Record/Range radio buttons
  ======================================================}

interface

uses
  SysUtils, WinTypes, WinProcs, Messages, Classes, Graphics, Controls,
  Forms, Dialogs, StdCtrls, Global, Fileio, shared, printers, plotlib ;

type
  TQuantFrm = class(TForm)
    ControlGrp: TGroupBox;
    bDoAnalysis: TButton;
    bClose: TButton;
    GroupBox2: TGroupBox;
    Label5: TLabel;
    cbEvokedType: TComboBox;
    MiinGrp: TGroupBox;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    cbMiniType: TComboBox;
    rbMiniEventsAvailable: TRadioButton;
    rbUserEntered: TRadioButton;
    EdMiniAmplitude: TEdit;
    EdMiniStDev: TEdit;
    GroupBox3: TGroupBox;
    rbCurrent: TRadioButton;
    rbPotentials: TRadioButton;
    edVRest: TEdit;
    Label1: TLabel;
    rbPoisson: TRadioButton;
    rbBinomial: TRadioButton;
    edVRev: TEdit;
    Label6: TLabel;
    edCorrectionFactor: TEdit;
    Label7: TLabel;
    mmResults: TMemo;
    GroupBox8: TGroupBox;
    edRecRange: TEdit;
    rbAllRecords: TRadioButton;
    rbRange: TRadioButton;
    cbChannel: TComboBox;
    Label8: TLabel;
    procedure bDoAnalysisClick(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
    procedure rbUserEnteredClick(Sender: TObject);
    procedure rbMiniEventsAvailableClick(Sender: TObject);
    procedure rbCurrentClick(Sender: TObject);
    procedure rbPotentialsClick(Sender: TObject);
    procedure bCloseClick(Sender: TObject);
  private
    { Private declarations }
    function UseRecord ( const RecH : TRecHeader ; RecType : string ) : Boolean ;
  public
    { Public declarations }
    procedure PrintResults ;
  end;

var
  QuantFrm: TQuantFrm;

implementation
{$R *.DFM}
uses MDIForm ;


procedure TQuantFrm.FormCreate(Sender: TObject);
{ ---------------
  Initialise form
  ---------------}
begin
     { Disable "Quantal Content" item of "Analysis" menu }
     Main.mnQuantalContent.enabled := false ;

     { Type of record containing evoked signals }
     cbEvokedType.items := RecordTypes ;
     cbEvokedType.itemIndex := cbEvokedType.items.IndexOf('TEST') ;
     { Type of record containing spontaneous miniature signals }
     cbMiniType.items := RecordTypes ;
     cbMiniType.itemIndex := cbMiniType.items.IndexOf('MINI') ;

     edRecRange.text := format(' 1-%d',[fH.NumRecords] ) ;

     cbChannel.items := ChannelNames ;
     if cbChannel.ItemIndex < 0 then cbChannel.ItemIndex := 0 ;
     end;


procedure TQuantFrm.bDoAnalysisClick(Sender: TObject);
{---------------------------
 Do quantal content analysis
 ---------------------------}
type
TSignal = record
        Sum : single ;
        Mean : single ;
        SD : single ;
        Variance : single ;
        Num : LongInt ;
        end ;
TQuantalContent = record
                Direct : single ;
                Variance : single ;
                Failures : single ;
                end ;

var
   Rec,RecStart,RecEnd : LongInt ;
   UseCh,i : Integer ;
   Evoked : TSignal ;
   QuantalContent : TQuantalContent ;
   Mini : TSignal ;
   NumFailures : LongInt ;
   x,VRest,VReversal,VDrive,CorrectionFactor : single ;
   ReleaseProbability,PoolSize : single ;
   Units : string ;
   rH : ^TRecHeader ;
begin

     New(rh) ;

     { Determine record range/channels to be plotted }
     if rbAllRecords.Checked then begin
        RecStart := 1 ;
        RecEnd := FH.NumRecords ;
        end
     else begin
        GetIntRangeFromEditBox(edRecRange,RecStart,RecEnd,1,fH.NumRecords) ;
        end ;

     UseCh := cbChannel.ItemIndex ;
     Units := Channel[UseCh].ADCUnits ;
     mmResults.Clear ;
     mmResults.Lines.Add('Quantal Analysis') ;

     { ** Calculate average of peak evoked and mini signals ** }

     Evoked.Sum := 0. ;
     Evoked.Num := 0 ;
     Mini.Sum := 0. ;
     Mini.Num := 0 ;
     NumFailures := 0 ;
     { Calculate driving force if potentials }
     if rbPotentials.checked then begin
        VRest := ExtractFloat( edVRest.text, -90. ) ;
        VReversal := ExtractFloat( edVRev.text, 0. ) ;
        VDrive := VRest - VReversal ;
        CorrectionFactor := ExtractFloat( edCorrectionFactor.text, 1. ) ;
        end ;

     for Rec := RecStart to RecEnd do begin
         { Read record analysis block from file }
         GetRecordHeaderOnly( fH, rH^, Rec ) ;
         { Add evoked peak to sum }
         if UseRecord( rH^, cbEvokedType.text ) then begin
            x := rH^.Analysis.Value[UseCh,vPeak] ;
            { Correct potentials for non-linear summation }
            if rbPotentials.checked then x := x/(1. - Abs((CorrectionFactor*x)/VDrive)) ;
            Evoked.Sum := Evoked.Sum + x ;
            Inc(Evoked.Num) ;
            end ;

         { Add spontaneous peak to sum }
         if UseRecord( rH^, cbMiniType.text ) then begin
            Mini.Sum := Mini.Sum + rH^.Analysis.Value[UseCh,vPeak] ;
            Inc(Mini.Num) ;
            end ;
         { Number of stimulations which failed to evoked a post-synaptic signal }
         if UseRecord( rH^, 'FAIL' ) then begin
            Inc(NumFailures) ;
            Inc(Evoked.Num) ;
            end ;

         end ;
     if Evoked.Num > 0 then Evoked.Mean := Evoked.Sum/Evoked.Num ;
     if Mini.Num > 0 then Mini.Mean := Mini.Sum/Mini.Num ;

     { ** Calculate standard dev. of peak evoked and mini signals ** }

     Evoked.Sum := 0. ;
     Mini.Sum := 0. ;
     for Rec := RecStart to RecEnd do begin
         { Read record analysis block from file }
         GetRecordHeaderOnly( fH, rH^, Rec ) ;
         { Add evoked residual to S.D. summation }
         if UseRecord( rH^, cbEvokedType.text ) then begin
            x := rH^.Analysis.Value[UseCh,vPeak] - Evoked.Mean ;
            Evoked.Sum := Evoked.Sum + x*x ;
            end ;

         { Add spontaneous residual to sum }
         if UseRecord( rH^, cbMiniType.text ) then begin
            x := rH^.Analysis.Value[UseCh,vPeak] - Mini.Mean ;
            Mini.Sum := Mini.Sum + x*x ;
            end ;
         end ;
     if Evoked.Num > 1 then Evoked.SD := Sqrt( Evoked.Sum/(Evoked.Num -1. ) )
                       else Evoked.SD := 0. ;
     Evoked.Variance := Evoked.SD*Evoked.SD ;
     if Mini.Num > 1 then Mini.SD := Sqrt( Mini.Sum/(Mini.Num -1. ) )
                     else Mini.SD := 0. ;
     Mini.Variance := Mini.SD*Mini.SD ;

     { ** Report mean/st.dev of evoked signal amplitudes ** }



     if Evoked.Num > 0 then begin

        if rbPotentials.checked then begin
           mmResults.Lines.Add(format('Evoked potentials (n=%d)',[Evoked.Num])) ;
           mmResults.Lines.Add(
           format('VRest = %.3g %s  VRev. = %.3g %s  Cor. factor (f)= %.3g',
                  [VRest,Units,VReversal,Units,CorrectionFactor] ) ) ;
           end
        else mmResults.Lines.Add(format('Evoked currents (n=%d)',[Evoked.Num])) ;

        mmResults.Lines.add( format( 'Mean = %.3g %s',[Evoked.Mean,Units] )) ;
        if Evoked.Num > 1 then
           mmResults.Lines.add( format( 'Standard deviation = %.3g %s',
                                        [Evoked.SD,Units])) ;
        end ;

       { Report mean/st.dev of miniature signal amplitudes }


     if rbMiniEventsAvailable.checked then begin
        mmResults.Lines.add( ' ' ) ;
        if rbPotentials.checked then
             mmResults.Lines.Add(format('Spontaneous miniature potentials (n=%d)',
                                         [Mini.Num]))
        else mmResults.Lines.Add(format('Spontaneous miniature currents (n=%d)',
                                       [Mini.Num])) ;

        if Mini.Num > 0 then begin
           mmResults.Lines.add( format( 'Mean = %.3g %s',[Mini.Mean,Units] )) ;
           if Mini.Num > 1 then
               mmResults.Lines.add( format( 'Standard deviation = %.3g %s',
                                            [Mini.SD,Units] )) ;
           end
           else mmResults.Lines.add( 'No minis found!' ) ;
        end
     else begin
          mmResults.Lines.add( ' ' ) ;
          if rbPotentials.checked then
               mmResults.Lines.Add('Spontaneous miniature potentials' )
          else mmResults.Lines.Add('Spontaneous miniature currents ') ;

          Mini.Mean := ExtractFloat( edMiniAmplitude.text, 0. ) ;
          Mini.SD := ExtractFloat( edMiniStDev.text, 0. ) ;
          Mini.Num := 1 ;

          if Mini.Mean <> 0. then begin
             mmResults.Lines.add( format( 'Mean = %.3g %s (user entered)',
                                         [Mini.Mean,Units] )) ;
             mmResults.Lines.add( format( 'Standard deviation = %.3g %s (user entered)',
                                         [Mini.SD,Units] )) ;
             Mini.Num := 1 ;
             end
          else mmResults.Lines.add( 'Invalid amplitude!' ) ;
          end ;

     mmResults.Lines.add( ' ' ) ;
     if (Evoked.Num > 0) and (Mini.Num > 0) then begin
        QuantalContent.Direct := Evoked.Mean / Mini.Mean ;
        mmResults.Lines.add( format( 'Quantal content = %.3g (direct method)',
                                     [QuantalContent.Direct] ) );
        end ;

     if (Evoked.Num > 0) and rbPoisson.checked then begin
        { Calculate quantal content by variance method
          (NB only valid for poisson distributions) }

        QuantalContent.Variance := (Evoked.Mean*Evoked.Mean)/Evoked.Variance ;
        mmResults.Lines.add( format( 'Quantal content = %.3g (variance method)',
                                     [QuantalContent.Variance] ) );

        if NumFailures > 0 then  begin
           QuantalContent.Failures := ln( (Evoked.Num+NumFailures)/NumFailures ) ;
           mmResults.Lines.add( format( 'Quantal content = %.3g (failures method)',
                                         [QuantalContent.Failures] ) );
           end ;
        end ;

     if (Evoked.Num > 0) and (Mini.Num > 0) and (not rbPoisson.checked) and
        (Abs(Mini.Mean) > 0. ) then begin
        ReleaseProbability := 1. - Evoked.Variance /
                                   (QuantalContent.Direct*Mini.Mean*Mini.Mean)
                                 + Mini.Variance/(Mini.Mean*Mini.Mean) ;
        PoolSize := QuantalContent.Direct / ReleaseProbability ;
        mmResults.Lines.add( ' ' ) ;
        mmResults.Lines.add( 'Binomial Analysis' ) ;
        mmResults.Lines.add( format( 'Release Probability = %.3g',
                                      [ReleaseProbability]));
        mmResults.Lines.add( format( 'Pool size = %.3g',[PoolSize]));

        end ;

     if (Evoked.Num = 0) and (Mini.Num = 0) then
        mmResults.Lines.add( 'No records available for analysis!' ) ;

     { Copy results to log file }
     for i := 1 to mmResults.Lines.Count do WriteToLogFile( mmResults.Lines[i] ) ;

     Dispose(rh) ;

     Main.Print1.Enabled := True ;
     end;


function TQuantFrm.UseRecord ( const RecH : TRecHeader ;
                                 RecType : string ) : Boolean ;
{ -----------------------------------------------------
  Select record for inclusion in quantal analysis
  -----------------------------------------------------}
begin
     if (RecH.Status = 'ACCEPTED') and RecH.Analysis.Available
        and ( (RecH.RecType = RecType) or ( RecType = 'ALL') ) then
        UseRecord := True
     else UseRecord := False ;
     end ;


procedure TQuantFrm.FormDestroy(Sender: TObject);
begin
     QuantChildExists := False ;
     end;

procedure TQuantFrm.FormClose(Sender: TObject; var Action: TCloseAction);
begin
     { Enable "Quantal Content" item of "Analysis" menu }
     Main.mnQuantalContent.enabled := True ;
     Main.Print1.Enabled := False ;
     Action := caFree ;
     end;

procedure TQuantFrm.rbUserEnteredClick(Sender: TObject);
begin
     edMiniAmplitude.enabled := True ;
     edMiniStDev.enabled := True ;
     cbMiniType.enabled := False ;
     end;

procedure TQuantFrm.rbMiniEventsAvailableClick(Sender: TObject);
begin
     edMiniAmplitude.enabled := False ;
     edMiniAmplitude.text := ' ' ;
     edMiniStDev.enabled := False  ;
     edMiniStDev.text := ' ' ;
     cbMiniType.enabled := True ;
     end;

procedure TQuantFrm.rbCurrentClick(Sender: TObject);
begin
     edVRev.enabled := False ;
     edVRest.enabled := False ;
     edCorrectionFactor.enabled := False ;
     end;

procedure TQuantFrm.rbPotentialsClick(Sender: TObject);
begin
     edVRev.enabled := True ;
     edVRest.enabled := True ;
     edCorrectionFactor.enabled := True ;
     end;

procedure TQuantFrm.bCloseClick(Sender: TObject);
begin
     close ;
     end;

procedure TQuantFrm.PrintResults ;
{ -----------------------
  Print the results table
  -----------------------}
var
   CharWidth,CharHeight,Row : Integer ;
   PageLeft,PageTop,PageBottom,Line : Integer ;
   FontScale : Integer ;
begin

     Screen.Cursor := crHourglass ;

     { Set print font and size }
     Printer.Canvas.font.name := Settings.Plot.FontName ;
     FontScale := PrinterPointsToPixels(10) ;
     Printer.Canvas.font.Height := FontScale ;

     CharWidth := Printer.canvas.TextWidth('X') ;
     CharHeight := Printer.canvas.TextHeight('X') ;
     PageTop := CharHeight*5 ;
     PageBottom := printer.PageHeight - PageTop ;
     PageLeft := CharWidth*8 ;

     Printer.BeginDoc ;

     { Print resultys
       ===========}

     Line := PageTop ;
     printer.canvas.textout(PageLeft,Line, 'File ... ' + fH.FileName ) ;
     Line := Line + CharHeight ;
     printer.canvas.textout(PageLeft,Line, fH.IdentLine) ;
     Line := Line + CharHeight*2 ;

     for Row := 0 to mmResults.Lines.Count-1 do begin
         printer.canvas.textout( PageLeft, Line, mmResults.Lines[Row] ) ;
         Line := Line + CharHeight ;
         end ;

     Printer.EndDoc ;

     Screen.Cursor := crDefault ;

     end ;


end.
