unit Maxill_Utils;

{
/* Copyright 1998-2011 Research Foundation State University of New York */

/* This file is part of QuB.                                            */

/* QuB is free software; you can redistribute it and/or modify          */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or    */
/* (at your option) any later version.                                  */

/* QuB is distributed in the hope that it will be useful,               */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of       */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        */
/* GNU General Public License for more details.                         */

/* You should have received a copy of the GNU General Public License,   */
/* named LICENSE.txt, in the QuB program directory.  If not, see        */
/* <http://www.gnu.org/licenses/>.                                      */
}
(*
 begin_html
 See also:
 <ul>
 <li><a href="maxill_iface.pas.html">maxill_iface.pas</a></li>
 <li><a href="Alloc.pas.html">Alloc.pas</a></li>
 <li><a href="AllocTypes.pas.html">AllocTypes.pas</a></li>
 </ul>
 <b>Up:</b> <a href="index.html">Index</a>
 end_html  
*)

interface

uses

Alloc, AllocTypes, Classes, SysUtils,

QUB_Node,
maxill_iface;

type

TMaxill_Model = class(TObject)
private
  procedure prAssignVariables;

public
  Model            :INode;
  Ns               :Integer;
  Clazz            :TPAInteger;
  P0               :TPADouble;
  K0, K1           :TPMDouble;
  Ligand, Voltage  :TPMInteger;
  Variables        :TStringList;

  constructor Create(Model :INode);
  destructor Destroy; override;

  procedure UpdateModelRates;
end;

TMaxill_Interval_Data = class(TObject)
private
  prDataSet              :INode;

public
  ExpCond                :INode;
  Constants              :TPADouble;
  Tdead_sec              :Double;
  Nseg                   :Integer;
  DwellCounts            :TPAInteger;
  Classeses              :TPAPointer;
  Durationses            :TPAPointer;

  constructor Create(DataSet :INode; Model :TMaxill_Model);
  destructor Destroy; override;
end;

function inter_ll_obj(Model :TMaxill_Model; Data :TMaxill_Interval_Data): Double;


type

TMaxill_Subi_Data = class(TObject)
private
  prSignames             :TStringList;
  prSignals              :array of INode;
  prExpCond              :INode;
  prConstSig             :INode;

  procedure prFree;
  procedure prAddSignal(Nm :String; Idl :INode);

public
  Nsignal_in             :Integer;
  Tdead_sec              :Double;
  Nseg, Nsignal          :Integer;
  DwellCountses          :TPMInteger;
  Classeseses, Durationseses, Ampseses :TPAPointer;
  DwellCounts            :TPAInteger;
  Classeses              :TPAPointer;
  Durationses            :TPAPointer;
  Nplex, Nstim           :Integer;
  Plexiclasses           :TPAInteger;
  Stimclasses            :TPADouble;
  // output from Subi_Run
  SegLL                  :TPADouble;

  constructor Create(DataSet :INode);
  destructor Destroy; override;

  procedure Multiplex(Model :TMaxill_Model);
end;


function subi_ll_obj(Model :TMaxill_Model; Data :TMaxill_Subi_Data; UsePeq :Boolean;
                     MsgOut :TMsgOutFunc; MsgOutData :Pointer;
                     StopFlag :TPAInteger): Double;

(*

Model := ...
Data  := ...
...
Data.Multiplex(Model);

try
  LL := subi_ll_obj(Model, Data);
except ex :Exception
  ...
end;

*)

type

TMaxill_Subi_Run = class(TObject)
private
  procedure prFree;
  procedure prFreeData;
  function prOutOfBounds: Boolean;

public
  Model          :TMaxill_Model;
  ModelBk        :TMaxill_Model;
  Datas          :Array of TMaxill_Subi_Data;
  Nk, Ncns, Np   :Integer;
  K, P, S        :TPADouble;
  Acns, Ascal    :TPMDouble;
  AcnsI, AscalI  :TPMDouble;
  Bcns, Bscal    :TPADouble;

  // options for you to set
  UsePeq         :Boolean;     // True
  SearchLimit    :Double;      // 1000.0
  FreeData       :Boolean;     // True
  MsgOut         :TMsgOutFunc; // nil
  MsgOutData     :Pointer;     // nil
  StopFlag       :TPAInteger;  // nil (point to your stop flag)

  constructor Create(Model :TMaxill_Model);
  destructor Destroy; override;

  procedure AddData(Data :TMaxill_Subi_Data);
  function Setup: Boolean; // True on fail // after all data added

  function SetPars(NewS :TPADouble): Boolean; // True on fail
  function LL(NewS :TPADouble): Double;
end;

implementation

uses

QUB_Global;


constructor TMaxill_Model.Create(Model :INode);
var
  X, K, V    :INode;
  i, j       :Integer;
  VarSort    :TStringList;
begin
  Self.Model := Model;
  Clazz      := nil;
  P0         := nil;
  K0         := nil;
  K1         := nil;
  Ligand     := nil;
  Voltage    := nil;
  Variables  := nil;

  if Model['ChannelCount'].Data[0] > 1 then begin
    raise Exception.Create('Channel count > 1 is not (yet) supported.');
  end;

  VarSort    := TStringList.Create;
  VarSort.Sorted := True;
  VarSort.Duplicates := dupIgnore;

  Ns := 0;
  X  := Model.Find('States').Find('State');
  while not X.IsNull do begin
    Inc(Ns);
    X := X.SiblingSameName;
  end;

  Clazz   := IAlloc1(Ns);
  P0      := DAlloc1(Ns);
  K0      := DAlloc2(Ns, Ns);
  K1      := DAlloc2(Ns, Ns);
  Ligand  := IAlloc2(Ns, Ns);
  Voltage := IAlloc2(Ns, Ns);

  Ns      := 0;
  X       := Model.Find('States').Find('State');
  while not X.IsNull do begin
    Clazz[Ns] := X['Class'].Data[0];
    P0[Ns]    := X['Pr'].Data[0];
    Inc(Ns);
    X := X.SiblingSameName;
  end;

  for i := 0 to Ns - 1 do begin
    for j := 0 to Ns - 1 do begin
      K0[i,j]      := 0.0;
      K1[i,j]      := 0.0;
      Ligand[i,j]  := 0;
      Voltage[i,j] := 0;
    end;
  end;

  X       := Model.Find('Rates').Find('Rate');
  while not X.IsNull do begin
    i              := X['States'].Data[0];
    j              := X['States'].Data[1];
    K              := X['k0'];
    K0[i,j]        := K.Data[0];
    K0[j,i]        := K.Data[1];
    K              := X['k1'];
    K1[i,j]        := K.Data[0];
    K1[j,i]        := K.Data[1];
    K              := X['P'];
    V              := X.Find('PNames').Find('PName');
    if 0<>K.Data[0] then VarSort.Add(V.DataAsString);
    if 0<>K.Data[1] then VarSort.Add(V.Sibling.DataAsString);
    K              := X['Q'];
    V              := X.Find('QNames').Find('QName');
    if 0<>K.Data[0] then VarSort.Add(V.DataAsString);
    if 0<>K.Data[1] then VarSort.Add(V.Sibling.DataAsString);

    X              := X.SiblingSameName;
  end;

  Variables  := TStringList.Create;
  Variables.Add('Current');
  for i := 0 to VarSort.Count - 1 do begin
    Variables.Add(VarSort[i]);
  end;
  VarSort.Free;

  prAssignVariables;
end;


destructor TMaxill_Model.Destroy;
begin
  try
    Variables.Free;
    Free2(TPMChar(Voltage));
    Free2(TPMChar(Ligand));
    Free2(TPMChar(K1));
    Free2(TPMChar(K0));
    FreeMem(P0);
    FreeMem(Clazz);
  finally
    inherited;
  end;
end;

procedure TMaxill_Model.prAssignVariables;
var
  X, K, V    :INode;
  Nm         :String;
  i, j       :Integer;
begin
  X       := Model.Find('Rates').Find('Rate');
  while not X.IsNull do begin
    i               := X['States'].Data[0];
    j               := X['States'].Data[1];
    K               := X['P'];
    V               := X.Find('PNames').Find('PName');
    Nm              := V.DataAsString;
    if (K.Data[0] <> 0) and (Nm <> '') and (Variables.IndexOf(Nm) >= 0) then begin
      Ligand[i,j]   := Variables.IndexOf(Nm);
    end else begin
      Ligand[i,j]   := 0;
    end;
    Nm              := V.Sibling.DataAsString;
    if (K.Data[1] <> 0) and (Nm <> '') and (Variables.IndexOf(Nm) >= 0) then begin
      Ligand[j,i]   := Variables.IndexOf(Nm);
    end else begin
      Ligand[j,i]   := 0;
    end;
    K               := X['Q'];
    V               := X.Find('QNames').Find('QName');
    Nm              := V.DataAsString;
    if (K.Data[0] <> 0) and (Nm <> '') and (Variables.IndexOf(Nm) >= 0) then begin
      Voltage[i,j]   := Variables.IndexOf(Nm);
    end else begin
      Voltage[i,j]   := 0;
    end;
    Nm              := V.Sibling.DataAsString;
    if (K.Data[1] <> 0) and (Nm <> '') and (Variables.IndexOf(Nm) >= 0) then begin
      Voltage[j,i]   := Variables.IndexOf(Nm);
    end else begin
      Voltage[j,i]   := 0;
    end;

    X               := X.SiblingSameName;
  end;
end;


procedure TMaxill_Model.UpdateModelRates;
var
  X, K       :INode;
  i, j       :Integer;
begin
  X       := Model.Find('Rates').Find('Rate');
  while not X.IsNull do begin
    i              := X['States'].Data[0];
    j              := X['States'].Data[1];
    K              := X['k0'];
    K.Data[0]      := K0[i,j];
    K.Data[1]      := K0[j,i];
    K              := X['k1'];
    K.Data[0]      := K1[i,j];
    K.Data[1]      := K1[j,i];

    X              := X.SiblingSameName;
  end;
end;



constructor TMaxill_Interval_Data.Create(DataSet :INode; Model :TMaxill_Model);
var
  Seg                    :INode;
  i                      :Integer;
begin
  prDataSet              := DataSet; // to reserve the arrays in memory
  ExpCond                := DataSet['ExpCond'];
  Tdead_sec              := 1e-3 * ExpCond['DeadTime'].Data[0];
  Nseg                   := 0;
  Constants              := nil;
  DwellCounts            := nil;
  Classeses              := nil;
  Durationses            := nil;

  // first count
  Seg                    := DataSet.Find('Segment');
  while not Seg.IsNull do begin
    if (Seg.Find('DwellCount').DataCount > 0) and (Seg['DwellCount'].Data[0] > 0) then begin
      Inc(Nseg);
    end;
    Seg                  := Seg.SiblingSameName;
  end;

  // then alloc and fill pointers into tree
  DwellCounts            := IAlloc1(Nseg);
  Classeses              := PAlloc1(Nseg);
  Durationses            := PAlloc1(Nseg);
  Seg                    := DataSet.Find('Segment');
  Nseg                   := 0;
  while not Seg.IsNull do begin
    if (Seg.Find('DwellCount').DataCount > 0) and (Seg['DwellCount'].Data[0] > 0) then begin
      DwellCounts[Nseg]  := Seg['DwellCount'].Data[0];
      Classeses[Nseg]    := Seg['Classes'].DataPtr;
      Durationses[Nseg]  := Seg['Durations'].DataPtr;
      Inc(Nseg);
    end;
    Seg                  := Seg.SiblingSameName;
  end;

  // and apply dead time etc.
  mil_prep_segments(Nseg, DwellCounts, TPMInteger(Classeses), TPMSingle(Durationses), Tdead_sec * 1e3);

  // now figure out constants
  Constants            := DAlloc1(Model.Variables.Count);
  // first constant is unused, corresponding to Current channel
  for i := 1 to Model.Variables.Count - 1 do begin
    try
      Constants[i]   := ExpCond.Find(Model.Variables[i]).Data[0];
    except
      raise Exception.Create(Format('%s is not defined in Data Properties: Data: Experimental conditions', [Model.Variables[i]]));
    end;
  end;
end;

destructor TMaxill_Interval_Data.Destroy;
begin
  try
    FreeMem(Constants);
    FreeMem(DwellCounts);
    FreeMem(Classeses);
    FreeMem(Durationses);
  finally
    inherited;
  end;
end;



function inter_ll_obj(Model :TMaxill_Model; Data :TMaxill_Interval_Data): Double;
var
  ResultCode   :Integer;
begin
  try
    ResultCode := inter_ll(Model.Ns, Model.Clazz, Model.P0, Model.K0, Model.K1,
                           Model.Ligand, Model.Voltage, Data.Constants, Data.Tdead_sec,
                           Data.Nseg, Data.DwellCounts, TPMInteger(Data.Classeses),
                           TPMSingle(Data.Durationses), Result);
  except
    ResultCode := -878;
  end;

  if ResultCode <> 0 then begin
    raise Exception.Create(Format('%d', [ResultCode]));
  end;
end;



constructor TMaxill_Subi_Data.Create(DataSet :INode);
var
  EC, Signal            :INode;
  Mapped                :array of Integer;
  ECName                :array of String;
  i                     :Integer;
begin
  prConstSig            := TNode.Create('');
  prExpCond             := DataSet['ExpCond'];
  Tdead_sec             := prExpCond['DeadTime'].Data[0] * 1e-3;
  Nseg                  := 0;
  NSignal               := 0;
  DwellCountses         := nil;
  Classeseses           := nil;
  Durationseses         := nil;
  Ampseses              := nil;
  DwellCounts           := nil;
  Nplex                 := 0;
  Nstim                 := 0;
  Classeses             := nil;
  Durationses           := nil;
  Plexiclasses          := nil;
  Stimclasses           := nil;
  SegLL                 := nil;

  // current, then only assigned signals, in order of channel index (not model order)
  prSigNames            := TStringList.Create;
  SetLength(prSignals, 0);
  NSignal_in            := 0;

  // count signals
  i            := 0;
  Signal       := DataSet.Find('Signal');
  while not Signal.IsNull do begin
    Inc(i);
    Signal     := Signal.SiblingSameName;
  end;
  // map signals
  SetLength(Mapped, i);
  SetLength(ECName, i);
  while i >= 1 do begin
    Dec(i);
    Mapped[i] := 0;
  end;
  EC           := prExpCond.Child;
  while not EC.IsNull do begin
    if (EC.Find('ChannelIndex').DataCount > 0) then begin
      i := EC['ChannelIndex'].Data[0];
      if i >= 0 then begin
        Mapped[i] := 1;
        ECName[i] := EC.Name;
      end;
    end;
    EC         := EC.Sibling;
  end;
  // add current signal
  i            := 0;
  Signal       := DataSet.Find('Signal');
  while not Signal.IsNull do begin
    if i = DataSet['ActiveChannel'].Data[0] then begin
      prAddSignal('Current', Signal);
      Break;
    end;
    Inc(i);
    Signal     := Signal.SiblingSameName;
  end;
  // add mapped signals
  i            := 0;
  Signal       := DataSet.Find('Signal');
  while not Signal.IsNull do begin
    if 0<>Mapped[i] then begin
      prAddSignal(ECName[i], Signal);
      Break;
    end;
    Inc(i);
    Signal     := Signal.SiblingSameName;
  end;
end;

destructor TMaxill_Subi_Data.Destroy;
begin
  try
    prFree;
    prSigNames.Free;
  finally
    inherited;
  end;
end;

procedure TMaxill_Subi_Data.prFree;
var
  s       :Integer;
begin
  for s := 0 to Nseg - 1 do begin
    FreeMem(Classeseses[s]);
    FreeMem(Durationseses[s]);
    FreeMem(Ampseses[s]);

    FreeMem(Classeses[s]);
    FreeMem(Durationses[s]);
  end;
  Free2(TPMChar(DwellCountses));
  FreeMem(Classeseses);
  FreeMem(Durationseses);
  FreeMem(Ampseses);
  FreeMem(Classeses);
  FreeMem(Durationses);
  FreeMem(DwellCounts);
  FreeMem(Plexiclasses);
  FreeMem(StimClasses);
  FreeMem(SegLL);
  DwellCountses := nil;
  Classeseses   := nil;
  Durationseses := nil;
  Ampseses      := nil;
  Classeses     := nil;
  Durationses   := nil;
  DwellCounts   := nil;
  Plexiclasses  := nil;
  StimClasses   := nil;
  SegLL         := nil;
  Nseg          := 0;
  Nsignal       := 0;
  Nplex         := 0;
  Nstim         := 0;
end;

procedure TMaxill_Subi_Data.prAddSignal(Nm :String; Idl :INode);
begin
  prSigNames.Add(Nm);
  SetLength(prSignals, Nsignal_in+1);
  prSignals[Nsignal_in]  := Idl;
  Inc(Nsignal_in);
end;

function BuildConstSignal(Parent :INode; Nm :String; Value :Double; Template :INode): INode;
var
  iSeg          :INode;
  oSeg          :INode;
  Durs          :TPASingle;
  Dur           :Single;
  i             :Integer;
begin
  Result        := Parent[Nm];
  iSeg          := Template.Find('Segment');
  oSeg          := TNode.CreateNull;
  while not iSeg.IsNull do begin
    Dur         := 0;
    Durs        := TPASingle(iSeg['Durations'].DataPtr);
    for i := 0 to iSeg['DwellCount'].Data[0] - 1 do begin
      Dur       := Dur + Durs[i];
    end;

    oSeg        := Result.InsertClone(oSeg, iSeg, False);
    oSeg['DwellCount'].SetData(1);
    oSeg['Classes'].SetData(0);
    oSeg['Durations'].SetupNumData(NODE_DATATYPE_SINGLE, 1, 1);
    oSeg['Durations'].Data[0] := Dur;
    oSeg['amp'].SetData(Value);

    iSeg        := iSeg.SiblingSameName;
  end;
end;

procedure TMaxill_Subi_Data.Multiplex(Model :TMaxill_Model);
var
  v, s, g      :Integer;
  X            :INode;
  Bounds       :array of TQUB_FirstLast;
begin
  prFree;

  Nsignal      := Model.Variables.Count;

  // count signal 0 segments and bounds
  X            := prSignals[0]['Idealization'].Find('Segment');
  while not X.IsNull do begin
    Inc(Nseg);
    SetLength(Bounds, Nseg);
    Bounds[Nseg-1].First := X.Data[0];
    Bounds[Nseg-1].Last  := X.Data[1];
    X          := X.SiblingSameName;
  end;

  if Nseg = 0 then begin
    raise Exception.Create('No idealized data.  Check the Data Source.');
  end;

  // alloc per-segseses
  // fill pointers[0]
  SegLL        := DAlloc1(Nseg);
  DwellCountses   := IAlloc2(Nseg, Nsignal);
  Classeseses     := PAlloc1(Nseg);
  Durationseses   := PAlloc1(Nseg);
  Ampseses        := PAlloc1(Nseg);
  X               := prSignals[0]['Idealization'].Find('Segment');
  for s := 0 to Nseg - 1 do begin
    Classeseses[s]      := PAlloc1(Nsignal);
    Durationseses[s]    := PAlloc1(Nsignal);
    Ampseses[s]         := PAlloc1(Nsignal);
    DwellCountses[s][0] := X['DwellCount'].Data[0];
    TPAPointer(Classeseses[s])[0]   := X['Classes'].DataPtr;
    TPAPointer(Durationseses[s])[0] := X['Durations'].DataPtr;
    TPAPointer(Ampseses[s])[0]      := X['amp'].DataPtr;
    X                   := X.SiblingSameName;
  end;

  for v := 1 to Model.Variables.Count - 1 do begin
    g     := prSigNames.IndexOf(Model.Variables[v]);
    if g >= 0 then begin
      X   := prSignals[g]['Idealization'].Find('Segment');
    end else begin
      X := prConstSig.Find(Model.Variables[v]);
      if X.IsNull then begin
        try
          X := BuildConstSignal(prConstSig, Model.Variables[v],
                                prExpCond.Find(Model.Variables[v]).Data[0],
                                prSignals[0]['Idealization']);
        except
          raise Exception.Create(Format('%s is not defined in Data Properties: Data: Experimental conditions', [Model.Variables[v]]));
        end;
      end;
      X := X.Find('Segment');
    end;

    for s := 0 to Nseg - 1 do begin
      //      confirm Nseg and bounds or exception
      if X.IsNull or (X.Data[0] <> Bounds[s].First) or (X.Data[1] <> Bounds[s].Last) then begin
        raise Exception.Create(Format('%s signal needs to be idealized; try the Stimulus preprocessor.', [Model.Variables[v]]));
      end;
      //      fill pointers[v+1]
      DwellCountses[s][v] := X['DwellCount'].Data[0];
      TPAPointer(Classeseses[s])[v]   := X['Classes'].DataPtr;
      TPAPointer(Durationseses[s])[v] := X['Durations'].DataPtr;
      TPAPointer(Ampseses[s])[v]      := X['amp'].DataPtr;
      X                     := X.SiblingSameName;
    end;

  end;

  DwellCounts             := IAlloc1(Nseg);

  msl_multiplex_bounds(Nseg, NSignal, DwellCountses, Classeseses, Durationseses,
                       Ampseses, DwellCounts, Nplex, Nstim);

  // alloc classeses[Nseg][Nd], durationses
  //       plexiclasses[2*Nplex], stimclasses[NsigConst*NStim]
  Classeses               := PAlloc1(Nseg);
  Durationses             := PAlloc1(Nseg);
  for s := 0 to Nseg - 1 do begin
    Classeses[s]          := IAlloc1(DwellCounts[s]);
    Durationses[s]        := SAlloc1(DwellCounts[s]);
  end;
  Plexiclasses            := IAlloc1(2*Nplex);
  Stimclasses             := DAlloc1(Nstim*NSignal);

  msl_multiplex(Nseg, Nsignal, DwellCountses, Classeseses, Durationseses, Ampseses,
                TDead_Sec*1e3, DwellCounts, TPMInteger(Classeses), TPMSingle(Durationses),
                Nplex, Plexiclasses, Nstim, Stimclasses);
end;


function subi_ll_obj(Model :TMaxill_Model; Data :TMaxill_Subi_Data; UsePeq :Boolean;
                     MsgOut :TMsgOutFunc; MsgOutData :Pointer; StopFlag :TPAInteger): Double;
var
  ResultCode   :Integer;
  P0           :TPADouble;
begin
  if   UsePeq
  then P0      := nil
  else P0      := Model.P0;
  try
    ResultCode := subi_ll(Model.Ns, Model.Clazz, P0, Model.K0, Model.K1,
                          Model.Ligand, Model.Voltage, Data.Tdead_sec, Data.Nseg,
                          Data.DwellCounts, TPMInteger(Data.Classeses),
                          TPMSingle(Data.Durationses),
                          Data.Nsignal, Data.Nplex, Data.Plexiclasses,
                          Data.Nstim, Data.Stimclasses, Result, nil, Data.SegLL,
                          MsgOut, MsgOutData, StopFlag);
  except
    ResultCode := -878;
  end;

  if ResultCode = -33 then begin
    // Stopped by StopFlag
  end else if ResultCode = -1 then begin
    raise Exception.Create('Not enough classes in the model.');
  end else if ResultCode = -2 then begin
    raise Exception.Create('Bad eigen decomposition.');
  end else if ResultCode = -878 then begin
    raise Exception.Create('Unknown exception in maxill.dll');
  end else if ResultCode <> 0 then begin
    raise Exception.Create(Format('Error %d in maxill.dll', [ResultCode]));
  end;
end;



constructor TMaxill_Subi_Run.Create(Model :TMaxill_Model);
begin
  Self.Model        := Model;
  UsePeq            := True;
  SearchLimit       := 1000.0;
  FreeData          := True;
  MsgOut            := nil;
  MsgOutData        := nil;
  StopFlag          := nil;
  SetLength(Datas, 0);
  Nk                := 0;
  Ncns              := 0;
  Np                := 0;
  K                 := nil;
  P                 := nil;
  S                 := nil;
  Acns              := nil;
  Ascal             := nil;
  AcnsI             := nil;
  AscalI            := nil;
  Bcns              := nil;
  Bscal             := nil;
  ModelBk           := nil;
end;

destructor TMaxill_Subi_Run.Destroy;
begin
  try
    prFree;
    if FreeData then prFreeData;
  finally
    inherited;
  end;
end;

procedure TMaxill_Subi_Run.prFree;
begin
  FreeMem(K);             K          := nil;
  FreeMem(P);             P          := nil;
  FreeMem(S);             S          := nil;
  Free2(TPMChar(Acns));   Acns       := nil;
  Free2(TPMChar(Ascal));  Ascal      := nil;
  Free2(TPMChar(AcnsI));  AcnsI      := nil;
  Free2(TPMChar(AscalI)); AscalI     := nil;
  FreeMem(Bcns);          Bcns       := nil;
  FreeMem(Bscal);         Bscal      := nil;
  ModelBk.Free;           ModelBk    := nil;
end;

procedure TMaxill_Subi_run.prFreeData;
var
  i          :Integer;
begin
  for i := 0 to Length(Datas) - 1 do begin
    Datas[i].Free;
  end;
  SetLength(Datas, 0);
end;

procedure TMaxill_Subi_Run.AddData(Data :TMaxill_Subi_Data);
var
  i          :Integer;
begin
  i          := Length(Datas);
  SetLength(Datas, i+1);
  Datas[i] := Data;
end;

function TMaxill_Subi_Run.Setup: Boolean; // True on fail // after all data added
var
  Ain            :TPMDouble;
  Bin            :TPADouble;
  Cns            :INode;
  Cname          :String;
  Coeffs         :INode;
  Rate           :INode;
  Aentry         :Array of Double;
  i, ci          :Integer;
  a, b           :Integer;
begin
  prFree;

  Nk             := Maxill_CountK(Model.Ns);
  Ain            := nil;
  Bin            := nil;

  try

    // count user constraints to alloc an upper bound of storage
    Ncns           := 0;
    Cns            := Model.Model.Find('Constraints').Child;
    while not Cns.IsNull do begin
      Inc(Ncns);
      Cns          := Cns.Sibling;
    end;
    Ncns           := 2*Ncns + Nk;

    Ain            := DAlloc2(Ncns, Nk);
    Bin            := DAlloc1(Ncns);
    Ncns           := 0; // now it's the count of filled-in constraints

    Maxill_AddAutoConstraints(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                              Ncns, Ain, Bin);

    Cns            := Model.Model.Find('Constraints').Child;
    while not Cns.IsNull do begin
      Cname        := Cns.Name;
      if Cname = 'FixRate' then begin
        Maxill_AddFixRateConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                    Ncns, Ain, Bin, Cns.Data[0], Cns.Data[1]);
      end else if Cname = 'FixExp' then begin
        Maxill_AddFixExpConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                   Ncns, Ain, Bin, Cns.Data[0], Cns.Data[1]);
      end else if Cname = 'ScaleRate' then begin
        Maxill_AddScaleRateConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                      Ncns, Ain, Bin, Cns.Data[0], Cns.Data[1], Cns.Data[2], Cns.Data[3]);
      end else if Cname = 'ScaleExp' then begin
        Maxill_AddScaleExpConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                     Ncns, Ain, Bin, Cns.Data[0], Cns.Data[1], Cns.Data[2], Cns.Data[3]);
      end else if Cname = 'LoopBal' then begin
        if 0<>Maxill_AddLoopBalConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                    Ncns, Ain, Bin, Cns.DataCount, Cns.DataPtr,
                                    MsgOut, MsgOutData) then begin
          Result := True;
          Exit;
        end;
      end else if Cname = 'LoopImbal' then begin
        if 0<>Maxill_AddLoopImbalConstraint(Model.Ns, Model.K0, Model.K1, Model.Ligand, Model.Voltage,
                                      Ncns, Ain, Bin, Cns.DataCount, Cns.DataPtr,
                                      MsgOut, MsgOutData) then begin
          Result := True;
          Exit;
        end;
      end else if Cname = 'Generalized' then begin
        Coeffs := Cns['Coeffs'];
        if Coeffs.DataCount > 0 then begin
          SetLength(Aentry, Maxill_CountK(Model.Ns));
          for i := Maxill_CountK(Model.Ns)-1 downto 0 do begin
            Aentry[i] := 0.0;
          end;
          // re-order from QUB_Q_Model order (interleaved k0, k1 in Rates order) to maxill order (all k0, then k1, in canonical order)
          ci   := 0;
          Rate := Model.Model.Find('Rates').Find('Rate');
          while not Rate.IsNull do begin
            a := Rate['States'].Data[0];
            b := Rate['States'].Data[1];
            Aentry[ Maxill_k_ix_of(Model.Ns, a, b)]        := Coeffs.Data[ci];
            Aentry[ Maxill_k_ix_of(Model.Ns, a, b, True) ] := Coeffs.Data[ci+1];
            Aentry[ Maxill_k_ix_of(Model.Ns, b, a)]        := Coeffs.Data[ci+2];
            Aentry[ Maxill_k_ix_of(Model.Ns, b, a, True) ] := Coeffs.Data[ci+3];
            ci   := ci + 4;
            Rate := Rate.SiblingSameName;
          end;
          Maxill_AddGeneralizedConstraint(Model.Ns, Ncns, Ain, Bin, @(Aentry[0]), Cns['Value'].Data[0]);
        end;
      end;
      Cns          := Cns.Sibling;
    end;

    Maxill_ReduceConstraints(Model.Ns, Ncns, Ain, Bin);

    Np             := Nk - Ncns;
    K              := DAlloc1(Nk);
    P              := DAlloc1(Np);
    S              := DAlloc1(Np);
    Acns           := DAlloc2(Nk, Np);
    AcnsI          := DAlloc2(Np, Nk);
    Ascal          := DAlloc2(Np, Np);
    AscalI         := DAlloc2(Np, Np);
    Bcns           := DAlloc1(Nk);
    Bscal          := DAlloc1(Np);

    Maxill_K01toK(Model.Ns, Model.K0, Model.K1, K);
    Result := 0<>Maxill_SetupLinearConstraints(Ncns, Ain, Bin, Nk, K, P, Acns, Bcns, AcnsI, MsgOut, MsgOutData);
    Maxill_SetupStartAtOnes(Np, P, Ascal, Bscal, AscalI);
    for i := 0 to Np - 1 do begin
      S[i] := 1.0;
    end;

    ModelBk        := TMaxill_Model.Create(Model.Model.Clone(True));

  finally
    Free2(TPMChar(Ain));
    FreeMem(Bin);
  end;
end;



function TMaxill_Subi_Run.SetPars(NewS :TPADouble): Boolean;
var
  i              :Integer;
  Msg            :String;
begin
  Result         := False;
  for i := 0 to Np - 1 do begin
    S[i]         := NewS[i];
  end;
  Maxill_ParsToK(Np, Np, S, Ascal, AscalI, Bscal, P);
  Maxill_ParsToK(Np, Nk, P, Acns, AcnsI, Bcns, K);
  i := Maxill_KtoK01(Model.Ns, K, Model.K0, Model.K1);
  if i <> 0 then begin
    Result       := True;
    if Assigned(MsgOut) then begin
      Msg        := Format('MSL: Parameter %d out of bounds: %f', [(-i)-1, K[(-i)-1]]);
      MsgOut(PChar(Msg), MsgOutData);
    end;
  end;
end;

function TMaxill_Subi_Run.LL(NewS :TPADouble): Double;
var
  d              :Integer;
begin
  Result         := 0.0;

  if SetPars(NewS)
  or prOutOfBounds
  then Exit;

  for d := 0 to Length(Datas) - 1 do begin
    Result     := Result + subi_ll_obj(Model, Datas[d], UsePeq, MsgOut, MsgOutData, StopFlag);
  end;
end;

function TMaxill_Subi_Run.prOutOfBounds: Boolean;
var
  i, j     :Integer;
begin
  Result   := False;

  for i := 0 to Model.Ns - 1 do begin
    for j := 0 to Model.Ns - 1 do begin
      if (i <> j) then begin
        if (Abs(Model.K0[i,j]) < Abs(ModelBk.K0[i,j] / SearchLimit))
        or (Abs(Model.K0[i,j]) > Abs(ModelBk.K0[i,j] * SearchLimit))
        or (Abs(Model.K1[i,j]) < Abs(ModelBk.K1[i,j] / SearchLimit))
        or (Abs(Model.K1[i,j]) > Abs(ModelBk.K1[i,j] * SearchLimit)) then begin
          Result := True;
          Exit;
        end;
      end;
    end;
  end;
end;





end.