;;; Carlo C. Maley   April 12, 95

;;; This is a simple model of biodiversity.


;;; *** update env
;;; *** update orgs
;;; *** save state/stats
;;; *** output organisms.

;;; *** I have to decide if I'm going to organize live organisms into
;;; stacks of different amounts of generalism (superimpossed on the
;;; array), or if I'm going to just loosely structure the array.  The
;;; first provides more flexibility in space allocation, but enforces
;;; the eating priorities rigidly.  Maybe that would make it easier to
;;; analyze. However, since I only have 20 live organisms in each
;;; location at the moment, it hardly seems worth it.

;;; *** Like noah's ark, I need to put pairs of organisms into the
;;; initial environment.
;;; *** I have to handle plant life cycles (feeding and reproduction)
;;; *** mate finding should be restricted to your type (plant or
;;; animal) and perhaps even generalism level-1.

;;;Debug:
;;; start-placement
;;; find-mate-in-loc
;;; reproduce


(in-package 'USER)
(setf *print-array* t)



;(load "~/lib/tools")
;(cd "~/diversity")
;(load "macros")
;(load "parameters")
;(load "poisson")
;(load "stats")
;(load "debug")


;--------------------------Structures-----------------------
;;; This is the data structure for an environment location.

(defstruct location
  (orgs (make-array CAPACITY :element-type 'organism)) ; The array of
						       ; organisms. 
  (num-orgs 0) ; The number of non-food organisms alive in this location.
  (num-food 0) ; The number of basic food organisms added each time
	       ; step.  This is never changed after initialization.
  (barrier BASE-TERRAIN) ; The chance an organism cannot enter this
			 ; location (used to simulate topology).
  (climate 0) ; A bit pattern to abstract climate and types of
	      ; resources for matching with the autotroph's prey patterns.
  (x 0) ; The x-coordinate.
  (y 0) ; The y-coordinate.
  (neighborhood nil)) ; A list of pointers to the neighbor locations
		      ; in the environment.  This isn't altered after
		      ; initialization. 


;;; This is the environment structure of an organism.

(defstruct organism
  (alive nil)
  (prey 0)
  (generalism 0)
  (phenotype 0)
  (photosynth 0) ; only for plants.
  (last-meal 0)
  (num-meals 0))


; ----------------------Global Variables------------------------------

;;; A temporary variable used to construct a new organism
(defvar *fetus* (make-organism))

;;; An array of probability values for the poisson distribution.
(defvar *probabilities* (make-array GENESIZE
				  :initial-element 0.0))


;;; The environment.
(defvar *env* (make-array (list XDIM YDIM)))

;;; Engergy units (meals) are given to the autotrophs each time step.
;;; This may change with seasons, and also with match to the climate
;;; of the environment location.  To determine the number of energy
;;; units a plant gets we do
;;; (round (* ENERGY-CONVERSION (match plant-prey/generalism climate))
;;;           (* *productivity-reduction* GENESIZE))
;;;
;;; So *productivity-reduction* should have a minimum of 1 and a maximum of
;;; ENERGY-CONVERSION.

(defvar *productivity-reduction* 1)



;(defvar *plant-phenos* ()) ;an array of the phenotypes of the
			   ;autotrophs.


(defvar *clock* 0)


; ========================Main======================================

(defun run-model (num-trials max-time &optional (run-num -1))
  (when (>= run-num 0) (setf num-trials 1))
  (init-random run-num)
  (initialize-environment)
  (dotimes (trial num-trials)
     (reset-model) ; clock, biota, topology.
     (run-trial max-time)))

(defun run-trial (max-time)
  (dotimes (i max-time)
     (incf *clock*) ; increment clock
     (update-env)
     (update-orgs)
     (save-state))) ; in stats.lisp


(defun alt-setup ()
  (init-random -1)
  (initialize-environment)
  (reset-model)
  ;;; Shelley's code for counting species goes here ***
  )
; ------------------Initialization-------------------------

;;; Intialize-environment first allocates all the structures in the
;;; environment (locations and organisms) and then intializes the
;;; topology (barriers), neighborhoods, and starting organisms.

(defun initialize-environment ()
  (format t "I am starting to intialize the environment.~%")
  (init-poisson)
  (allocate-space)
  (setup-neighborhoods)
  (format t "I have finished intializing the environment.~%"))

(defun reset-model ()
  (setf *clock* 0)
  (init-topology)  
  (init-food-sources)
  (setup-population))
 
;;; This allocates the location and organism structures once and for
;;; all.

(defun allocate-space ()
  (dotimes (y YDIM)
	   (dotimes (x XDIM)
		    (let ((new-loc (make-location)))
		      (setf (location-x new-loc) x)
		      (setf (location-y new-loc) y)
		      (setf (aref *env* x y) new-loc))
		    (let ((orgs (location-orgs (aref *env* x y))))
		      (dotimes (org-ptr CAPACITY)
			       (setf (svref orgs org-ptr) (make-organism)))))))




(defun setup-neighborhoods ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(setf (location-neighborhood (aref *env* x y))
	      (remove-if #'null
			 (list (aref *env* x y)
			       (when (> x 0) (aref *env* (1- x) y))
			       (when (< x XBORDER) (aref *env* (1+ x) y))
			       (when (> y 0) (aref *env* x (1- y)))
			       (when (< y YBORDER) (aref *env* x (1+ y)))
			       (when (and (> x 0) (> y 0))
				     (aref *env* (1- x) (1- y)))
			       (when (and (< x XBORDER) (> y 0))
				     (aref *env* (1+ x) (1- y)))
			       (when (and (> x 0) (< y YBORDER))
				     (aref *env* (1- x) (1+ y)))
			       (when (and (< x XBORDER) (< y YBORDER))
				     (aref *env* (1+ x) (1+ y)))))))))



;;; This prints out some aspect of an environment.  Note: do not try
;;; to print out the structure for an environment location (by
;;; printing out the neighborhood or some such silliness).

(defun print-env (proc)
  (dotimes (y YDIM)
     (format t "~%")
     (dotimes (x XDIM)
	      (format t "~a " (funcall proc (aref *env* x y)))))
     (format t "~%"))

(defun print-loc (loc)
  (format t " Location [~a, ~a]~%" (location-x loc) (location-y loc))
  (format t "   Number of organisms: ~a~%" (location-num-orgs loc))
  (format t "   Number of food orgs: ~a~%" (location-num-food loc))
  (format t "   Barrier: ~a~%" (location-barrier loc)))
  

;;; This sets up mountains on the even rows and columns.  It also
;;; initializes the climate of a location to be a slight modification
;;; of the barrier.

(defun init-topology ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(let ((loc (aref *env* x y)))    
	  (when (or (evenp x) (evenp y))
		(setf (location-barrier loc)
		      (round (* 0.98 MAX-TERRAIN))))
	  (setf (location-climate loc)
		(mutate-bits (location-barrier loc)
			     (random MAX-CLIMATE-DIFF)
			     (if *increasing-adaptive-space*
				 GENESIZE ADAPTZONE)))))))


;;; Patchiness emerges out of the plant's prey/generalism genes
;;; matching the climate of the location.  Initially they are all
;;; random.  Except that I put in NOAH-FACTOR (3 in my case) copies of
;;; each genotype so that they can find a mate for reproduction.
;;; Otherwise no one would be able to reproduce.  The mates are the
;;; same except for the number and time of their last meal.

(defun init-food-sources ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(let* ((loc (aref *env* x y))
	       (orgs (location-orgs loc)))
	  (dotimes (i MAX-FOOD-CAPACITY)
	     (let ((org (svref orgs i)))
	       (if (= (rem i NOAH-FACTOR) 0)
		   (progn 
		     (init-org org)
		     (setf (organism-photosynth org)
			   (* (match (logior (organism-prey org)
					     (organism-generalism org))
				     (logior (organism-generalism org)
					     (location-climate loc)))
			      ENERGY-CONVERSION)))
		 (progn ; mates for the previous org.
		   (copy-org i (* NOAH-FACTOR         ; a macro
				  (floor i NOAH-FACTOR))
			     orgs)
		   (setf (organism-last-meal org)
			 (- *clock* (random MAX-FAST)))
		   (setf (organism-num-meals org) (random REPRO-THRESHOLD))))))
	  (setf (location-num-food loc) MAX-FOOD-CAPACITY)))))



;;; I would like to make this patchy.  Autotrophs ought to have
;;; phenotypes across the entire spectrum.  Seconly, the autotrophs
;;; should be optimally mixed, so that restrictions on the
;;; replinishment of the food sources (modelling seasons) doesn't
;;; dramatically favor one plant over another.  Or perhaps it should,
;;; simulating different plant tolerances for seasons.  The extent of
;;; the reduction of a plant would depend on its mostly random
;;; placement in the location.
  
;  (setf *plant-phenos* (make-array *num-autotrophs*))
;  (dotimes (i *num-autotrophs*)
;     (setf (svref *plant-phenos* i) (random (expt 2 ADAPTZONE))))
;  (dotimes (j *num-colonies*)
;     (dotimes (i *num-autotrophs*)
;          (plant-autotroph (svref *plant-phenos* i))))
;  (mix-seeds))

;;; An alternative to simple random phenotype autotrophs would be
;;; plants evenly spaced across hamming space.  Say we want 2^n
;;; autotrophs.  Then we could split up the 29 bits into n sections.
;;; Each section would either be all 1's or all 0's.  Constructing all
;;; combinations gives you 2^n autotrophs with each autotroph at least
;;; n bits away in hamming space.  This would have the effect that
;;; only an n+ bit generalist could feed on more than 1 autotroph.

;;; Note also that plants never have a pheontype outside the adaptive
;;; zone.  Only predators (evolving orgs) of one sort or another can
;;; explore that zone.

(defun plant-autotroph (pheno)
  (let* ((amount (1+ (random *patch-quantity*)))
	 (xloc (- (random (1+ (ash (* 3 XDIM) -1))) (/ XDIM 2)))
	 (yloc (- (random (1+ (ash (* 3 YDIM) -1))) (/ YDIM 2)))
	 (xend (+ xloc (1+ (random XDIM))))
	 (yend (+ yloc (1+ (random YDIM)))))
    (when (and (>= xend 0) (>= yend 0))
       (setf xloc (max 0 xloc))
       (setf yloc (max 0 yloc))
       (setf xend (min XBORDER (max 0 xend)))
       (setf yend (min YBORDER (max 0 yend)))
       (dotimes (y (1+ (- yend yloc)))
	  (dotimes (x (1+ (- xend xloc)))
	     (add-autotroph (aref *env* (+ xloc x) (+ yloc y))
			    pheno amount))))))

(defun add-autotroph (loc pheno amount)
  (let ((orgs (location-orgs loc))
	(base (location-num-food loc)))
    (if (>= (+ amount base) MAX-FOOD-CAPACITY)
	(setf amount (- MAX-FOOD-CAPACITY 1 base))
      (setf *done-planting* ()))
    (dotimes (i amount)
       (setf (organism-phenotype (svref orgs (+ i base))) pheno)
       (setf (organism-alive (svref orgs (+ i base))) t))
    (incf (location-num-food loc) amount)))

(defun mix-seeds ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(mix-array (location-orgs (aref *env* x y)) 0
		   MAX-FOOD-CAPACITY))))

;;; I have to decide whether or not to construct the initial
;;; population so as to be viable, or just make them random.  It would
;;; be nice to see a food web form rather than impose it.
;;; MadDog says make it random - least arbitrary.

(defmacro populate-location (loc)
  `(let* ((orgs (location-orgs ,loc)))
     (dotimes (i MAX-ORG-CAPACITY)
       (let ((org (svref orgs (+ i MAX-FOOD-CAPACITY))))
	 (if (= (rem i NOAH-FACTOR) 0)
	     (init-org org)
	   (progn ; mates for the previous org.
	     (copy-org (+ i MAX-FOOD-CAPACITY)  ; a macro
		       (+ MAX-FOOD-CAPACITY
			  (* NOAH-FACTOR         
			     (floor i NOAH-FACTOR)))
		       orgs)
	     (setf (organism-last-meal org)
		   (- *clock* (random MAX-FAST)))
	     (setf (organism-num-meals org) (random REPRO-THRESHOLD))))))
     (setf (location-num-orgs ,loc) MAX-ORG-CAPACITY)))


(defun setup-population ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(populate-location (aref *env* x y)))))



(defun init-org (org)
  (setf (organism-alive org) t)
  (setf (organism-prey org) (generate-prey-gene)) ;macros.
  (setf (organism-generalism org) (generate-generalism-gene))
  (setf (organism-phenotype org) (generate-phenotype-gene))
  (setf (organism-last-meal org) (- *clock* (random MAX-FAST)))
  (setf (organism-num-meals org) (random REPRO-THRESHOLD)))
  

 
;;; Initializes the *probabilities* array with a poisson distribution
;;; for the number of bits expected to flip in a single gene.
(defun init-poisson ()
  (let* ((lambda (* GENESIZE *mutation-rate*))
	 (numer (exp (- lambda)))
	 (denom 1.0))
    (dotimes (i GENESIZE)
       (setf (svref *probabilities* i)
	     (if (= i 0) (/ numer denom)
	       (+ (svref *probabilities* (1- i))
		  (/ numer denom))))
       (setf numer (* numer lambda))
       (setf denom (* denom (1+ i))))))


;;; Determines the number of bits to flip in a gene.
(defun gene-poisson ()
  (let ((rand (uniform)))
    (dotimes (i (* NUMGENES GENESIZE) i)
       (when (< rand (svref *probabilities* i))
	     (return i)))))



;-----------------------#5,#6,#7Update Environment----------------------

;;; This must accomplish several things, depending on the parameters:
;;; 1. Replenish autotrophs (perhaps with seasonality) #7
;;; 2. Change barriers if physical env. is fragmenting #5, #6

(defun update-env ()
  ;(replenish-autotrophs)
  (when (or *endemnicity* *fragmentation* *fluctuating-env*)
	(alter-topography)))

(defun replenish-autotrophs ()
  (dotimes (y YDIM)
     (dotimes (x XDIM)
	(let* ((loc (aref *env* x y))
	       (num-food (modify-food-quantity (location-num-food loc)))
	       (orgs (location-orgs loc)))
	  (dotimes (i num-food)
	     (setf (organism-alive (svref orgs i)) t)))))))
  

(defun alter-topography ())

;-----------------------Update Organisms-----------------------

(defun update-orgs ()
  (let ((product-divisor (round (* *productivity* GENESIZE))))
    (dotimes (y YDIM)
       (dotimes (x XDIM)
	  (let* ((loc (aref *env* x y))
		 (orgs (location-orgs loc)))
	    (dotimes (i MAX-FOOD-CAPACITY)
		     (update-an-autotroph loc (svref orgs i) product-divisor))
	    (do ((i MAX-FOOD-CAPACITY (1+ i)))
		((>= i CAPACITY))
		(update-an-org loc (svref orgs i))))))))
  

;;; (round (organism-photosynth org)
;;;           divisor)
(defun update-an-autotroph (loc plant-org divisor))



  
(defun update-an-org (loc org))

;----------------------------Organisms--------------------------

(defun print-org (org)
  (unless (organism-alive org)
	  (format t "   *DEAD*  :( ~%"))
  (format t "   Phenotype: ~29,'0B~%" (organism-phenotype org))
  (format t "   Preferred Prey: ~29,'0B~%" (organism-prey org))
  (format t "   Generalism:     ~29,'0B~%" (organism-generalism org))
  (format t "   Number of meals eaten: ~a~%" (organism-num-meals org))
  (format t "   Time of last meal: ~a~%" (organism-last-meal org)))





;;; Edible? checks to see if the predator can eat the prey, and
;;; returns a boolean.

(defun edible? (pred prey)
  (when (and (organism-alive pred) (organism-alive prey))
	(let ((generalism (organism-generalism pred)))
	  (= (logior generalism (organism-prey pred))
	     (logior generalism (organism-phenotype prey))))))


;;; Eat assumes that the prey is a viable food source (edible? = true)
;;; for the predator.  It then updates the environment and organisms.
(defun eat (pred pred-loc prey-index prey-loc)
  (kill prey-index prey-loc)
  (setf (organism-last-meal pred) *clock*)
  (if (= (organism-num-meals pred) REPRO-THRESHOLD)
      (progn (setf (organism-num-meals pred) 0)
	     (reproduce pred pred-loc *dispersal*))
    (incf (organism-num-meals pred))))


(defun kill (org org-loc)
  (setf (organism-alive org) nil)
  (decf (location-num-orgs org-loc)))



;;; Reproduce first looks for a mate in the neighborhood.  If one is
;;; found, it then takes a random walk of length "disperse."  Once it has
;;; settled on a location, it looks for an empty slot in the organism
;;; array.  If that is found, it forms the new organism from the two
;;; parents and places it in that slot.
(defun reproduce (org org-loc disperse &optional (plant ()))
  (let ((mate (find-mate org (location-neighborhood org-loc) plant)))
    (unless (null mate)
       (let ((new-loc (random-walk (location-x org-loc)
				   (location-y org-loc)
				   (random disperse))))
	 (unless (>= (location-num-orgs new-loc) MAX-ORG-CAPACITY)
	    (let* ((fetus-generalism (crossover (organism-generalism org)
						(organism-generalism mate)))
		   (new-org (svref new-loc
				   (find-empty-org (location-organisms new-loc)
						   (start-placement
						    (generalism-rank
						     fetus-generalism)
						    plant)
						   plant))))
	      (unless (null new-org)
		      (place-org new-loc (make-fetus org mate
						     fetus-generalism)
				 new-org plant))))))))





(defmacro find-mate-in-loc (org loc plant)
  `(let ((mates (location-orgs ,loc)))
    (do ((i (start-placement (max 0 (- (generalism-rank
					(organism-generalism ,org))
				       *species-radius*)
			     ,plant)
	    (1+ i))))
	((if ,plant (>= i MAX-FOOD-CAPACITY) (>= i CAPACITY)) nil)
	(when (feasible-mate? ,org (svref mates i))
	      (return (svref mates i))))))


;;; This looks in the neighborhood for an organism that is a viable
;;; mate for the organism.  However, selfing is not allowed. 
(defun find-mate (org neighborhood plant)
  (dolist (loc neighborhood)
     (let ((potential (find-mate-in-loc org loc plant)))
       (when potential (return potential)))))





;;;***
(defun make-fetus (parent mate generalism)
  (setf (organism-generalism *fetus*) generalism)
  (setf (organism-prey *fetus*)
	(crossover (organism-prey parent)
		   (organism-prey mate)))
  (setf (organism-phenotype *fetus*)
	(crossover (organism-phenotype parent)
		   (organism-phenotype mate)))
  *fetus*)


;;; This carries out a crossover that is relatively unbiased.  2/3 of
;;; the time, there is only one crossover point in the gene, but 1/3
;;; of the time there are two crossover points.  This means that the
;;; further a bit is from another bit, the more likely they will be
;;; separated, but the endpoints are not necessarily separated as in
;;; one point crossover.
(defun crossover (gene1 gene2)
  (if (<= *species-radius* 1)
      (if (random 2) gene1 gene2)
    (let* ((pt1 (random (1+ GENESIZE)))
	   (pt2 (max GENESIZE (min 0 (- (random (* 3 GENESIZE)) GENESIZE)))))
      (if (< pt1 pt2)
	  (2p-crossover gene1 gene2 pt1 pt2)
	(2p-crossover gene1 gene2 pt2 pt1)))))
  

;;; This gives unbiased crossover.
(defun 2p-crossover (gene1 gene2 pt1 pt2)
  (let ((mask1 (logand GENEMASK (ash GENEMASK pt1))) ; 1111...11000...
	(mask2 (ash GENEMASK (- pt2 GENESIZE)))); 000111111111...
;    (format t "Crossing over at points ~a and ~a~%" pt1 pt2)
;    (print-bin mask1)
;    (print-bin mask2)
    (logand GENEMASK ; cut off high order bits.
	    (logior  ; paste together two halves.
	     (logior (logand gene1 (logandc2 GENEMASK mask1)) ;gene1's
		     (logand gene1 (logandc2 GENEMASK mask2)));contribution.
	     (logand mask1 (logand mask2 gene2)))))) ; the middle of gene2.



(defun start-placement (generalist-level &optional (plant ()))
  (if plant
      (if *specialist-advantage*
	  (* PLANT-RANK-CAPACITY generalist-level)
	0)
    (if *specialist-advantage*
	(+ MAX-FOOD-CAPACITY
	   (* RANK-CAPACITY generalist-level))
      MAX-FOOD-CAPACITY)))


;;; Place-org places initializes a dead organism structure with the
;;; values in the structure "org" and updates the environment.  This
;;; should only be called once the "barrier" has been checked and
;;; passed. 

(defun place-org (loc org new-org &optional (plant ()))
  (if plant
      (progn
	(incf (location-num-food loc))
	(when plant
	   (setf (organism-photosynth new-org)
		 (* (match (logior (organism-prey org)
				   (organism-generalism org))
			   (logior (organism-generalism org)
				   (location-climate loc)))
		    ENERGY-CONVERSION))))
    (incf (location-num-orgs loc)))
  (setf (organism-alive new-org) t)
  (setf (organism-last-meal new-org) *clock*)
  (setf (organism-num-meals new-org) 0)
  (setf (organism-phenotype new-org) (organism-phenotype org))
  (setf (organism-prey new-org) (organism-prey org))
  (setf (organism-generalism new-org) (organism-generalism org)))


;(defun find-empty-org (orgs index)
;  (do ((i index (1+ i)))
;      ((>= i CAPACITY) nil)
;      (when (not (organism-alive (svref orgs i)))
;            (return i))))

(defun find-empty-org (orgs index &optional (plant ()))
  (cond ((and plant (>= index MAX-FOOD-CAPACITY)) nil)
	((>= index CAPACITY) nil)
        ((not (organism-alive (svref orgs index))) index)
        (t (find-empty-org orgs (1+ index)))))

	
;;; This computes the effective genotype of an organism by masking the
;;; prey gene with the generalism gene.  This is used to calculate
;;; species identity.
(defun genotype (org)
  (let ((pheno (ash (organism-phenotype org) GENESIZE))
	(prey (logior (organism-prey org)
		      (organism-generalism org))))
    (+ pheno prey)))


;--------------------#8 Genetic Reduction in Gene Flow-----------------

(defun feasible-mate? (me org)
  (when (and (organism-alive me)
	     (organism-alive org)
	     (not (eq me org)))  ; selfing not allowed.
	(<= (+ (logcount (logxor (organism-phenotype me)
				 (organism-phenotype org)))
	       (logcount (logxor (organism-prey me)
				 (organism-prey org)))
	       (logcount (logxor (organism-generalism me)
				 (organism-generalism org))))
	    *species-radius*)))

; -------------------#8 Behavioral Reduction in Gene Flow--------------

(defun random-walk (x y steps)
  (format t "I'm at ~a ~a with this many steps left: ~a~%" x y steps)
  (if (<= steps 0) (aref *env* x y)
    (let ((new-x (min XBORDER (max 0 (+ x (1- (random 3))))))
	  (new-y (min YBORDER (max 0 (+ y (1- (random 3)))))))
      (if (< (random MAX-TERRAIN)
	     (location-barrier (aref *env* new-x new-y)))
	  (random-walk x y (1- steps))
	(random-walk new-x new-y (1- steps))))))



;------------------#9 and #10 Adaptive Space---------------------------

;;; I will want to do something fancy with this to set aside some
;;; range of bits to be "unoccupied" adaptive space.

;;; Mutates each of the three genes of the organism.
(defun mutate-org (org)
  (setf (organism-phenotype org)
	(mutate-gene (organism-phenotype org) 
		     (if *increasing-adaptive-space*
			 GENESIZE ADAPTZONE)))
  (setf (organism-prey org)
	(mutate-gene (organism-prey org)
		     (if *increasing-adaptive-space*
			 GENESIZE ADAPTZONE)))
  (setf (organism-generalism org)
	(mutate-generalism (organism-generalism org)
			   (if *increasing-adaptive-space*
			       GENESIZE ADAPTZONE))))

(defun mutate-gene (gene range)
  (let ((num-bits (gene-poisson)))
    (mutate-bits gene num-bits range))) ;a macro, flips this many bits.


;;; We are assuming that in the non-*niche-subdivision* case, all
;;; organisms will have 1's in the NICHEZONE of their generalism
;;; gene. 
(defun mutate-generalism (gene range)
  (let ((num-bits (gene-poisson)))
    (if *niche-subdivision*
	(dotimes (i num-bits gene)
            (let ((locus (random range)))
	      (setf gene (logxor (ash 1 locus) gene))))
      (dotimes (i num-bits gene)
	 (let ((locus (+ SUBDIVOFFSET (random NICHEZONE))))
	   (setf gene (logxor (ash 1 locus) gene)))))))