unit Simsyn;
{ ==================================================================
  WinWCP ... Synaptic current simulation module (c) J. Dempster 1996-97
  ==================================================================}

interface

uses
  SysUtils, WinTypes, WinProcs, Messages, Classes, Graphics, Controls,
  Forms, Dialogs, ExtCtrls, StdCtrls, Shared, FileIo, Global, maths ;

type
  TSynapseSim = class(TForm)
    pbDisplay: TPaintBox;
    lbTMin: TLabel;
    lbTMax: TLabel;
    Timer: TTimer;
    REleasegrp: TGroupBox;
    Label6: TLabel;
    lbP: TLabel;
    edPoolSize: TEdit;
    edReleaseProbability: TEdit;
    QuantumGrp: TGroupBox;
    GroupBox1: TGroupBox;
    Label7: TLabel;
    edNumRecords: TEdit;
    GroupBox2: TGroupBox;
    Label8: TLabel;
    rbPotentials: TRadioButton;
    rbCurrents: TRadioButton;
    EdVRest: TEdit;
    bStart: TButton;
    bAbort: TButton;
    bClose: TButton;
    edProgress: TEdit;
    GroupBox3: TGroupBox;
    edTau1: TEdit;
    Label3: TLabel;
    edRatio: TEdit;
    Label5: TLabel;
    edTau2: TEdit;
    Label4: TLabel;
    ck2Exponential: TCheckBox;
    GroupBox4: TGroupBox;
    edQuantumAmplitude: TEdit;
    edQuantumStDev: TEdit;
    edNoiseRMS: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    Label9: TLabel;
    GroupBox5: TGroupBox;
    edTauRise: TEdit;
    Label11: TLabel;
    edLatency: TEdit;
    Label10: TLabel;
    procedure TimerTimer(Sender: TObject);
    procedure bStartClick(Sender: TObject);
    procedure bAbortClick(Sender: TObject);
    procedure FormDeactivate(Sender: TObject);
    procedure bCloseClick(Sender: TObject);
    procedure ck2ExponentialClick(Sender: TObject);
    procedure rbCurrentsClick(Sender: TObject);
    procedure rbPotentialsClick(Sender: TObject);
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
    procedure FormShow(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormActivate(Sender: TObject);
    procedure edQuantumAmplitudeKeyPress(Sender: TObject; var Key: Char);
    procedure edLatencyKeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
    procedure UpdateAmplitudes ;
    procedure HeapBuffers( Operation : THeapBufferOp ) ;
  public
    { Public declarations }
  end;
function binomial( pIn, n : Single ) : Single ;
function gammln( xx : Single ) : Single ;

var
  SynapseSim: TSynapseSim;

implementation
uses
    SSQUnit,MDIForm ;
{$R *.DFM}
type
    TState = (Idle, DoSimulation, ClearDisplay) ;
var
   State : TState ;
   NumRecordsDone,NumRecordsRequired,ChSim : LongInt ;
   FirstRecord : Boolean ;
   RH : ^TRecHeader ; { Record header }
   ADC : ^TIntArray ; { Digitised signal buffer }
   A1,A2,Tau1,Tau2,Ratio,Amplitude,MEPC,NoiseRMS,AmpTemp,n,p : Single  ;
   QuantumAmplitude,QuantumStDev,xStart,ySig : Single ;
   TauRise,Latency : single ;
   Units : string ;
   BuffersAllocated : boolean ;{ Indicates if memory buffers have been allocated }


procedure TSynapseSim.HeapBuffers( Operation : THeapBufferOp ) ;
{ -----------------------------------------------
  Allocate/deallocation dynamic buffers from heap
  -----------------------------------------------}
begin
     case Operation of
          Allocate : begin
             if not BuffersAllocated then begin
                New(ADC) ;
                New(RH) ;
                BuffersAllocated := True ;
                end ;
             end ;
          Deallocate : begin
             if BuffersAllocated then begin
                Dispose(ADC) ;
                Dispose(RH) ;
                BuffersAllocated := False ;
                end ;
             end ;
          end ;
     end ;


procedure TSynapseSim.FormCreate(Sender: TObject);
begin
     { Disable "Synaptic Current" item in "Simulation" menu }
     Main.Synapse.enabled := false ;
     BuffersAllocated := False ;
     end;

procedure TSynapseSim.FormShow(Sender: TObject);
{ --------------------------------------
  Initialisations when form is displayed
  --------------------------------------}
const
     ReleaseProbability = 0.5 ;
     PoolSize = 40 ;
var
   Tau : single ;
   QuantalAmplitude : single ;
begin
     { Create buffers }
     HeapBuffers( Allocate ) ;

     Timer.Enabled := True ;
     State := ClearDisplay ;
     bStart.Enabled := True ;
     bAbort.Enabled := False ;

     edPoolSize.text := format( ' %d ', [Trunc(PoolSize)]) ;
     edReleaseProbability.text := format( ' %.3f ', [ReleaseProbability] ) ;
     rbCurrents.checked := True ;
     Units := 'nA' ;
     UpdateAmplitudes ;
     edRatio.text := ' ' ;
     ck2Exponential.checked := false ;
     edtau2.text := ' ' ;

     end ;


procedure TSynapseSim.TimerTimer(Sender: TObject);
var
   i,j,iStartEPC,xPix,yPix,ch,ChOffset,MaxCh, nOffScreen : LongInt ;
  x,y,dx,NoiseTemp,QuantalContent,Rise,Val,VRest : Single  ;
begin

     case State of
          ClearDisplay : Begin
                       { Erase display }
                       lbTMin.Left := pbDisplay.Left ;
                       lbTMin.Top := pbDisplay.Top + pbDisplay.Height + 2 ;
                       lbTMax.Left := pbDisplay.Left + pbDisplay.Width
                                      - canvas.TextWidth(lbTMax.caption) ;
                       lbTMax.Top := lbTMin.Top ;
                       lbTMin.caption := '' ;
                       lbTMax.caption := '' ;
                       EraseDisplay(pbDisplay) ;
                       State := Idle
                       end ;

          DoSimulation : Begin

              { Initialisations when first record is done }
              if FirstRecord then begin

                 { Set number of channels = 1 (only of file is empty }
                 if RawFH.NumRecords = 0 then begin
                    RawFH.NumChannels := 1 ;
                    Channel[0].ChannelOffset := 0 ;
                    end ;

                 ChSim := 0 ;
                 Channel[ChSim].ADCUnits := Units ;
                 if Units = 'mV' then Channel[ChSim].ADCName := 'Vm'
                                 else Channel[ChSim].ADCName := 'Im' ;

                 WriteToLogFile( 'Synaptic current simulation') ;

                 { Rise time constant of quantal event }
                 TauRise := GetFromEditBox(edTauRise,0.1,0.0,10000.,'%.2f',
                         Settings.TUnits)*Settings.TUnScale ;
                 WriteToLogFile( 'Tau(rise) = ' + edTauRise.text ) ;

                 { Decay time constant of quantal event }
                 Tau1 := GetFromEditBox(edTau1,1.,0.01,10000.,'%.2f',
                         Settings.TUnits)*Settings.TUnScale ;
                 WriteToLogFile( 'Tau(1) = ' + edTau1.text ) ;

                 { Latency of quantal event }
                 Latency := GetFromEditBox(edLatency,0.,0.0,10000.,'%.2f',
                         Settings.TUnits)*Settings.TUnScale ;
                 WriteToLogFile( 'Latency = ' + edLatency.text ) ;

                 { Sampling interval }
                 RawFH.dt := (Tau1*10. ) / RawFH.NumSamples ;
                 A1 := 1. ;
                 if ck2Exponential.checked then begin
                    { Time constant for 2nd decay exponential }
                    Tau2 := GetFromEditBox(edTau2,1.,0.01,10000.,'%.3g',
                            Settings.TUnits)*Settings.TUnScale ;
                    WriteToLogFile( 'Tau(2) = ' + edTau2.text ) ;
                    { Ratio of 1st/2nd exponential amplitudes }
                    Ratio := GetFromEditBox(edRatio,1.,0.1,100.,'%.3g','') ;
                    WriteToLogFile( 'Ratio = ' + edRatio.text ) ;
                    RawFH.dt := MaxFlt([(Tau2*4. )/RawFH.NumSamples,RawFH.dt]) ;
                    A2 := 1. / (1. + Ratio) ;
                    A1 := 1. - A2 ;
                    end ;

                 { Amplitude of quantal event }
                 QuantumAmplitude := GetFromEditBox(edQuantumAmplitude,
                                     1.,-1E4,1E4, '%.1f', Channel[ChSim].ADCUnits ) ;
                 WriteToLogFile( 'Quantum amplitude = ' + edQuantumAmplitude.text ) ;
                 { Standard deviation of quantum amplitude }
                 QuantumStDev := GetFromEditBox(edQuantumStDev,
                                     0.1,0.,1E4, '%.2f', Channel[ChSim].ADCUnits ) ;
                 WriteToLogFile( 'Quantum st. dev. = ' + edQuantumStDev.text ) ;
                 { Size of pool of releasable quanta in nerve terminal }
                 n := GetFromEditBox(edPoolSize, 1., 1., 1E4, '%.0f', '' ) ;
                 WriteToLogFile( 'Pool Size = ' + edPoolSize.text ) ;
                 { Prob. of release for a quantum }
                 p := GetFromEditBox(edReleaseProbability,0.5,0.,1.,'%.2f','') ;
                 WriteToLogFile( 'Release Probability = ' + edReleaseProbability.text ) ;
                 { Background noise RMS }
                 NoiseRMS := GetFromEditBox(edNoiseRMS,
                             0., 0., 1E4, '%.2f', Channel[ChSim].ADCUnits) ;
                 WriteToLogFile( 'Noise RMS = ' + edNoiseRMS.text ) ;
                 { Temp variable for amplitude random number generator }
                 AmpTemp := 1. ;

                 { Set channel parameters }
                 { Estimate maximal quantal content of 100 trials }
                 QuantalContent := 0. ;
                 for i := 1 to 100 do
                     QuantalContent := MaxFlt([QuantalContent,Binomial(p,n)]);

                 { If this is the first record in the file create
                   an appropriate scaling factor }
                 if RawfH.NumRecords = 0 then begin
                    RawfH.NumChannels := 1 ;
                    Channel[ChSim].ADCCalibrationFactor := 1. ;
                    Channel[ChSim].ADCScale := Abs(QuantalContent*QuantumAmplitude*1.2)
                                               / MaxADCValue ;
                    rH^.ADCVoltageRange[ChSim] :=  Channel[ChSim].ADCCalibrationFactor
                                * ( Channel[ChSim].ADCScale * (MaxADCValue+1) ) ;
                    end
                 else begin
                      { Adjust record scaling factor if signal is too small/big for record }
                      Val := QuantalContent*QuantumAmplitude/Channel[ChSim].ADCScale ;
                      if (Abs(Trunc(Val*20. )) < MaxADCValue)
                         or (Abs(Trunc(Val)) > MaxADCValue) then begin
                         Channel[ChSim].ADCScale := Abs(QuantalContent*QuantumAmplitude*1.2)
                                               / MaxADCValue ;
                         end ;
                      rH^.ADCVoltageRange[ChSim] :=  Channel[ChSim].ADCCalibrationFactor
                                * ( Channel[ChSim].ADCScale * (MaxADCValue+1) ) ;
                      end ;

                 Channel[ChSim].ADCZero := 0 ;
                 Channel[ChSim].xMin := 0 ;
                 Channel[ChSim].xMax := RawfH.NumSamples-1 ;
                 Channel[ChSim].yMin := MinADCValue ;
                 Channel[ChSim].yMax := MaxADCValue ;
                 Channel[ChSim].InUse := True ;

                 FirstRecord := False ;
                 end ;

              { Clear all channels }
              for ch := 0 to RawFH.NumChannels-1 do
                  for i := 0 to RawFH.NumSamples-1 do ADC^[i] := 0 ;

              { Create simulated synaptic current }

              { Set EPC amplitude }
              QuantalContent := Binomial( p, n ) ;
              if QuantalContent > 0. then begin
                 if rbCurrents.checked then begin
                    Amplitude := (QuantumAmplitude*QuantalContent) +
                        GaussianRandom(AmpTemp) * (QuantumStDev*sqrt(QuantalContent))
                    end
                 else begin
                    { Calculate amplitude accouting for non-linear summation}
                    VRest := Abs(ExtractFloat( edVRest.text, -90. )) ;
                    Amplitude := VRest -
                                 (VRest / ( 1.+ QuantalContent*(
                                         (VRest/(VRest-Abs(QuantumAmplitude)))-1. ))) ;
                    { Add variance scaled for non-linear sum.}
                    Amplitude := Amplitude
                                 + Abs(Amplitude/(QuantalContent*QuantumAmplitude))
                                   *QuantumStDev*GaussianRandom(AmpTemp) ;
                    end ;
                 end
              else begin
                   Amplitude := 0. ;
                   end ;

              ChOffset := RawFH.NumChannels - 1 - ChSim ;
              NoiseTemp := 1. ;
              x := 0. ;
              iStartEPC := RawFH.NumSamples div 10 + Trunc(-Latency*ln(random)/RawFH.dt);
              for i := 0 to RawFH.NumSamples-1 do begin

                  { Create background noise with gaussian distribution }
                  y := GaussianRandom(NoiseTemp)*NoiseRMS ;

                  if i >= iStartEPC then begin
                     ySig := Amplitude*(A1*exp(-x/Tau1))  ;
                     if ck2Exponential.checked then
                        ySig := ySig + Amplitude*(A2*exp(-x/Tau2)) ;
                     if TauRise > 0. then ySig := ySig*0.5*(1.+erf( (x/TauRise)-3. ) ) ;
                     y := y + ySig ;
                     x := x + Rawfh.dt ;
                     end ;
                  j := RawFH.NumChannels*i + ChOffset ;
                  ADC^[j] := Trunc(y/Channel[ChSim].ADCScale) + Channel[ChSim].ADCZero ;
                  end ;

              { Display simulated record }

              InitializeDisplay( Channel, RawFH.NumChannels, rH^,
                                 lbTMin,lbTMax, pbDisplay) ;

              { Erase display }
              EraseDisplay(pbDisplay) ;

              MaxCh := RawFH.NumChannels - 1;
              for ch := 0 to RawFH.NumChannels-1 do begin
                  if Channel[ch].InUse then begin
                     dx := Channel[ch].xScale ;
                     x := -Channel[ch].xScale*Channel[ch].xMin + Channel[ch].Left ;
                     nOffScreen := 0 ;
                     for i := 0 to RawFH.NumSamples-1 do begin

                         y := ADC^[(i*RawFH.NumChannels) + (MaxCh-ch) ] ;
                         xPix := Trunc(x) ;
                         yPix := Trunc(Channel[ch].Bottom
                                 - Channel[ch].yScale*(y - Channel[ch].yMin));

                         if (xPix < Channel[ch].Left)
                            or (Channel[ch].Right < xPix )
                            or (yPix < Channel[ch].Top)
                            or (Channel[ch].Bottom < yPix ) then begin
                            xPix := IntLimitTo( xPix, Channel[ch].Left,
                                                Channel[ch].Right ) ;
                            yPix := IntLimitTo( yPix, Channel[ch].Top,
                                                Channel[ch].Bottom ) ;
                            nOffScreen := nOffScreen + 1 ;
                            end
                         else nOffScreen := 0 ;

                         if (nOffScreen > 1) or (i=0) then
                              pbDisplay.canvas.moveto(xPix,yPix)
                         else pbDisplay.canvas.lineto(xPix,yPix);

                         x := x + dx ;
                         end ;

                     { Display Channel Name }
                     pbDisplay.Canvas.TextOut( Channel[ch].Left,
                              (Channel[ch].Top + Channel[ch].Bottom) div 2,
                              ' ' + Channel[ch].ADCName ) ;

                     { Draw cursors }
                     DrawHorizontalCursor(pbDisplay,Channel[ch],
                                          Channel[ch].ADCZero) ;
                     end ;
                  end ;

              { Save Record to file }
              Inc(RawFH.NumRecords) ;
              rH^.Status := 'ACCEPTED' ;
              if (n = 1. ) and (p = 1. ) then
                   rH^.RecType := 'MINI'
              else rH^.RecType := 'EVOK' ;

              rH^.Number := RawfH.NumRecords ;
              rH^.Time := rH^.Number ;
              rH^.dt := RawfH.dt ;
              rH^.Equation.Available := False ;
              rH^.Analysis.Available := False ;
              rH^.Ident := ' ' ;
              PutRecord( RawfH, rH^, RawfH.NumRecords, ADC^ ) ;

              { Terminate when all records done }
              Inc(NumRecordsDone) ;
              edProgress.text := format( 'Rec. %d/%d Done', [NumRecordsDone,
                                          NumRecordsRequired] ) ;
              if NumRecordsDone >= NumRecordsRequired then bABort.click ;
              end ;
          end ;
     end;


procedure TSynapseSim.UpdateAmplitudes ;
begin
     edQuantumAmplitude.text := format( ' %.3g %s',
                             [ExtractFloat(edQuantumAmplitude.text,1. ),Units]);
     edQuantumStDev.text := format( ' %.3g %s',
                             [ExtractFloat(edQuantumStDev.text,0. ),Units]);
     edNoiseRMS.text := format( ' %.3g %s',[ExtractFloat(edNoiseRMS.text,0.01 ),Units]);
     end ;


procedure TSynapseSim.bStartClick(Sender: TObject);
begin
     State := DoSimulation ;
     bStart.Enabled := False ;
     bAbort.Enabled := True ;
     NumRecordsRequired := Trunc(GetFromEditBox(edNumRecords,10.,1.,1E4,'%.0f','')) ;
     NumRecordsDone := 0 ;
     FirstRecord := True ;
     end;


procedure TSynapseSim.bAbortClick(Sender: TObject);
begin
     State := Idle ;
     bStart.Enabled := True ;
     bAbort.Enabled := False ;
     if RawFH.NumRecords > 0 then SaveHeader( RawFH ) ;
     WriteToLogFile( format('(%d records)',[RawFH.NumRecords]));
     end;


procedure TSynapseSim.FormDeactivate(Sender: TObject);
begin
     Timer.Enabled := False ;
     end;


procedure TSynapseSim.bCloseClick(Sender: TObject);
begin
     HeapBuffers( Deallocate ) ; 
     close ;
     end;


function binomial( pIn, n : Single ) : Single ;
{ --------------------------------------------
  Binomial random number generator
  Returns a number from the distribution B(pIn,n)
  where pIn = probability, n = number of items
  (Base on Numerical Recipes code)
  --------------------------------------------}
var
   p,mean,r,em,g,t,oldg,pc,pclog,y,plog,sq,zz : Single ;
   i : LongInt ;
   quit : Boolean ;
begin

	if pIn > 0.5 then p := 1. - pIn
                     else p := pIn ;

	mean := n*p ;
	if n <= 25.  then begin
	    r := 0. ;
	    for i := 1 to Trunc(n) do if random < p then r := r + 1. ;
            end
	else if mean < 1. then begin
	    g := exp(-mean) ;
	    t := 1. ;
	    r := 0. ;
	    while( (r<n) and (t<g) ) do begin
		t := t*random ;
		r := r + 1. ;
                end ;
            end
	else begin
	   oldg := gammln(n+1. ) ;
	   pc := 1. - p ;
	   plog := ln(p) ;
	   pclog := ln(pc) ;
	   sq := sqrt(2.*mean*pc) ;

	   quit := False ;
	   while( not quit ) do begin
               { Make sure TAN(infinity) is not calculated }
               repeat zz := random until (Abs(zz-0.5)>0.001) ;
	       y := sin(zz*Pi) / cos(zz*Pi) ;
	       em := sq*y + mean ;
	       if (em >= 0. ) and (em < n+1. ) then begin
		   em := Int(em) ;
		   t := 1.2*sq*(1.+y*y)*exp(oldg-gammln(em+1. ) -
     		        gammln(n-em+1. ) + em*plog + (n-em)*pclog) ;
		   if( random <= t ) then quit := True ;
                   end ;
               end ;
	    r := em ;
	    end ;

	if ( p <> pIn ) then r := n - r ;
	binomial := r ;
        end ;

function gammln( xx : Single ) : Single ;
var
   stp,x,tmp,ser : Double ;
   cof : Array[1..7] of Double ;
   i : LongInt ;
begin
	cof[1] := 76.18009173 ;
	cof[2] := -86.50532033 ;
	cof[3] := 24.01409822 ;
	cof[4] := -1.231739516 ;
	cof[5] := 0.120858003E-2 ;
	cof[6] := -0.536382E-5 ;
	stp := 2.50662827465 ;

	x := (xx - 1. ) ;
	tmp := x + 5.5 ;
	tmp := ( x + 0.5)*ln(tmp) - tmp;
	ser := 1. ;
	for i := 1 to 6 do begin
	    x := x + 1. ;
	    ser := ser + cof[i]/x ;
	    end ;
	tmp := tmp + ln(stp*ser) ;
	gammln := tmp ;
        end ;


procedure TSynapseSim.ck2ExponentialClick(Sender: TObject);
var
   Tau : single ;
begin
     Tau := ((RawFH.dt*RawFH.NumSamples) / 3. ) * Settings.TScale ;
     edTau2.text := format( ' %.3g %s', [Tau,Settings.TUnits]) ;
     if edTau2.text = ' ' then
        edTau2.text := format( ' %.3g %s', [Tau,Settings.TUnits]) ;
     if edRatio.text = ' ' then edRatio.text := ' 2 ' ;
     end;

procedure TSynapseSim.rbCurrentsClick(Sender: TObject);
begin
     Units := 'nA' ;
     UpdateAmplitudes ;
     end;

procedure TSynapseSim.rbPotentialsClick(Sender: TObject);
begin
     Units := 'mV' ;
     UpdateAmplitudes ;
     end;
procedure TSynapseSim.FormClose(Sender: TObject; var Action: TCloseAction);
begin
     if RawFH.NumRecords > 0 then begin
        Main.ShowRaw.Enabled := True ;
        Main.ShowRaw.checked := True ;
        Main.ShowAveraged.checked :=  False ;
        Main.ShowAveraged.Enabled :=  False ;
        Main.ShowLeakSubtracted.checked := False ;
        Main.ShowLeakSubtracted.Enabled := False ;
        FH := RawFH ;
        Main.UpdateWindows ;
        end ;

     { Enable "Synaptic Current" item in "Simulation" menu }
     Main.Synapse.enabled := true ;
     Action := caFree ;
     end;


procedure TSynapseSim.FormActivate(Sender: TObject);
begin
     Timer.enabled := True ;
     end;


procedure TSynapseSim.edQuantumAmplitudeKeyPress(Sender: TObject;
  var Key: Char);
begin
     if key = chr(13) then UpdateAmplitudes ;
     end;

procedure TSynapseSim.edLatencyKeyPress(Sender: TObject; var Key: Char);
var
   r : single ;
begin
     if key = chr(13) then begin
        r := ExtractFloat( Tedit(Sender).text, 1.0 ) ;
        Tedit(Sender).text := format( ' %.3g %s', [r,Settings.TUnits] ) ;
        end ;
     end;

end.
