(***********************************************************************
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.
***********************************************************************)

(* Tilte: 
    Nonlinear Regression Analysis for Segregation Distortion in F2 Intercross 
      by TwoZTwo model *)

(* :Context: Application`TwoZTwo`   *)

(* :Name:TwoZTwo` *)

(* :Author: Yoshiaki Harushima *)

BeginPackage["SegAnalysis`TwoZTwo`", "SegAnalysis`TwoZTwoA`", 
  "SegAnalysis`TwoZTwoB`", "SegAnalysis`TwoZTwoH`", 
  "SegAnalysis`jaco2Z2`Ax1`", "SegAnalysis`jaco2Z2`Bx1`", 
  "SegAnalysis`jaco2Z2`Hx1`", "SegAnalysis`jaco2Z2`At1`", 
  "SegAnalysis`jaco2Z2`Bt1`", "SegAnalysis`jaco2Z2`Ht1`", 
  "SegAnalysis`jaco2Z2`Ax2`", "SegAnalysis`jaco2Z2`Bx2`", 
  "SegAnalysis`jaco2Z2`Hx2`", "SegAnalysis`jaco2Z2`At2`", 
  "SegAnalysis`jaco2Z2`Bt2`", "SegAnalysis`jaco2Z2`Ht2`", 
  "SegAnalysis`jaco2Z2`Axv1`", "SegAnalysis`jaco2Z2`Bxv1`", 
  "SegAnalysis`jaco2Z2`Hxv1`", "SegAnalysis`jaco2Z2`Avb1`", 
  "SegAnalysis`jaco2Z2`Bvb1`", "SegAnalysis`jaco2Z2`Hvb1`", 
  "SegAnalysis`jaco2Z2`Avh1`", "SegAnalysis`jaco2Z2`Bvh1`", 
  "SegAnalysis`jaco2Z2`Hvh1`", "SegAnalysis`jaco2Z2`Axv2`", 
  "SegAnalysis`jaco2Z2`Bxv2`", "SegAnalysis`jaco2Z2`Hxv2`", 
  "SegAnalysis`jaco2Z2`Avb2`", "SegAnalysis`jaco2Z2`Bvb2`", 
  "SegAnalysis`jaco2Z2`Hvb2`", "SegAnalysis`jaco2Z2`Avh2`", 
  "SegAnalysis`jaco2Z2`Bvh2`", "SegAnalysis`jaco2Z2`Hvh2`",
	(* needed for graph *)
	"Graphics`MultipleListPlot`"]

(*    *)
TwoZTwo::usage = 
  "TwoZTwo[dataA, dataB, dataH, renge, guess, {opts}] searches for positions \
and strengths of the segregation distortion genes by a least-squares fitting \
to lists of data within the renge according to the modele.  The model assumes \
two gametephyte genes on the same gamete and two zygotic viability genes to \
explain the marker segregation distortions on the chromosome.  The guess \
should be {{gametephyte gene position1, transmittance of A gamete}, \
{gametephyte gene position2, transmittance of A gamete}, {gene position 1, \
relative viability of B homozygote to A homozygote, relative viability of H \
homozygote to A homozygote}, {gene position 2, relative viability of B \
homozygote to A homozygote, relative viability of H homozygote to A \
homozygote}}.  The dataA, dataB and dataH can have the form {{x1, A1}, {x2, \
A2}, ...},  {{x1, B1}, {x2, B2}, ...},  and {{x1, H1}, {x2, H2}, ...}.  "

Options[TwoZTwo] = 
  Sort[Join[
		{AccuracyGoal -> Automatic, 
     Gradient -> Automatic, 
		MaxIterations -> 30, 
     PrecisionGoal -> Automatic, 
		Weights -> Equal, 
     WorkingPrecision -> $MachinePrecision,
		FixedParameters ->{}}, 
    Options[LinearSolve]]]

Begin["`Private`"]

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

TwoZTwo::baddataA = 
  "The dataA argument of TwoZTwo must be a vector."

TwoZTwo::baddataB = 
  "The dataB argument of TwoZTwo must be a vector."

TwoZTwo::baddataH = 
  "The dataH argument of TwoZTwo must be a vector."

TwoZTwo::badrenge = 
  "The renge argument of TwoZTwo must be `1` < `2`."

TwoZTwo::degrees = 
  "Warning: The data set is smaller than the number of parameters."

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

TwoZTwo::badits =
"Warning: MaxIterations -> `1` is not set to a positive integer;
setting to 30."

TwoZTwo::badwork =
"Warning: The value `1` given for the WorkingPrecision is not an integer or \
is
an integer less than $MachinePrecision.  Setting to $MachinePrecision."

TwoZTwo::badprec =
"Warning: The value `1` given for the PrecisionGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

TwoZTwo::badacc =
"Warning: The value `1` given for the AccuracyGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

(* ===  TwoZTwo fatal error messages === *)

TwoZTwo::fitfail = "The fitting algorithm failed."

TwoZTwo::badguess = 
  "The guess should be {{gametephyte gene position\:30001, transmittance of A \
gamete}, {gametephyte gene position\:30002, transmittance of A gamete}, {gene \
position\:30001, relative viability of B homozygote to A homozygote, relative \
viability of H homozygote to A homozygote}, {gene position\:30002, relative \
viability of B homozygote to A homozygote, relative viability of H homozygote \
to A homozygote}}. Your guess is ``"

(* =========    TwoZTwo     =========== *)

(* TwoZTwo arguments;
	   dataA - matrix of numbers;
	   dataB - matrix of numbers;
	   dataH - matrix of numbers;
	   renge - the renge to fit modele equations;
	   guess - the initial guess;
	   opts - the options
		*)

TwoZTwo[dataA_List, dataB_List, dataH_List, renge_List, 
    guess_List, (opts___)?OptionQ]:=Module[{result, check}, 
    (
			result 
		)			/; (check = Check[
					 Module[{irenge, iguess}, 
					irenge = N[renge]; iguess = N[guess];
          result = oneoneztwo[dataA, dataB, dataH,irenge, iguess, opts]; 
           result =!= $Failed
						], False,TwoZTwo::fitfail];
				 check)
	] /; If[MatrixQ[N[dataA], NumberQ] || VectorQ[N[dataA], NumberQ], True, 
     Message[TwoZTwo::baddataA]; False] && 
    If[MatrixQ[N[dataB], NumberQ] || VectorQ[N[dataB], NumberQ], True, 
     Message[TwoZTwo::baddataB]; False] && 
    If[MatrixQ[N[dataH], NumberQ] || VectorQ[N[dataH], NumberQ], True, 
     Message[TwoZTwo::baddataH]; False]

oneoneztwo[dataA_, dataB_, dataH_, renge_, guess_, opts___] := 
  Module[{fixparams, maxits, workprec, accgoal, precgoal, modulus, zerotest, 
      minrenge, maxrenge, rdataA, rdataB, rdataH, ptsA, ptsB, ptsH, 
      responseA, responseB, responseH, prec, nptsA, nptsB, nptsH, x1, t1, x2, 
      t2, xv1, vb1, vh1, xv2, vb2, vh2, nresponseA, nresponseB, nresponseH, 
      nx1, nt1, nx2, nt2, nxv1, nvb1, nvh1, nxv2, nvb2, nvh2},
		(* Get options *) 
   {fixparams, maxits, workprec, accgoal, precgoal, modulus, zerotest} = 
     {FixedParameters, MaxIterations, WorkingPrecision, AccuracyGoal, 
            PrecisionGoal, Modulus, ZeroTest} /. 
				{opts} /.  Options[TwoZTwo]; 
		(* Get and check renge *)
		 {minrenge, maxrenge} = renge;
    If[!(minrenge ===-Infinity || NumberQ[minrenge]) || 
      !(maxrenge === Infinity || NumberQ[maxrenge]) || 
      minrenge >= maxrenge, 
     Message[TwoZTwo::badrenge, minrenge, maxrenge]; Return[$Failed]];
		(*  - check guess *)
		{{x1, t1}, {x2, t2}, {xv1, vb1, vh1}, {xv2, vb2, vh2}} = guess;
		If[MatrixQ[guess, !NumberQ], Message[TwoZTwo::badguess, guess]; 
      Return[$Failed]];
		(* Select data by renge *)
		rdataA = Select[dataA, minrenge < #1[[1]] < maxrenge & ];
		rdataB = Select[dataB, minrenge < #1[[1]] < maxrenge & ];
		rdataH = Select[dataH, minrenge < #1[[1]] < maxrenge & ];
		(*  - assume univariate response *)
		responseA = Last /@ rdataA;
		ptsA = Flatten[(Take[#1, 1] & ) /@ rdataA];
		responseB = Last /@ rdataB;
		ptsB = Flatten[(Take[#1, 1] & ) /@ rdataB];
		responseH = Last /@ rdataH;
		ptsH =  Flatten[(Take[#1, 1] & ) /@ rdataH];
		(* - check degrees of freedom *)
		If[ Length[ptsA] + Length[ptsB] + Length[ptsH] < 10-Length[fixparams], 
      Message[TwoZTwo::degrees]];
		(*    *)
		If[ !(IntegerQ[maxits] && Positive[maxits]), 
     Message[TwoZTwo::badits, maxits]; maxits = 30];
		If[ !(IntegerQ[workprec] && 
        workprec >= $MachinePrecision), 
     Message[TwoZTwo::badwork, workprec]; 
      workprec = $MachinePrecision];
		If[precgoal === Automatic, precgoal = workprec - 10];
		If[accgoal === Automatic, accgoal = workprec - 10];
		If[ !(IntegerQ[precgoal] && Positive[precgoal]), 
     Message[TwoZTwo::badprec, precgoal]; 
      precgoal = $MachinePrecision - 10];
		If[ !(IntegerQ[accgoal] && Positive[accgoal]), 
     Message[TwoZTwo::badacc, accgoal]; 
      accgoal = $MachinePrecision - 10];
		(* - set working precision of inputs *)
		(* Note: 
        need the precision of the argument function to be at least as large 
          as WorkingPrecision *)
		nptsA = N[ptsA, workprec];
		nresponseA = N[responseA, workprec];
		nptsB = N[ptsB, workprec];
		nresponseB = N[responseB, workprec];
		nptsH = N[ptsH, workprec];
		nresponseH = N[responseH, workprec];
		(* - call proper fitting routine *)
		nx1 = N[x1, workprec]; nt1 = N[t1, workprec];
		nx2 = N[x2, workprec]; nt2 = N[t2, workprec];
		nxv1 = N[xv1, workprec];
		nvb1 = N[vb1, workprec]; nvh1 = N[vh1, workprec];
		nxv2 = N[xv2, workprec];
		nvb2 = N[vb2, workprec]; nvh2 = N[vh2, workprec];
		out = mnlr[nptsA, nptsB, nptsH, nresponseA, nresponseB, nresponseH, nx1, 
        nt1, nx2, nt2, nxv1, nvb1, nvh1, nxv2, nvb2, nvh2, fixparams, maxits, 
        workprec, precgoal, accgoal, modulus, zerotest];
		If[out =!= $Failed, Return[MakeOutputList[out]]];
]

(* This uses Guess-Newton method for finding the least-
    squares solution for two gametophyte genes and two zygotic viability 
      genes model  *)

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, x1_, t1_, 
    x2_, t2_, xv1_, vb1_, vh1_, xv2_, vb2_, vh2_, fp_, maxits_, wp_, pg_, 
    ag_, mdls_, zt_] := 
  Module[{predA, predB, predH, chilistA, chilistB, chilistH, tmpchilistA, 
      tmpchilistB, tmpchilistH, tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, 
      chisq,derAx1, derBx1, derHx1, derx1, derAt1, derBt1, derHt1, dert1, 
      derAx2, derBx2, derHx2, derx2, derAt2, derBt2, derHt2, dert2, derAxv1, 
      derBxv1, derHxv1, derxv1, derAvb1, derBvb1, derHvb1, dervb1,derAvh1, 
      derBvh1, derHvh1, dervh1, derAxv2, derBxv2, derHxv2, derxv2, derAvb2, 
      derBvb2, derHvb2, dervb2, derAvh2, derBvh2, derHvh2, dervh2, 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, guessxv1 = xv1,
      guessvb1 = vb1, guessvh1 = vh1, guessxv2 = xv2,guessvb2 = vb2, 
      guessvh2 = vh2, oldx1, oldt1, oldx2, oldt2, oldxv1, oldvb1, oldvh1, 
      oldxv2, oldvb2, oldvh2, delta, deltax1, deltat1, deltaxv1, deltavb1, 
      deltavh1, deltaxv2, deltavb2, deltavh2, tmpx1, tmpt1, tmpx2, tmpt2, 
      tmpxv1, tmpvb1, tmpvh1, tmpxv2, tmpvb2, tmpvh2, dits, dmaxits, 
      newpredA, newpredB, newpredH, newchilistA, newchilistB, newchilistH, 
      accgoal = N[10^(-ag)], precgoal = N[10^(-pg)], paramerr, errorDOF},
		
		(* - calculation chi-squars *)
		{{chilistA, chilistB, chilistH}, chisq} = 
      calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB,  nresponseH}, {
          guessx1, guesst1, guessx2, guesst2, guessxv1,guessvb1, guessvh1, 
          guessxv2,guessvb2, guessvh2}, wp];
		tmpchi = chisq; newchi = -1; its = 1;
		oldx1 = guessx1; oldt1 = guesst1;
		oldx2 = guessx2; oldt2 = guesst2; 
		oldxv1 = guessxv1;
		oldvb1 = guessvb1; oldvh1 = guessvh1; 
		oldxv2 = guessxv2;
		oldvb2 = guessvb2; oldvh2 = guessvh2;

		(* 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;
					oldx2 = guessx2; oldt2 = guesst2; 
					oldxv1 = guessxv1; 
					oldvb1 = guessvb1; oldvh1 = guessvh1;
					oldxv2 = guessxv2; 
					oldvb2 = guessvb2; oldvh2 = guessvh2;
					
					(* === find deltax1, deltat1, deltax2, deltat2, deltaxv1, deltavb1, 
            deltavh1, deltaxv2, deltavb2, 
            and deltavh2 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]];
					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, guessxv1]];
					tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessxv1]];
					tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessxv1]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessxv1]];
					tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessxv1]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessxv1]];
					tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessxv2]];
					tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessxv2]];
					tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessxv2]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessxv2]];
					tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessxv2]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessxv2]];
					tmpchilist = Join[tmpchilistA , tmpchilistB, tmpchilistH];
					(* calculate jacobian *)
					
					der={};
					delta={};
					If[!(MemberQ[fp, "x1"]),
						derAx1 =
              SegAnalysis`jaco2Z2`Ax1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBx1 =
              SegAnalysis`jaco2Z2`Bx1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHx1 =
              SegAnalysis`jaco2Z2`Hx1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						derx1 = Join[derAx1, derBx1, derHx1];
						der=Append[der, derx1]];
					If[!(MemberQ[fp, "t1"]),
						derAt1 =
              SegAnalysis`jaco2Z2`At1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBt1 =
              SegAnalysis`jaco2Z2`Bt1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHt1 =
              SegAnalysis`jaco2Z2`Ht1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dert1 = Join[derAt1, derBt1, derHt1];
						der=Append[der, dert1]];
					If[!(MemberQ[fp, "x2"]),
						derAx2 =
              SegAnalysis`jaco2Z2`Ax2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBx2 =
              SegAnalysis`jaco2Z2`Bx2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHx2 =
              SegAnalysis`jaco2Z2`Hx2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						derx2 = Join[derAx2, derBx2, derHx2];
						der=Append[der, derx2]];
					If[!(MemberQ[fp, "t2"]),
						derAt2 =
              SegAnalysis`jaco2Z2`At2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBt2 =
              SegAnalysis`jaco2Z2`Bt2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHt2 =
              SegAnalysis`jaco2Z2`Ht2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dert2 = Join[derAt2, derBt2, derHt2];
						der=Append[der, dert2]];
					If[!(MemberQ[fp, "xv1"]),
						derAxv1 =
              SegAnalysis`jaco2Z2`Axv1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBxv1 =
              SegAnalysis`jaco2Z2`Bxv1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHxv1 =
              SegAnalysis`jaco2Z2`Hxv1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						derxv1 = Join[derAxv1, derBxv1, derHxv1];
						der=Append[der, derxv1]];
					If[!(MemberQ[fp, "vb1"]),
						derAvb1 =
              SegAnalysis`jaco2Z2`Avb1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBvb1 =
              SegAnalysis`jaco2Z2`Bvb1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHvb1 =
              SegAnalysis`jaco2Z2`Hvb1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dervb1 = Join[derAvb1, derBvb1, derHvb1];
						der=Append[der, dervb1]];
					If[!(MemberQ[fp, "vh1"]),
						derAvh1 =
              SegAnalysis`jaco2Z2`Avh1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBvh1 =
              SegAnalysis`jaco2Z2`Bvh1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHvh1 =
              SegAnalysis`jaco2Z2`Hvh1`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dervh1 = Join[derAvh1, derBvh1, derHvh1];
						der=Append[der, dervh1]];
					If[!(MemberQ[fp, "xv2"]),
						derAxv2 =
              SegAnalysis`jaco2Z2`Axv2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBxv2 =
              SegAnalysis`jaco2Z2`Bxv2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHxv2 =
              SegAnalysis`jaco2Z2`Hxv2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						derxv2 = Join[derAxv2, derBxv2, derHxv2];
						der=Append[der, derxv2]];
					If[!(MemberQ[fp, "vb2"]),
						derAvb2 =
              SegAnalysis`jaco2Z2`Avb2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBvb2 =
              SegAnalysis`jaco2Z2`Bvb2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHvb2 =
              SegAnalysis`jaco2Z2`Hvb2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dervb2 = Join[derAvb2, derBvb2, derHvb2];
						der=Append[der, dervb2]];
					If[!(MemberQ[fp, "vh2"]),
						derAvh2 =
              SegAnalysis`jaco2Z2`Avh2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsA, wp];
						derBvh2 =
              SegAnalysis`jaco2Z2`Bvh2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsB, wp];
						derHvh2 =
              SegAnalysis`jaco2Z2`Hvh2`Private`der[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
                guessvh2, tmpnptsH, wp];
						dervh2 = Join[derAvh2, derBvh2, derHvh2];
						der=Append[der, dervh2]];
					
					
					(* 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, "xv1"]),deltaxv1=First[delta];delta=Rest[delta], 
            deltaxv1=0];
					If[!(MemberQ[fp, "vb1"]),deltavb1=First[delta]; delta=Rest[delta],  
            deltavb1=0];
					If[!(MemberQ[fp, "vh1"]),deltavh1=First[delta];delta=Rest[delta],  
            deltavh1=0];
					If[!(MemberQ[fp, "xv2"]),deltaxv2=First[delta];delta=Rest[delta], 
            deltaxv2=0];
					If[!(MemberQ[fp, "vb2"]),deltavb2=First[delta]; delta=Rest[delta],  
            deltavb2=0];
					If[!(MemberQ[fp, "vh2"]),deltavh2=First[delta],  deltavh2=0];
					
					{tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, 
              tmpvh2} = 
            N[{oldx1, oldt1, oldx2, oldt2, oldxv1, oldvb1, oldvh1, oldxv2, 
                  oldvb2, oldvh2} +  {deltax1, deltat1, deltax2, deltat2, 
                  deltaxv1, deltavb1, deltavh1, deltaxv2, deltavb2, 
                  deltavh2}, wp];
					
					{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                nresponseH}, {tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, 
                tmpvh1, tmpxv2,tmpvb2, tmpvh2}, wp];
					
					(* find newguess for smaller chisq  *)
					dits = 1; dmaxits = 15;
					While[newchi - chisq >= 0 && dits <= dmaxits,
						{tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, 
                tmpvh2} = 
              N[{oldx1, oldt1, oldx2, oldt2, oldxv1, oldvb1, oldvh1, oldxv2, 
                    oldvb2, 
                    oldvh2} + (2^(-1*dits))*{deltax1, deltat1, deltax2, 
                      deltat2, deltaxv1, deltavb1, deltavh1, deltaxv2, 
                      deltavb2, deltavh2}, wp];
						{{newchilistA, newchilistB, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                  nresponseH}, {tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, 
                  tmpvh1, tmpxv2,tmpvb2, tmpvh2}, wp];
						dits++];
					If[dits > dmaxits, Message[TwoZTwo::lmincrease, dmaxits]];
					tmpchi = newchi;
					{guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, guessvh1, 
              guessxv2, guessvb2, guessvh2} = {tmpx1, tmpt1, tmpx2, tmpt2, 
              tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2};
					chilistA = newchilistA;
					chilistB = newchilistB;
					chilistH = newchilistH;
					printprogress[its, 
            newchi, {guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, 
              guessvh1, guessxv2, guessvb2, guessvh2}];
					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[chilistH] - 10+Length[fp];
		paramerr = Sqrt[ ( newchi/errorDOF) DiagonalElements[Inverse[alpha]]];
		{newchi, {guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, guessvh1, 
        guessxv2, guessvb2, guessvh2}, paramerr, errorDOF, fp}
	]  (* end TwoZTwo *)


(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2, ptsA, ptsB, ptsH, 
      responseA, responseB, responseH, predA, predB, predH, chilistA, 
      chilistB, chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{ptsA, ptsB, ptsH} = pts;
		{x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2} =N[guess,wp];
		predA = 
      N[TwoZTwoA[#1, x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2] & /@ ptsA, 
        wp];
		predB = 
      N[TwoZTwoB[#1, x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2] & /@ ptsB, 
        wp];
		predH = 
      N[TwoZTwoH[#1, x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2] & /@ 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}
		]

(* ============  MakeOutPutList  ============== *)

MakeOutputList[out_] := 
  Block[
		{x, chisq, param, x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2, paramerr, 
      errx1, errt1, errx2, errt2, errxv1, errvb1, errvh1, errxv2, errvb2, 
      errvh2, errorDOF},
		{chisq, param, paramerr, errorDOF, fp} = out;
		{x1, t1, x2, t2, xv1, vb1, vh1, xv2, vb2, vh2} = param;
		If[!(MemberQ[fp, "x1"]),errx1=First[paramerr];
				paramerr=Rest[paramerr], errx1=0];
		If[!(MemberQ[fp, "t1"]),errt1=First[paramerr];paramerr=Rest[paramerr], 
      errt1=0];
		If[!(MemberQ[fp, "x2"]),errx2=First[paramerr];
				paramerr=Rest[paramerr], errx2=0];
		If[!(MemberQ[fp, "t2"]),errt2=First[paramerr];paramerr=Rest[paramerr], 
      errt2=0];
		If[!(MemberQ[fp, "xv1"]),errxv1=First[paramerr];paramerr=Rest[paramerr], 
      errxv1=0];
		If[!(MemberQ[fp, "vb1"]),errvb1=First[paramerr];
				paramerr=Rest[paramerr], errvb1=0];
		If[!(MemberQ[fp, "vh1"]),errvh1=First[paramerr];paramerr=Rest[paramerr], 
      errvh1=0];
		If[!(MemberQ[fp, "xv2"]),errxv2=First[paramerr];paramerr=Rest[paramerr], 
      errxv2=0];
		If[!(MemberQ[fp, "vb2"]),errvb2=First[paramerr];
				paramerr=Rest[paramerr], errvb2=0];
		If[!(MemberQ[fp, "vh2"]),errvh2=First[paramerr], errvh2=0];
		Print["The model is two gametophyte genes on the same gamete and two \
zygotic viability genes."];
		If[fp=!={},Print["The fixed parameters are ",fp]];
		Print[StringForm[
        "One gametophyte gene position is at `` \[PlusMinus] `` .", 
        NumberForm[x1, 6], NumberForm[errx1, 6]]];
		Print[StringForm[
        "The transmittance of the A gamete by the gene is `` \[PlusMinus] `` \
.", NumberForm[t1, 6], NumberForm[errt1, 6]]];
		Print[StringForm[
        "Another gametophyte gene position is at `` \[PlusMinus] `` .", 
        NumberForm[x2, 6], NumberForm[errx2, 6]]];
		Print[StringForm[
        "The transmittance of the A gamete by the gene is `` \[PlusMinus] `` \
.", NumberForm[t2, 6], NumberForm[errt2, 6]]];
		Print[StringForm[
        "One zygotic viability gene position is at `` \[PlusMinus] `` .", 
        NumberForm[xv1, 6], NumberForm[errxv1, 6]]];
		Print[StringForm[
        "The relative viability of B homozygote is `` \[PlusMinus] `` .", 
        NumberForm[vb1, 6], NumberForm[errvb1, 6]]];
		Print[StringForm[
        "The relative viability of heterozygote is `` \[PlusMinus] `` .", 
        NumberForm[vh1, 6], NumberForm[errvh1, 6]]];
		Print[StringForm[
        "Another zygotic viability gene position is at `` \[PlusMinus] `` .", 
        NumberForm[xv2, 6], NumberForm[errxv2, 6]]];
		Print[StringForm[
        "The relative viability of B homozygote is `` \[PlusMinus] `` .", 
        NumberForm[vb2, 6], NumberForm[errvb2, 6]]];
		Print[StringForm[
        "The relative viability of heterozygote is `` \[PlusMinus] `` .", 
        NumberForm[vh2, 6], NumberForm[errvh2, 6]]];
		Print[StringForm[
        "The number of analysed data is `1`.  The root mean square is `2`.", 
        NumberForm[errorDOF + 10-Length[fp], 6], 
        NumberForm[Sqrt[chisq/errorDOF], 6]]]
	]

(*   *)
printprogress[iteration_, chisq_, params_] := 
     Print[
    StringForm[
      "Iteration:`1` ChiSquared:`2`  {x1, t1, x2, t2, xv1, vb1, vh1, xv2, \
vb2, vh2}:`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[]