(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialiation Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

BeginPackage["SegAnalysis`ThreeOne`", "SegAnalysis`jaco311`", 
  "SegAnalysis`jaco312`"]

ThreeOneA::usage = 
  "ThreeOneA[x, x1, t1, x2, t2, x3, t3, x4, t4] gives the expected fraction \
of A genotype at x, when three segregation distortion genes, {x1, t1} {x2, \
t2} {x3, t3}, are on the same gamete, and another gene, {x4, t4}, are on the \
other gamate."

ThreeOneB::usage = 
  "ThreeOneB[x, x1, t1, x2, t2, x3, t3, x4, t4] gives the expected fraction \
of B genotype at x, when three segregation distortion genes, {x1, t1} {x2, \
t2} {x3, t3}, are on the same gamete, and another gene, {x4, t4}, are on the \
other gamate."

ThreeOneH::usage = 
  "ThreeOneH[x, x1, t1, x2, t2, x3, t3, x4, t4] gives the expected fraction \
of H genotype at x, when three segregation distortion genes, {x1, t1} {x2, \
t2} {x3, t3}, are on the same gamete, and another gene, {x4, t4}, are on the \
other gamate."

Begin["`Private`"]

(* === Warning messages ====  *)

SegregAnalysis::lmincrease = 
  "Warning: after `` decreases in the step size, NonlinearFit failed to find \
a new parameter estimate that will decrease the sum of squares."

(* ==================    ThreeOne   ================= *)
(* This uses Guess-Newton method for finding the least-
    squares solution for segregation distortion for four genes model *)

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, x1_, t1_, 
    x2_, t2_, x3_, t3_, x4_, t4_, fp_, maxits_, wp_, pg_, ag_, mdls_, zt_]:=
  Module[{predA, predB, predH, chilistA, chilistB, chilistH, tmpchilistA, 
      tmpchilistB, tmpchilistH, tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, 
      chisq,derx1, dert1, derx2, dert2,derx3, dert3, derx4, dert4, der, 
      alpha, beta, chisq0, newchi,tmpchi, its, accDelta, NAcc, precDelta, 
      NPrec, invaliddigits, accflag0, accflag1 = False, accflag2 = False, 
      precflag0, precflag1 = False, precflag2 = False, ovflag = False, 
      guessx1 = x1, guesst1 = t1, guessx2 = x2, guesst2 = t2, guessx3 = x3, 
      guesst3 = t3, guessx4 = x4, guesst4 = t4, oldx1, oldt1, oldx2, oldt2, 
      oldx3, oldt3, oldx4, oldt4, delta, deltax1, deltat1, deltax2, deltat2, 
      deltax3, deltat3, deltax4, deltat4, tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, 
      tmpt3, tmpx4, tmpt4, dits, dmaxits, newpredA, newpredB, newpredH, 
      newchilistA, newchilistB, newchilistH, accgoal = N[10^(-ag)], 
      precgoal = N[10^(-pg)], paramerr, errorDOF},
			
			(* calculating chisq *)
			{{chilistA, chilistB, chilistH}, chisq} = 
      calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, nresponseH}, {
          guessx1, guesst1, guessx2, guesst2, guessx3, guesst3, guessx4, 
          guesst4}, wp];
			tmpchi = chisq; newchi = -1;
			its = 1; oldx1 = guessx1; oldt1 = guesst1;
			oldx2 = guessx2; oldt2 = guesst2;
			oldx3 = guessx3; oldt3 = guesst3;
			oldx4 = guessx4; oldt4 = guesst4;
			
			(* main loop for twoone *)
			ovflag = Check[
					While[((accDelta = Abs[chisq - newchi];
								NAcc = -Floor[Log[10, accDelta]] - 1;
								precDelta = Abs[chisq - newchi]/(newchi + chisq)/2;
								NPrec = -Floor[Log[10, precDelta]] - 1;
								invaliddigits = 
                  TrueQ[NAcc =!= Infinity && Accuracy[accDelta] < NAcc && 
                      NPrec =!= Infinity && Precision[precDelta] < Nprec];
								!invaliddigits ) &&
								( accflag0 = ( 0 <= accDelta < accgoal);
									precflag0 =( 0 <= precDelta < precgoal);
									!((accflag0 && accflag1 && accflag2) || 
											(precflag0 && precflag1 && precflag2))  ) &&
								its <= maxits),
						{accflag2, precflag2} = {accflag1, precflag1};
						{accflag1, precflag1} = {accflag0, precflag0};
						chisq = tmpchi;
						oldx1 = guessx1; oldt1 = guesst1;
						oldx2 = guessx2; oldt2 = guesst2;
						oldx3 = guessx3; oldt3 = guesst3;
						oldx4 = guessx4; oldt4 = guesst4;
						(* find deltax1, deltat1, deltax2, deltat2, deltax3, deltat3, deltax4, 
            deltat4 using LineaSolve for narmal equation *)
						tmpchilistA = Delete[chilistA, Position[nptsA, guessx1]];
						tmpnptsA = Delete[nptsA, Position[nptsA, guessx1]];
						tmpchilistB =Delete[chilistB, Position[nptsB, guessx1]];
						tmpnptsB =Delete[nptsB, Position[nptsB, guessx1]];
						tmpchilistH =Delete[chilistH, Position[nptsH, guessx1]];
						tmpnptsH =Delete[nptsH, Position[nptsH, guessx1]];
						tmpchilistA =Delete[tmpchilistA, Position[tmpnptsA, guessx2]];
						tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessx2]];
						tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessx2]];
						tmpnptsB = Delete[tmpnptsB, Position[tmpnptsB, guessx2]];
						tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessx2]];
						tmpnptsH = Delete[tmpnptsH, Position[tmpnptsH, guessx2]];
						tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessx3]];
						tmpnptsA = Delete[tmpnptsA, Position[tmpnptsA, guessx3]];
						tmpchilistB =Delete[tmpchilistB, Position[tmpnptsB, guessx3]];
						tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessx3]];
						tmpchilistH =Delete[tmpchilistH, Position[tmpnptsH, guessx3]];
						tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessx3]];
						tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessx4]];
						tmpnptsA = Delete[tmpnptsA, Position[tmpnptsA, guessx4]];
						tmpchilistB =Delete[tmpchilistB, Position[tmpnptsB, guessx4]];
						tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessx4]];
						tmpchilistH =Delete[tmpchilistH, Position[tmpnptsH, guessx4]];
						tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessx4]];
						tmpchilist=Join[tmpchilistA, tmpchilistB, tmpchilistH];
						(* calculate jacobian *)
					
					der={};
					delta={};
					If[!(MemberQ[fp, "x1"]),
						derx1 = 
              SegAnalysis`jaco311`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,derx1]];
					If[!(MemberQ[fp, "t1"]),
						dert1 = 
              SegAnalysis`jaco311`Private`b[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,dert1]];
					If[!(MemberQ[fp, "x2"]),
						derx2 = 
              SegAnalysis`jaco311`Private`c[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,derx2]];
					If[!(MemberQ[fp, "t2"]),
						dert2 = 
              SegAnalysis`jaco311`Private`d[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,dert2]];
					If[!(MemberQ[fp, "x3"]),
						derx3 = 
              SegAnalysis`jaco312`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,derx3]];
					If[!(MemberQ[fp, "t3"]),
						dert3 = 
              SegAnalysis`jaco312`Private`b[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,dert3]];
					If[!(MemberQ[fp, "x4"]),
						derx4 = 
              SegAnalysis`jaco312`Private`c[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,derx4]];
					If[!(MemberQ[fp, "t4"]),
						dert4 = 
              SegAnalysis`jaco312`Private`d[guessx1, guesst1, guessx2, 
                guesst2, guessx3, guesst3, guessx4, guesst4, tmpnptsA, 
                tmpnptsB, tmpnptsH, wp];
						der=Append[der,dert4]];
			
					
				(* normal equation: alpha . deltaguess === beta   *)
						{alpha, beta} = 
            N[{der . Transpose[der], Map[tmpchilist . #1& , der]}, wp];
					
						delta = N[LinearSolve[alpha, beta], wp];
					If[!(MemberQ[fp, "x1"]),deltax1=First[delta];delta=Rest[delta], 
            deltax1=0];
					If[!(MemberQ[fp, "t1"]),deltat1=First[delta]; delta=Rest[delta],  
            deltat1=0];
					If[!(MemberQ[fp, "x2"]),deltax2=First[delta];delta=Rest[delta], 
            deltax2=0];
					If[!(MemberQ[fp, "t2"]),deltat2=First[delta]; delta=Rest[delta],  
            deltat2=0];
					If[!(MemberQ[fp, "x3"]),deltax3=First[delta];delta=Rest[delta], 
            deltax3=0];
					If[!(MemberQ[fp, "t3"]),deltat3=First[delta]; delta=Rest[delta],  
            deltat3=0];
					If[!(MemberQ[fp, "x4"]),deltax4=First[delta];delta=Rest[delta], 
            deltax4=0];
					If[!(MemberQ[fp, "t4"]),deltat4=First[delta],  deltat4=0];
					
						{tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, tmpt3, tmpx4, tmpt4} = 
            N[{oldx1, oldt1, oldx2, oldt2, oldx3, oldt3, oldx4, oldt4} + {
                  deltax1, deltat1, deltax2, deltat2, deltax3, deltat3, 
                  deltax4, deltat4}, wp];
						(* calculation chi-squares *)
						{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                nresponseH},  {tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, tmpt3, 
                tmpx4, tmpt4}, wp];
			
							(* find newguess for smaller chisq  *)
						dits = 1; dmaxits = 15;
						While[newchi - chisq >= 0 && dits <= dmaxits,
							{tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, tmpt3, tmpx4, tmpt4} = 
              N[{oldx1, oldt1, oldx2, oldt2, oldx3, oldt3, oldx4, 
                    oldt4} + (0.5 ^ dits)  {deltax1, deltat1, deltax2, 
                      deltat2, deltax3, deltat3, deltax4, deltat4}, wp];
							{{newchilistA, newchilistB, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                  nresponseH},  {tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, tmpt3, 
                  tmpx4, tmpt4}, wp];
							dits++];  (* end While[newchi     *)
						If[dits > dmaxits, Message[SegregAnalysis::lmincrease, dmaxits]];
						tmpchi = newchi;
						{guessx1, guesst1, guessx2, guesst2, guessx3, guesst3, guessx4, 
              guesst4} = {tmpx1, tmpt1, tmpx2, tmpt2, tmpx3, tmpt3, tmpx4, 
              tmpt4};
						chilistA = newchilistA; chilistB = newchilistB;
						chilistH = newchilistH;
						printprogress[its, 
            newchi, {guessx1, guesst1, guessx2, guesst2, guessx3, guesst3, 
              guessx4, guesst4}];
						its++], (* end While *)
					True, General::ovfl];  (* end Check *)
			(* Adjust the accuracy and precision of chisq *)
			If[Accuracy[newchi] > NAcc,
			 newchi = SetAccuracy[newchi, NAcc]];
			If[Precision[newchi] > NPrec,
			newchi = SetPrecision[newchi, NPrec]];
			errorDOF = 
      Length[chilistA] + Length[chilistB] + Length[chilistH]- 8+Length[fp];
			paramerr = Sqrt[(newchi/errorDOF) DiagonalElements[Inverse[alpha]]];
			{newchi, {guessx1, guesst1, guessx2, guesst2, guessx3, guesst3, guessx4, 
        guesst4}, paramerr, errorDOF, fp}
			]

(*   *)
ThreeOneA[x_, x1_, t1_, x2_, t2_, x3_, t3_, x4_, t4_]:=Which[
		x <= x1 && x<= x4,
		1/4*(2*t1 + (-1 + 2*t1)*Tanh[2*(x - x1)])*
  (2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)]),
		x1 < x  && x<=x2 && x <= x4,
		(((Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - Sinh[4*x1])*
       ((-1 + 4*t2)*(Cosh[12*x] + Sinh[12*x]) + 
         3*(1 + 4*t2)*(Cosh[4*(2*x + x1)] + Sinh[4*(2*x + x1)]) + 
         (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
          (Cosh[8*x2] + Sinh[8*x2]) + 
         4*(t2*(Cosh[4*x] + Sinh[4*x]) + 
            (-1 + 3*t2)*(Cosh[4*x1] + Sinh[4*x1]))*
          (Cosh[4*(x + x2)] + Sinh[4*(x + x2)])) - 
      16*t1*Cosh[2*(x - x1)]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
       (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
    (t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x <= x4,
		(((Cosh[4*x] - Cosh[4*x2] + Sinh[4*x] - Sinh[4*x2])*
       ((-1 + 4*t3)*(Cosh[12*x] + Sinh[12*x]) + 
         3*(1 + 4*t3)*(Cosh[4*(2*x + x2)] + Sinh[4*(2*x + x2)]) + 
         (Cosh[4*x] + Cosh[4*x2] + Sinh[4*x] + Sinh[4*x2])*
          (Cosh[8*x3] + Sinh[8*x3]) + 
         4*(t3*(Cosh[4*x] + Sinh[4*x]) + 
            (-1 + 3*t3)*(Cosh[4*x2] + Sinh[4*x2]))*
          (Cosh[4*(x + x3)] + Sinh[4*(x + x3)])) - 
      16*t2*Cosh[2*(x - x2)]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
       (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
    (t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x <= x4,
		1/4*(-2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
  (-2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]),
		x <= x1 && x4 < x,
		1/4*(-2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
  (-2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)]),
		x1 < x && x <= x2 && x4 < x,
		(((Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - Sinh[4*x1])*
       ((-1 + 4*t2)*(Cosh[12*x] + Sinh[12*x]) + 
         3*(1 + 4*t2)*(Cosh[4*(2*x + x1)] + Sinh[4*(2*x + x1)]) + 
         (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
          (Cosh[8*x2] + Sinh[8*x2]) + 
         4*(t2*(Cosh[4*x] + Sinh[4*x]) + 
            (-1 + 3*t2)*(Cosh[4*x1] + Sinh[4*x1]))*
          (Cosh[4*(x + x2)] + Sinh[4*(x + x2)])) - 
      16*t1*Cosh[2*(x - x1)]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
       (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
    (t4 + (1/2 - t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x4 < x,
		(((Cosh[4*x] - Cosh[4*x2] + Sinh[4*x] - Sinh[4*x2])*
       ((-1 + 4*t3)*(Cosh[12*x] + Sinh[12*x]) + 
         3*(1 + 4*t3)*(Cosh[4*(2*x + x2)] + Sinh[4*(2*x + x2)]) + 
         (Cosh[4*x] + Cosh[4*x2] + Sinh[4*x] + Sinh[4*x2])*
          (Cosh[8*x3] + Sinh[8*x3]) + 
         4*(t3*(Cosh[4*x] + Sinh[4*x]) + 
            (-1 + 3*t3)*(Cosh[4*x2] + Sinh[4*x2]))*
          (Cosh[4*(x + x3)] + Sinh[4*(x + x3)])) - 
      16*t2*Cosh[2*(x - x2)]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
       (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
    (t4 + (1/2 - t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x4 < x,
		1/4*(-2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
  (-2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)])]
		
(*   *)
ThreeOneB[x_, x1_, t1_, x2_, t2_, x3_, t3_, x4_, t4_]:=Which[
		x <= x1 && x<= x4,
		1/4*(2 - 2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
  (2 - 2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]),
		x1 < x  && x<=x2 && x <= x4,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x1] + Sinh[4*x1]))*
       (-(3*Cosh[4*x] + Cosh[4*x2] + 3*Sinh[4*x] + Sinh[4*x2])*
          (Cosh[8*x] - 3*Cosh[4*(x + x1)] + Cosh[4*(x + x2)] + 
            Cosh[4*(x1 + x2)] + Sinh[8*x] - 3*Sinh[4*(x + x1)] + 
            Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)]) + 
         16*t2*Cosh[2*(x - x2)]*Sinh[2*(x - x1)]*
          (Cosh[2*(4*x + x1 + x2)] + Sinh[2*(4*x + x1 + x2)])) + 
      16*t1*Cosh[2*(x - x1)]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
       (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
    (1 - t4 + 1/2*Tanh[2*(x - x4)] + t4*Tanh[2*(-x + x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x <= x4,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x2] + Sinh[4*x2]))*
       (-(3*Cosh[4*x] + Cosh[4*x3] + 3*Sinh[4*x] + Sinh[4*x3])*
          (Cosh[8*x] - 3*Cosh[4*(x + x2)] + Cosh[4*(x + x3)] + 
            Cosh[4*(x2 + x3)] + Sinh[8*x] - 3*Sinh[4*(x + x2)] + 
            Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)]) + 
         16*t3*Cosh[2*(x - x3)]*Sinh[2*(x - x2)]*
          (Cosh[2*(4*x + x2 + x3)] + Sinh[2*(4*x + x2 + x3)])) + 
      16*t2*Cosh[2*(x - x2)]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
       (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
    (1 - t4 + 1/2*Tanh[2*(x - x4)] + t4*Tanh[2*(-x + x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x <= x4,
		1/4*(2 - 2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
  (2 - 2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]),
		x <= x1 && x4 < x,
		1/4*(2 - 2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
  (2 - 2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)]),
		x1 < x && x <= x2 && x4 < x,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x1] + Sinh[4*x1]))*
       (-(3*Cosh[4*x] + Cosh[4*x2] + 3*Sinh[4*x] + Sinh[4*x2])*
          (Cosh[8*x] - 3*Cosh[4*(x + x1)] + Cosh[4*(x + x2)] + 
            Cosh[4*(x1 + x2)] + Sinh[8*x] - 3*Sinh[4*(x + x1)] + 
            Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)]) + 
         16*t2*Cosh[2*(x - x2)]*Sinh[2*(x - x1)]*
          (Cosh[2*(4*x + x1 + x2)] + Sinh[2*(4*x + x1 + x2)])) + 
      16*t1*Cosh[2*(x - x1)]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
       (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
    (1 - t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x4 < x,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x2] + Sinh[4*x2]))*
       (-(3*Cosh[4*x] + Cosh[4*x3] + 3*Sinh[4*x] + Sinh[4*x3])*
          (Cosh[8*x] - 3*Cosh[4*(x + x2)] + Cosh[4*(x + x3)] + 
            Cosh[4*(x2 + x3)] + Sinh[8*x] - 3*Sinh[4*(x + x2)] + 
            Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)]) + 
         16*t3*Cosh[2*(x - x3)]*Sinh[2*(x - x2)]*
          (Cosh[2*(4*x + x2 + x3)] + Sinh[2*(4*x + x2 + x3)])) + 
      16*t2*Cosh[2*(x - x2)]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
       (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
    (1 - t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x4 < x,
		1/4*(2 - 2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
  (2 - 2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)])]
	
(*   *)
ThreeOneH[x_, x1_, t1_, x2_, t2_, x3_, t3_, x4_, t4_]:=Which[
		x <= x1 && x<= x4,
		-(1/4)*(2 - 2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
   (-2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]) - 
  1/4*(-2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
   (2 - 2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]),
		x1 < x  && x<=x2 && x <= x4,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x1] + Sinh[4*x1]))*
        (-(3*Cosh[4*x] + Cosh[4*x2] + 3*Sinh[4*x] + Sinh[4*x2])*
           (Cosh[8*x] - 3*Cosh[4*(x + x1)] + Cosh[4*(x + x2)] + 
             Cosh[4*(x1 + x2)] + Sinh[8*x] - 3*Sinh[4*(x + x1)] + 
             Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)]) + 
          16*t2*Cosh[2*x - 2*x2]*Sinh[2*x - 2*x1]*
           (Cosh[2*(4*x + x1 + x2)] + Sinh[2*(4*x + x1 + x2)])) + 
       16*t1*Cosh[2*x - 2*x1]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
        (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
     (t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]) + 
    ((Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - Sinh[4*x1])*
        ((-1 + 4*t2)*(Cosh[12*x] + Sinh[12*x]) + 
          3*(1 + 4*t2)*(Cosh[4*(2*x + x1)] + Sinh[4*(2*x + x1)]) + 
          (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
           (Cosh[8*x2] + Sinh[8*x2]) + 
          4*(t2*(Cosh[4*x] + Sinh[4*x]) + 
             (-1 + 3*t2)*(Cosh[4*x1] + Sinh[4*x1]))*
           (Cosh[4*(x + x2)] + Sinh[4*(x + x2)])) - 
       16*t1*Cosh[2*x - 2*x1]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
        (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
     (1 - t4 + 1/2*Tanh[2*x - 2*x4] + t4*Tanh[2*(-x + x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x <= x4,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x2] + Sinh[4*x2]))*
        (-(3*Cosh[4*x] + Cosh[4*x3] + 3*Sinh[4*x] + Sinh[4*x3])*
           (Cosh[8*x] - 3*Cosh[4*(x + x2)] + Cosh[4*(x + x3)] + 
             Cosh[4*(x2 + x3)] + Sinh[8*x] - 3*Sinh[4*(x + x2)] + 
             Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)]) + 
          16*t3*Cosh[2*x - 2*x3]*Sinh[2*x - 2*x2]*
           (Cosh[2*(4*x + x2 + x3)] + Sinh[2*(4*x + x2 + x3)])) + 
       16*t2*Cosh[2*x - 2*x2]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
        (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
     (t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]) + 
    ((Cosh[4*x] - Cosh[4*x2] + Sinh[4*x] - Sinh[4*x2])*
        ((-1 + 4*t3)*(Cosh[12*x] + Sinh[12*x]) + 
          3*(1 + 4*t3)*(Cosh[4*(2*x + x2)] + Sinh[4*(2*x + x2)]) + 
          (Cosh[4*x] + Cosh[4*x2] + Sinh[4*x] + Sinh[4*x2])*
           (Cosh[8*x3] + Sinh[8*x3]) + 
          4*(t3*(Cosh[4*x] + Sinh[4*x]) + 
             (-1 + 3*t3)*(Cosh[4*x2] + Sinh[4*x2]))*
           (Cosh[4*(x + x3)] + Sinh[4*(x + x3)])) - 
       16*t2*Cosh[2*x - 2*x2]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
        (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
     (1 - t4 + 1/2*Tanh[2*x - 2*x4] + t4*Tanh[2*(-x + x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x <= x4,
		1/4*(-(-2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
     (2 - 2*t4 + (1 - 2*t4)*Tanh[2*(x - x4)]) + 
    (2 - 2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
     (2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)])),
		x <= x1 && x4 < x,
		1/4*(-(2 - 2*t1 + (1 - 2*t1)*Tanh[2*(x - x1)])*
     (-2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)]) + 
    (2*t1 + (-1 + 2*t1)*Tanh[2*(x - x1)])*
     (2 - 2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)])),
		x1 < x && x <= x2 && x4 < x,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x1] + Sinh[4*x1]))*
        (-(3*Cosh[4*x] + Cosh[4*x2] + 3*Sinh[4*x] + Sinh[4*x2])*
           (Cosh[8*x] - 3*Cosh[4*(x + x1)] + Cosh[4*(x + x2)] + 
             Cosh[4*(x1 + x2)] + Sinh[8*x] - 3*Sinh[4*(x + x1)] + 
             Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)]) + 
          16*t2*Cosh[2*x - 2*x2]*Sinh[2*x - 2*x1]*
           (Cosh[2*(4*x + x1 + x2)] + Sinh[2*(4*x + x1 + x2)])) + 
       16*t1*Cosh[2*x - 2*x1]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
        (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
     (t4 + (1/2 - t4)*Tanh[2*(x - x4)]) + 
    ((Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - Sinh[4*x1])*
        ((-1 + 4*t2)*(Cosh[12*x] + Sinh[12*x]) + 
          3*(1 + 4*t2)*(Cosh[4*(2*x + x1)] + Sinh[4*(2*x + x1)]) + 
          (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
           (Cosh[8*x2] + Sinh[8*x2]) + 
          4*(t2*(Cosh[4*x] + Sinh[4*x]) + 
             (-1 + 3*t2)*(Cosh[4*x1] + Sinh[4*x1]))*
           (Cosh[4*(x + x2)] + Sinh[4*(x + x2)])) - 
       16*t1*Cosh[2*x - 2*x1]*(-1 + Cosh[4*(x - x2)] + 2*Sinh[4*(x - x2)])*
        (Cosh[6*x + 6*x1 + 4*x2] + Sinh[6*x + 6*x1 + 4*x2]))*
     (1 - t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
      Sinh[4*(x1 + x2)])*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
      Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])),
		x2 < x && x <= x3 && x4 < x,
		((-(Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x2] + Sinh[4*x2]))*
        (-(3*Cosh[4*x] + Cosh[4*x3] + 3*Sinh[4*x] + Sinh[4*x3])*
           (Cosh[8*x] - 3*Cosh[4*(x + x2)] + Cosh[4*(x + x3)] + 
             Cosh[4*(x2 + x3)] + Sinh[8*x] - 3*Sinh[4*(x + x2)] + 
             Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)]) + 
          16*t3*Cosh[2*x - 2*x3]*Sinh[2*x - 2*x2]*
           (Cosh[2*(4*x + x2 + x3)] + Sinh[2*(4*x + x2 + x3)])) + 
       16*t2*Cosh[2*x - 2*x2]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
        (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
     (t4 + (1/2 - t4)*Tanh[2*(x - x4)]) + 
    ((Cosh[4*x] - Cosh[4*x2] + Sinh[4*x] - Sinh[4*x2])*
        ((-1 + 4*t3)*(Cosh[12*x] + Sinh[12*x]) + 
          3*(1 + 4*t3)*(Cosh[4*(2*x + x2)] + Sinh[4*(2*x + x2)]) + 
          (Cosh[4*x] + Cosh[4*x2] + Sinh[4*x] + Sinh[4*x2])*
           (Cosh[8*x3] + Sinh[8*x3]) + 
          4*(t3*(Cosh[4*x] + Sinh[4*x]) + 
             (-1 + 3*t3)*(Cosh[4*x2] + Sinh[4*x2]))*
           (Cosh[4*(x + x3)] + Sinh[4*(x + x3)])) - 
       16*t2*Cosh[2*x - 2*x2]*(-1 + Cosh[4*(x - x3)] + 2*Sinh[4*(x - x3)])*
        (Cosh[6*x + 6*x2 + 4*x3] + Sinh[6*x + 6*x2 + 4*x3]))*
     (1 - t4 + (-(1/2) + t4)*Tanh[2*(x - x4)]))/
  (2*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + Sinh[8*x] - 
      3*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + Sinh[4*(x + x3)] + 
      Sinh[4*(x2 + x3)])*(Cosh[8*x] + Cosh[4*(x + x3)] + Cosh[4*(x2 + x3)] + 
      Sinh[8*x] + 5*(Cosh[4*(x + x2)] + Sinh[4*(x + x2)]) + 
      Sinh[4*(x + x3)] + Sinh[4*(x2 + x3)])),
		x3 < x && x4 < x,
		-(1/4)*(2 - 2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
   (-2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)]) - 
  1/4*(-2*t3 + (-1 + 2*t3)*Tanh[2*(x - x3)])*
   (2 - 2*t4 + (-1 + 2*t4)*Tanh[2*(x - x4)])]

(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{x1,t1, x2, t2, x3, t3, x4, t4, ptsA, ptsB, ptsH, responseA, 
      responseB, responseH, predA, predB, predH, chilistA, chilistB, 
      chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{x1, t1, x2, t2, x3, t3, x4, t4} = N[guess,wp];
		{ptsA, ptsB, ptsH} = pts;
		predA = N[ThreeOneA[#1, x1, t1, x2, t2, x3, t3, x4, t4] & /@ ptsA, wp];
		predB = N[ThreeOneB[#1, x1, t1, x2, t2, x3, t3, x4, t4] & /@ ptsB, wp];
		predH = N[ThreeOneH[#1, x1, t1, x2, t2, x3, t3, x4, t4] & /@ ptsH, wp];
		chilistA = responseA - predA;
		chilistB = responseB - predB;
		chilistH = responseH - predH;
		chisq =N[ 
        Plus @@ (
            Join[Map[#1^2 &, chilistA], Map[#1^2 &, chilistB], 
              Map[#1^2&, chilistH]]), wp];
		{{chilistA, chilistB, chilistH}, chisq}
		]

(*    *)
printprogress[iteration_, chisq_, params_] := 
     Print[
    StringForm[
      "Iteration:`1` ChiSquared:`2`  {Position, Transmittance}:`3`", 
											NumberForm[iteration, 6], NumberForm[chisq, 6], 
      NumberForm[params, 6]]]

(*    *)
DiagonalElements[matrix_] := 
  Module[{k, size = Length[matrix]}, 
   Table[matrix[[k,k]], {k, 1, size}]
		]

(* =========================================== *)
End[]
EndPackage[]