(***********************************************************************
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 ZThree model *)

(* :Context: Application`ZThree`   *)

(* :Name:ZThree` *)

(* :Author: Yoshiaki Harushima *)

BeginPackage["SegAnalysis`ZThree`", "SegAnalysis`ZThreeA`", 
  "SegAnalysis`ZThreeB`", "SegAnalysis`ZThreeH`", "SegAnalysis`jacoZ3`Axv1`", 
  "SegAnalysis`jacoZ3`Bxv1`", "SegAnalysis`jacoZ3`Hxv1`", 
  "SegAnalysis`jacoZ3`Avb1`", "SegAnalysis`jacoZ3`Bvb1`", 
  "SegAnalysis`jacoZ3`Hvb1`", "SegAnalysis`jacoZ3`Avh1`", 
  "SegAnalysis`jacoZ3`Bvh1`", "SegAnalysis`jacoZ3`Hvh1`", 
  "SegAnalysis`jacoZ3`Axv2`", "SegAnalysis`jacoZ3`Bxv2`", 
  "SegAnalysis`jacoZ3`Hxv2`", "SegAnalysis`jacoZ3`Avb2`", 
  "SegAnalysis`jacoZ3`Bvb2`", "SegAnalysis`jacoZ3`Hvb2`", 
  "SegAnalysis`jacoZ3`Avh2`", "SegAnalysis`jacoZ3`Bvh2`", 
  "SegAnalysis`jacoZ3`Hvh2`", "SegAnalysis`jacoZ3`Axv3`", 
  "SegAnalysis`jacoZ3`Bxv3`", "SegAnalysis`jacoZ3`Hxv3`", 
  "SegAnalysis`jacoZ3`Avb3`", "SegAnalysis`jacoZ3`Bvb3`", 
  "SegAnalysis`jacoZ3`Hvb3`", "SegAnalysis`jacoZ3`Avh3`", 
  "SegAnalysis`jacoZ3`Bvh3`", "SegAnalysis`jacoZ3`Hvh3`",
	(* needed for graph *)
	"Graphics`MultipleListPlot`"]

(*    *)
ZThree::usage = 
  "ZThree[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 \
three zygotic viability genes to explain the marker segregation distortions \
on the chromosome.  The guess should be {{the first gene position, the \
relative viability of B homozygote to A homozygote, the relative viability of \
H homozygote to A homozygote}, {the second gene position, the relative \
viability of B homozygote to A homozygote, the relative viability of H \
homozygote to A homozygote}, {the third gene position, the relative viability \
of B homozygote to A homozygote, the 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[ZThree] = 
  Sort[Join[
		{AccuracyGoal -> Automatic, 
     Gradient -> Automatic, 
		MaxIterations -> 30, 
     PrecisionGoal -> Automatic, 
		Weights -> Equal, 
     WorkingPrecision -> $MachinePrecision,
		FixedParameters ->{}}, 
    Options[LinearSolve]]]

Begin["`Private`"]

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

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

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

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

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

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

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

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

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

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

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

(* ===  ZThree fatal error messages === *)

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

ZThree::badguess = 
  "The guess should be {{the first gene position, the relative viability of B \
homozygote to A homozygote, the relative viability of H homozygote to A \
homozygote}, {the second gene position, the relative viability of B \
homozygote to A homozygote, the relative viability of H homozygote to A \
homozygote}, {the third gene position, the relative viability of B homozygote \
to A homozygote, the relative viability of H homozygote to A homozygote}}. \
Your guess is ``"

(* =========    ZThree     =========== *)

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

ZThree[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 = zthree[dataA, dataB, dataH,irenge, iguess, opts]; 
           result =!= $Failed
						], False,ZThree::fitfail];
				 check)
	] /; If[MatrixQ[N[dataA], NumberQ] || VectorQ[N[dataA], NumberQ], True, 
     Message[ZThree::baddataA]; False] && 
    If[MatrixQ[N[dataB], NumberQ] || VectorQ[N[dataB], NumberQ], True, 
     Message[ZThree::baddataB]; False] && 
    If[MatrixQ[N[dataH], NumberQ] || VectorQ[N[dataH], NumberQ], True, 
     Message[ZThree::baddataH]; False]

zthree[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, xv1, vb1, 
      vh1, xv2, vb2, vh2, xv3, vb3, vh3, nresponseA, nresponseB, nresponseH, 
      nxv1, nvb1, nvh1, nxv2, nvb2, nvh2, nxv3, nvb3, nvh3},
		(* Get options *) 
   {fixparams, maxits, workprec, accgoal, precgoal, modulus, zerotest} = 
     {FixedParameters, MaxIterations, WorkingPrecision, AccuracyGoal, 
            PrecisionGoal, Modulus, ZeroTest} /. 
				{opts} /.  Options[ZThree]; 
		(* Get and check renge *)
		 {minrenge, maxrenge} = renge;
    If[!(minrenge ===-Infinity || NumberQ[minrenge]) || 
      !(maxrenge === Infinity || NumberQ[maxrenge]) || 
      minrenge >= maxrenge, 
     Message[ZThree::badrenge, minrenge, maxrenge]; Return[$Failed]];
		(*  - check guess *)
		{{xv1, vb1, vh1}, {xv2, vb2, vh2}, {xv3, vb3, vh3}} = guess;
		If[MatrixQ[guess, !NumberQ], Message[ZThree::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] < 9-Length[fixparams], 
      Message[ZThree::degrees]];
		(*    *)
		If[ !(IntegerQ[maxits] && Positive[maxits]), 
     Message[ZThree::badits, maxits]; maxits = 30];
		If[ !(IntegerQ[workprec] && 
        workprec >= $MachinePrecision), 
     Message[ZThree::badwork, workprec]; 
      workprec = $MachinePrecision];
		If[precgoal === Automatic, precgoal = workprec - 10];
		If[accgoal === Automatic, accgoal = workprec - 10];
		If[ !(IntegerQ[precgoal] && Positive[precgoal]), 
     Message[ZThree::badprec, precgoal]; 
      precgoal = $MachinePrecision - 10];
		If[ !(IntegerQ[accgoal] && Positive[accgoal]), 
     Message[ZThree::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];
		nxv1 = N[xv1, workprec];
		nvb1 = N[vb1, workprec]; nvh1 = N[vh1, workprec];
		nxv2 = N[xv2, workprec];
		nvb2 = N[vb2, workprec]; nvh2 = N[vh2, workprec];
		nxv3 = N[xv3, workprec];
		nvb3 = N[vb3, workprec]; nvh3 = N[vh3, workprec];
		out = mnlr[nptsA, nptsB, nptsH, nresponseA, nresponseB, nresponseH, nxv1, 
        nvb1, nvh1, nxv2, nvb2, nvh2, nxv3, nvb3, nvh3, fixparams, maxits, 
        workprec, precgoal, accgoal, modulus, zerotest];
		If[out =!= $Failed, Return[MakeOutputList[out]]];
]

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

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, xv1_, 
    vb1_, vh1_, xv2_, vb2_, vh2_, xv3_, vb3_, vh3_, fp_, maxits_, wp_, pg_, 
    ag_, mdls_, zt_] := 
  Module[{predA, predB, predH, chilistA, chilistB, chilistH, tmpchilistA, 
      tmpchilistB, tmpchilistH, tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, 
      chisq, derAxv1, derBxv1, derHxv1, derxv1, derAvb1, derBvb1, derHvb1, 
      dervb1, derAvh1, derBvh1, derHvh1, dervh1, derAxv2, derBxv2, derHxv2, 
      derxv2, derAvb2, derBvb2, derHvb2, dervb2, derAvh2, derBvh2, derHvh2, 
      dervh2, derAxv3, derBxv3, derHxv3, derxv3, derAvb3, derBvb3, derHvb3, 
      dervb3, derAvh3, derBvh3, derHvh3, dervh3, der, alpha, beta, chisq0, 
      newchi, tmpchi, its, accDelta, NAcc, precDelta, NPrec, invaliddigits, 
      accflag0, accflag1 = False, accflag2 = False, precflag0, 
      precflag1 = False, precflag2 = False, ovflag = False, guessxv1 = xv1,
      guessvb1 = vb1, guessvh1 = vh1, guessxv2 = xv2,guessvb2 = vb2, 
      guessvh2 = vh2, guessxv3 = xv3,guessvb3 = vb3, guessvh3 = vh3, oldxv1, 
      oldvb1, oldvh1, oldxv2, oldvb2, oldvh2, oldxv3, oldvb3, oldvh3, delta, 
      deltaxv1, deltavb1, deltavh1, deltaxv2, deltavb2, deltavh2, deltaxv3, 
      deltavb3, deltavh3, tmpxv1, tmpvb1, tmpvh1, tmpxv2, tmpvb2, tmpvh2, 
      tmpxv3, tmpvb3, tmpvh3, 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}, {
          guessxv1,guessvb1, guessvh1, guessxv2,guessvb2, guessvh2, guessxv3,
          guessvb3, guessvh3}, wp];
		tmpchi = chisq; newchi = -1; its = 1; 
		oldxv1 = guessxv1;
		oldvb1 = guessvb1; oldvh1 = guessvh1; 
		oldxv2 = guessxv2;
		oldvb2 = guessvb2; oldvh2 = guessvh2; 
		oldxv3 = guessxv3;
		oldvb3 = guessvb3; oldvh3 = guessvh3;

		(* 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;
					oldxv1 = guessxv1; 
					oldvb1 = guessvb1; oldvh1 = guessvh1;
					oldxv2 = guessxv2; 
					oldvb2 = guessvb2; oldvh2 = guessvh2;
					oldxv3 = guessxv3; 
					oldvb3 = guessvb3; oldvh3 = guessvh3;
					
					(* === find deltaxv1, deltavb1, deltavh1, deltaxv2, deltavb2, deltavh2, 
            deltaxv3, 
            deltavb3 and deltavh3 by LinearSolve normalmatrix ==== *)
					
					(* === delete not differentiable position ==== *)
					tmpchilistA = Delete[chilistA, Position[nptsA, guessxv1]];
					tmpnptsA =Delete[nptsA, Position[nptsA, guessxv1]];
					tmpchilistB = Delete[chilistB, Position[nptsB, guessxv1]];
					tmpnptsB =Delete[nptsB, Position[nptsB, guessxv1]];
					tmpchilistH = Delete[chilistH, Position[nptsH, guessxv1]];
					tmpnptsH =Delete[nptsH, Position[nptsH, 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]];
					tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessxv3]];
					tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessxv3]];
					tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessxv3]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessxv3]];
					tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessxv3]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessxv3]];
					tmpchilist = Join[tmpchilistA , tmpchilistB, tmpchilistH];
					(* calculate jacobian *)
					
					der={};
					delta={};
					If[!(MemberQ[fp, "xv1"]),
						derAxv1 =
              SegAnalysis`jacoZ3`Axv1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBxv1 =
              SegAnalysis`jacoZ3`Bxv1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp] ;
						derHxv1 =
              SegAnalysis`jacoZ3`Hxv1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						derxv1 = Join[derAxv1, derBxv1, derHxv1];
						der=Append[der, derxv1]];
					If[!(MemberQ[fp, "vb1"]),
						derAvb1 =
              SegAnalysis`jacoZ3`Avb1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvb1 =
              SegAnalysis`jacoZ3`Bvb1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvb1 =
              SegAnalysis`jacoZ3`Hvb1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						dervb1 = Join[derAvb1, derBvb1, derHvb1];
						der=Append[der, dervb1]];
					If[!(MemberQ[fp, "vh1"]),
						derAvh1 =
              SegAnalysis`jacoZ3`Avh1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvh1 =
              SegAnalysis`jacoZ3`Bvh1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvh1 =
              SegAnalysis`jacoZ3`Hvh1`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						dervh1 = Join[derAvh1, derBvh1, derHvh1];
						der=Append[der, dervh1]];
					If[!(MemberQ[fp, "xv2"]),
						derAxv2 =
              SegAnalysis`jacoZ3`Axv2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBxv2 =
              SegAnalysis`jacoZ3`Bxv2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHxv2 =
              SegAnalysis`jacoZ3`Hxv2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						derxv2 = Join[derAxv2, derBxv2, derHxv2];
						der=Append[der, derxv2]];
					If[!(MemberQ[fp, "vb2"]),
						derAvb2 =SegAnalysis`jacoZ3`Avb2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvb2 =
              SegAnalysis`jacoZ3`Bvb2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvb2 =
              SegAnalysis`jacoZ3`Hvb2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp] ;
						dervb2 = Join[derAvb2, derBvb2, derHvb2];
						der=Append[der, dervb2]];
					If[!(MemberQ[fp, "vh2"]),
						derAvh2 =
              SegAnalysis`jacoZ3`Avh2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvh2 =
              SegAnalysis`jacoZ3`Bvh2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvh2 =
              SegAnalysis`jacoZ3`Hvh2`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						dervh2 = Join[derAvh2, derBvh2, derHvh2];
						der=Append[der, dervh2]];
					If[!(MemberQ[fp, "xv3"]),
						derAxv3 =
              SegAnalysis`jacoZ3`Axv3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp] ;
						derBxv3 =
              SegAnalysis`jacoZ3`Bxv3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHxv3 =
              SegAnalysis`jacoZ3`Hxv3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						derxv3 = Join[derAxv3, derBxv3, derHxv3];
						der=Append[der, derxv3]];
					If[!(MemberQ[fp, "vb3"]),
						derAvb3 =
              SegAnalysis`jacoZ3`Avb3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvb3 =
              SegAnalysis`jacoZ3`Bvb3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvb3 =
              SegAnalysis`jacoZ3`Hvb3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						dervb3 = Join[derAvb3, derBvb3, derHvb3];
						der=Append[der, dervb3]];
					If[!(MemberQ[fp, "vh3"]),
						derAvh3 =
              SegAnalysis`jacoZ3`Avh3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsA, wp];
						derBvh3 =
              SegAnalysis`jacoZ3`Bvh3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsB, wp];
						derHvh3 =
              SegAnalysis`jacoZ3`Hvh3`Private`der[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, guessvb3, 
                guessvh3, tmpnptsH, wp];
						dervh3 = Join[derAvh3, derBvh3, derHvh3];
						der=Append[der, dervh3]];
				
					
					(* 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, "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];delta=Rest[delta],  
            deltavh2=0];
					If[!(MemberQ[fp, "xv3"]),deltaxv3=First[delta];delta=Rest[delta], 
            deltaxv3=0];
					If[!(MemberQ[fp, "vb3"]),deltavb3=First[delta]; delta=Rest[delta],  
            deltavb3=0];
					If[!(MemberQ[fp, "vh3"]),deltavh3=First[delta],  deltavh3=0];
					
					{tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2, tmpxv3,tmpvb3, tmpvh3} = 
            N[{oldxv1, oldvb1, oldvh1, oldxv2, oldvb2, oldvh2, oldxv3, 
                  oldvb3, oldvh3} +  {deltaxv1, deltavb1, deltavh1, deltaxv2, 
                  deltavb2, deltavh2, deltaxv3, deltavb3, deltavh3}, wp];
					
					{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                nresponseH}, {tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2, 
                tmpxv3,tmpvb3, tmpvh3}, wp];
					
					(* find newguess for smaller chisq  *)
					dits = 1; dmaxits = 15;
					While[newchi - chisq >= 0 && dits <= dmaxits,
						{tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2, tmpxv3,tmpvb3, tmpvh3} = 
              N[{oldxv1, oldvb1, oldvh1, oldxv2, oldvb2, oldvh2, oldxv3, 
                    oldvb3, 
                    oldvh3} + (2^(-1*dits))*{deltaxv1, deltavb1, deltavh1, 
                      deltaxv2, deltavb2, deltavh2, deltaxv3, deltavb3, 
                      deltavh3}, wp];
						{{newchilistA, newchilistB, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                  nresponseH}, {tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2, 
                  tmpxv3,tmpvb3, tmpvh3}, wp];
						dits++];
					If[dits > dmaxits, Message[ZThree::lmincrease, dmaxits]];
					tmpchi = newchi;
					{guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, guessvh2, guessxv3, 
              guessvb3, guessvh3} = {tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, 
              tmpvh2, tmpxv3,tmpvb3, tmpvh3};
					chilistA = newchilistA;
					chilistB = newchilistB;
					chilistH = newchilistH;
					printprogress[its, 
            newchi, {guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
              guessvh2, guessxv3, guessvb3, guessvh3}];
					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] - 9+Length[fp];
		paramerr = Sqrt[ ( newchi/errorDOF) DiagonalElements[Inverse[alpha]]];
		{newchi, {guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, guessvh2, 
        guessxv3, guessvb3, guessvh3}, paramerr, errorDOF, fp}
	]  (* end ZThree *)

(*   calculation of chi-list and chi-squares *)

calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3, ptsA, ptsB, ptsH, 
      responseA, responseB, responseH, predA, predB, predH, chilistA, 
      chilistB, chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{ptsA, ptsB, ptsH} = pts;
		{xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3} =N[guess,wp];
		predA = 
      N[ZThreeA[#1, xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3] & /@ ptsA, 
        wp];
		predB = 
      N[ZThreeB[#1, xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3] & /@ ptsB, 
        wp];
		predH = 
      N[ZThreeH[#1, xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3] & /@ 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, xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3, paramerr, 
      errxv1, errvb1, errvh1, errxv2, errvb2, errvh2, errxv3, errvb3, errvh3, 
      errorDOF},
		{chisq, param, paramerr, errorDOF, fp} = out;
		{xv1, vb1, vh1, xv2, vb2, vh2, xv3, vb3, vh3} = param;
		
		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];paramerr=Rest[paramerr], 
      errvh2=0];
		If[!(MemberQ[fp, "xv3"]),errxv3=First[paramerr];paramerr=Rest[paramerr], 
      errxv3=0];
		If[!(MemberQ[fp, "vb3"]),errvb3=First[paramerr];
				paramerr=Rest[paramerr], errvb3=0];
		If[!(MemberQ[fp, "vh3"]),errvh3=First[paramerr], errvh3=0];
		Print["The model is three zygotic viability genes."];
		If[fp=!={},Print["The fixed parameters are ",fp]];
		Print[StringForm[
        "The first 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[
        "The second 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 third zygotic viability gene position is at `` \[PlusMinus] `` \
.", NumberForm[xv3, 6], NumberForm[errxv3, 6]]];
		Print[StringForm[
        "The relative viability of B homozygote is `` \[PlusMinus] `` .", 
        NumberForm[vb3, 6], NumberForm[errvb3, 6]]];
		Print[StringForm[
        "The relative viability of heterozygote is `` \[PlusMinus] `` .", 
        NumberForm[vh3, 6], NumberForm[errvh3, 6]]];
		Print[StringForm[
        "The number of analysed data is `1`.  The root mean square is `2`.", 
        NumberForm[errorDOF + 9-Length[fp], 6], 
        NumberForm[Sqrt[chisq/errorDOF], 6]]]
	]

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