(***********************************************************************
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`Single`"]

SingleA::usage =
  "SingleA[x, x1, t1] gives the expected fraction of A genotype at x, when \
the segregation distortion gene trnsmits A genotype by t1 at the gene \
positition x1."

SingleB::usage =
  "SingleB[x, x1, t1] gives the expected fraction of B genotype at x, when \
the segregation distortion gene trnsmits A genotype by t1 at the gene \
positition x1."

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."

(* =========    single     =========== *)
(* This uses Guess-Newton method for finding the least-
    squares solution for a single gene segregation distortion model  *)

mnlr[nptsA_, nptsB_, nresponseA_, nresponseB_, nresponseH_, x1_, t1_, fp_, 
    maxits_, wp_, pg_, ag_, mdls_, zt_] := 
  Module[{predA, predB, chilistA, chilistB, tmpchilistA, 
    tmpchilistB, tmpchilist, tmpnptsA, tmpnptsB, chisq, 
    derAt1, derBt1, dert1, derAx1, derBx1, derx1, 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, oldx1, oldt1, delta, deltax1, 
    deltat1, tmpt1, tmpx1, dits, dmaxits, newpredA, 
    newpredB, newchilistA, newchilistB, 
    accgoal = N[10^(-ag)], precgoal = N[10^(-pg)], 
    paramerr, errorDOF},
		
			(* - calculation chi-squars *) 
  {{chilistA, chilistB}, chisq} = 
      calcchi[{nptsA, nptsB}, {nresponseA, nresponseB,  nresponseH}, {
          guessx1, guesst1}, wp];
    tmpchi = chisq; newchi = -1; its = 1; oldx1 = guessx1; 
    oldt1 = guesst1;
		(*  Print[chisq];  *)

		(* main loop *)
	 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; 
					
				(* === find deltax1 deltat1 by LinearSolve normalmatrix ==== *)	
				(* === delete not differentiable position ==== *)	
        tmpchilistA = Delete[chilistA, Position[nptsA, guessx1]]; 
        tmpnptsA =Delete[nptsA, Position[nptsA, guessx1]]; 
        tmpchilistB = Delete[chilistB, Position[nptsB, guessx1]]; 
        tmpnptsB =Delete[nptsB, Position[nptsB, guessx1]]; 
        tmpchilist = Join[tmpchilistA , tmpchilistB];
					der={};
					delta={};
					If[!(MemberQ[fp, "x1"]),
						derAx1 = N[Map[Which[
									#1 < guessx1, 
                          
                      1/2*(Sech[2*(-#1 + guessx1)]^2 - 
                            2*guesst1*Sech[2*(-#1 + guessx1)]^2), 
                     #1 > guessx1, 
                          
                      1/2*(-Sech[2*(#1 - guessx1)]^2 + 
                            2*guesst1*Sech[2*(#1 - guessx1)]^2)
								] &, tmpnptsA], wp];
						derBx1 = N[Map[Which[
									#1 < guessx1, 
                         
                      1/2* (-Sech[2*(-#1 + guessx1)]^2 + 
                            2*guesst1*Sech[2*(-#1 + guessx1)]^2), 
                     #1 > guessx1, 
                         
                      1/2*(Sech[2*(#1 - guessx1)]^2 - 
                            2*guesst1*Sech[2*(#1 - guessx1)]^2)
								] &, tmpnptsB], wp];
						derx1 = Join[derAx1, derBx1];
						der=Append[der,derx1]];
					If[!(MemberQ[fp, "t1"]),
						derAt1 = N[Map[1/2*(1 - Tanh[2*Abs[#1 - guessx1]]) &, tmpnptsA], wp];
						 derBt1 = N[Map[1/2*(-1 + Tanh[2*Abs[#1 - guessx1]]) &, tmpnptsB], wp];
						dert1 = Join[derAt1, derBt1];
						der=Append[der, dert1]];
					
			(* 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],  deltat1=0];
			{tmpx1, tmpt1} = N[{oldx1, oldt1} +  {deltax1, deltat1}, wp];
		
			 {{newchilistA, newchilistB}, newchi} = 
            calcchi[{nptsA, nptsB}, {nresponseA, nresponseB, nresponseH}, {
                tmpx1,tmpt1}, wp];
		(* Print[{newchi, chisq}]; *)
		(* find newguess for smaller chisq  *) 
			dits = 1; dmaxits = 15; 
       While[newchi - chisq >= 0 && dits <= dmaxits, 
         {tmpx1, tmpt1} = 
              N[{oldx1, oldt1} + (0.5 ^ dits)  {deltax1, deltat1}, wp];
			 (* Print[{tmpx1, tmpt1}];	*)
			 {{newchilistA, newchilistB}, newchi} = 
              calcchi[{nptsA, nptsB}, {nresponseA, nresponseB, nresponseH}, {
                  tmpx1,tmpt1}, wp];
			dits++]; 
       If[dits > dmaxits, Message[SegregAnalysis::lmincrease, dmaxits]]; 
       tmpchi = newchi; 
       {guessx1, guesst1} = {tmpx1, tmpt1}; 
       chilistA = newchilistA; chilistB = newchilistB; 
       printprogress[its, newchi, {guessx1, guesst1}]; 
      its++], (* end of While *)
		True, General::ovfl];  (* end of 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[nresponseH]- 2+Length[fp];
    paramerr = Sqrt[ ( newchi/errorDOF) DiagonalElements[Inverse[alpha]]]; 
		{newchi, {guessx1, guesst1}, paramerr, errorDOF, fp}
	]  (* end single *)

(*   *)
SingleA[x_, x1_, 
    t1_] := ((1/2) (t1 + (1/2) Tanh[2 Abs[x - x1]] - t1 Tanh[2 Abs[x - x1]]))
	
(*  *)
SingleB[x_, x1_, t1_] := 
    1/2*(1 - t1 - 1/2*Tanh[2*Abs[x - x1]] + t1*Tanh[2*Abs[x - x1]])

(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{x1,t1, ptsA, ptsB, ptsH, responseA, responseB, responseH, predA, 
      predB, predH, chilistA, chilistB, chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{ptsA, ptsB} = pts;
		{x1, t1} =N[guess,wp];
		predA = N[SingleA[#1, x1, t1] & /@ ptsA, wp];
		predB = N[SingleB[#1, x1, t1] & /@ ptsB, wp];
		chilistA = responseA - predA;
		chilistB = responseB - predB;
		chilistH = N[ (#1 - 0.5 )& /@ responseH, wp];
		chisq =N[ 
        Plus @@ (
            Join[Map[#1^2 &, chilistA], Map[#1^2 &, chilistB], 
              Map[#1^2 &, chilistH]]), wp];
		{{chilistA, chilistB}, 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[]