+# default number of generations to run. typically over ridden on the command
+# line
+GENS=1
+
+# how to connect to the database
+PSQL= psql "service=afpopgen" -e -1 -v ON_ERROR_STOP=1
+
+# Scenario to run. S to make the command line shorter
+# can be set via make override or environment
+#S=
+
+# scenario iteration, has to be 0 for updates, could be
+# non zero for image dumps and such
+N ?= 0
+
+# a target iteration number for finishing a run
+#T=
+
+export S
+
+# the default target
+list:
+ @echo Scenarios Available
+ @echo -------------------
+ @find scenarios -type f -name '*.sql' -printf '%P\n'
+
+images/$(M)/$(N)/$(M).$(N).mp4:
+ rm -f $@
+ ffmpeg -loglevel warning -framerate 24/1 -i images/$(M)/$(N)/ann%04d.png -c:v libx264 -pix_fmt yuv420p $@
+
+images:
+ utils/fetchimg $(S) $(N)
+
+annotate:
+ utils/annotate $(S) $(N)
+
+# prepare the database for running a model
+prepdb:
+
+# load a scenario into the database
+popmodel:
+ $(PSQL) -f scenarios/$(S).sql
+
+run: imclean popmodel nextgenimg images movies/$(S).$(N).mp4
+
+nextgen:
+ set -e ; for n in {1..$(GENS)} ; do $(PSQL) -v scen=$(S) -f model/generation.sql || exit 1 ; done
+
+nextgenimg:
+ for n in {1..$(GENS)} ; do $(PSQL) -v scen=$(S) -f model/generation.sql && $(PSQL) -v scen=$(S) -f model/dumpgenimg.sql ; done
+
+
+imclean:
+ rm -rf imgout/*
+ #$(PSQL) -c 'truncate ap.images'
+
+renumber:
+ # bail if N = T
+ # bail if any S N T blank
+ # just dump the image files, we can get them from the DB again if we need them
+ if test -d images/$(S)/$(N) ; then rm -r images/$(S)/$(N) ; fi
+ # These should be done together really
+ $(PSQL) -c "update ap.popgen set run = $(T) where run = $(N) and scenario = '$(S)'"
+ $(PSQL) -c "update ap.images set run = $(T) where run = $(N) and scenario = '$(S)'"
+ if test -f movies/$(S).$(N).mp4 ; then mv movies/$(S).$(N).mp4 movies/$(S).$(T).mp4 ; fi
+
+#more: nextgenimg images movies/$(S).$(N).mp4
+more: nextgen
+
+movies/$(S).$(N).mp4: images/$(S)/$(N)
+ mkdir -p movies
+ ffmpeg -y -loglevel warning -framerate 24/1 -i images/$(S)/$(N)/$(S)$(N)-%04d.png -c:v libx264 -pix_fmt yuv420p $@
+
+movie: genimages images movies/$(S).$(N).mp4
+
+update: images freq.mp4
+
+small0000.png: ann/ann0000.png
+ convert $< -background white -alpha remove -resize 200x200\> $@
+small1000.png: ann/ann1000.png
+ convert $< -background white -alpha remove -resize 200x200\> $@
+
+web/index.html: index.md
+ multimarkdown $< > $@
+
+# smallcap.mp4 largecap.mp4 layout.css smallcap0000.png smallcap1000.png largecap0000.png largecap1000.png
+
+website: web/index.html
+ rsync -v --progress $+ granicus.if.org:public_html/africa
+ #rsync -rv --progress --exclude '.*' --delete references/ granicus.if.org:public_html/africa/references
+
+freq.mp4:
+ rm -f $@
+ ffmpeg -loglevel warning -framerate 10/1 -i imgout/freq%04d.png -c:v libx264 -pix_fmt yuv420p $@
+
+ann.mp4:
+ rm -f $@
+ ffmpeg -loglevel warning -framerate 24/1 -i ann/ann%04d.png -c:v libx264 -pix_fmt yuv420p $@
+
+base:
+ $(PSQL) -f model/base.sql
+
+imgsetup:
+ $(PSQL) -f model/image.sql
+
+genimages:
+ $(PSQL) -f model/dumpgenimg.sql
+
+africadata:
+ $(PSQL) -f model/africa.sql
+
+colors:
+ $(PSQL) -f model/colors.sql
+
+reloadall: africadata imgsetup
+ $(PSQL) -f model/colors.sql
+ $(PSQL) -f model/base.sql
+
notes.pdf: notes.tex notes.bib
xelatex $<
bibtex notes
clean:
rm -f notes.pdf
+.PHONY: annotate images
--- /dev/null
+--select ap.xy2ijk(1,2);
+--select ap.ijk2xy(ap.xy2ijk(1,2));
+
+--select ap.ijk2xy(ap.adjhex(dir, ap.xy2ijk(1,2))) from generate_series(0,5) as dir;
+
+\timing on
+
+-- build an adjacency list for the hexes
+create temp table hexdata as
+with AL as (
+ select H.cantor, hex.adjacent(cantor) as adjacent
+ from ap.afrhex H
+)
+select
+AL.cantor, array_agg(AL.adjacent) as adjlist
+from AL, ap.afrhex AH
+where AL.adjacent = AH.cantor
+group by AL.cantor
+;
+
+
+\q
+select H.cantor, hex.adjacent(H.cantor) as adjacent
+from ap.afrhex H
+where H.cantor = 123;
+
+select * from ap.hexdata limit 10;
+select * from ap.hexdata where cantor = 123;
+
--- /dev/null
+--INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 102022, 'esri', 102022, '+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ', 'PROJCS["Africa_Albers_Equal_Area_Conic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",25],PARAMETER["Standard_Parallel_1",20],PARAMETER["Standard_Parallel_2",-23],PARAMETER["latitude_of_center",0],UNIT["Meter",1],AUTHORITY["EPSG","102022"]]');
+
+-- set hex width in km
+\set hexsize 100
+
+drop schema if exists ap cascade;
+create schema if not exists ap;
+
+set search_path to ap,gis,public;
+
+-- first we create a general climate information table
+create table climate (
+ climate text primary key,
+ basecapacity integer, -- in people per 100 km
+ capacity integer, -- adjusted for model
+ color integer[] -- rgb
+);
+
+-- this is in persons per 100 sq km
+insert into climate (climate, basecapacity, color) values
+('Af', 3, ARRAY[161,26,20]), -- tropical rain forest
+('Am', 45, ARRAY[255,33,26]), -- tropical monsoon
+('As', 50, ARRAY[255,163,150]), -- Savanna dry, 400
+('Aw', 125, ARRAY[255,205, 199]), -- Savanna wet, up to 1000
+('BSh', 12, ARRAY[211,139,27]), -- Steppe, hotter
+('BSk', 12, ARRAY[207,169,81]), -- Steppe, colder
+('BWk', 3, ARRAY[254,245,91]), -- Desert, colder, 0 - 120
+('BWh', 3, ARRAY[255,197,13]), -- Desert, hotter, 0 - 120
+('Cfa', 125, ARRAY[14,43,13]), -- Warm temperate
+('Cfb', 50, ARRAY[21,68,22]), -- Maritime, 300-600
+('Csa', 50, ARRAY[17,214,16]), -- Med
+('Csb', 50, ARRAY[140,230,18]), -- Med
+('Cwa', 125, ARRAY[184,103,24]), -- Warm temperate
+('Cwb', 125, ARRAY[154,100,25]) -- Maritime
+;
+
+-- project the africa polygon into the SRS we will use for the map
+-- africa is polygon 3 in the gshhs table
+create temp table afrp as select gid, st_transform(geom,102022) from gshhs.gshhs_intermediate where gid = 3;
+
+create table afrhex (
+ x integer,
+ y integer,
+ hex numeric primary key,
+ id text,
+ climate text references climate,
+ geom geometry,
+ adjacent numeric[]
+);
+
+insert into afrhex (x,y,hex,geom)
+select * from hex.hexgrid(0,0,0,(select geom from gshhs.gshhs_intermediate where gid = 3), 102022, NULL, :hexsize * 1000);
+
+create index afrhex_index on afrhex using gist(geom);
+create index afrhex_cantor on afrhex (hex);
+analyze afrhex;
+analyze koeppengeiger;
+
+-- hexgrid creates a rectangular grid
+-- so we remove hexes that don't intersect with the africa polygon
+delete from afrhex where not st_intersects(geom, (select st_transform from afrp));
+
+-- create a helper table for the climate calculation
+create temp table cli as
+select H.hex, C.climate, st_area(ST_intersection(H.geom, st_transform(C.geom,102022))) as area
+from
+afrhex H, koeppengeiger C
+where
+C.geom && (select geom from gshhs.gshhs_intermediate where gid = 3)
+and H.geom && st_transform(C.geom, 102022)
+and ST_intersects(H.geom, st_transform(C.geom,102022))
+;
+
+create temp table maxcli as
+select hex as hexc, climate, area, row_number() over (partition by hex order by area desc) from cli order by hex, area desc
+;
+
+-- create a climate table
+-- update afrhex with climate data
+update afrhex set climate = (select climate from maxcli where hex = hexc and row_number = 1);
+
+-- build an adjacency list for the hexes
+create temp table hexdata as
+with AL as (
+ select H.hex, hex.adjacent(hex) as adjacent
+ from afrhex H
+)
+select
+AL.hex, array_agg(AL.adjacent) as adjlist
+from AL, afrhex AH
+where AL.adjacent = AH.hex
+group by AL.hex
+;
+
+update afrhex H set adjacent = (select L.adjlist from hexdata L where L.hex = H.hex);
+
+\q
+-- create elevation data
+create table ap.afrelev as
+with mean as
+(select cantor as hexc,
+(st_summarystats(ST_clip(st_transform(G.rast,102022),H.hex))).*
+from
+ap.afrhex H inner join gtopo30 G
+on
+st_intersects(H.hex, st_transform(G.rast,102022))
+where H.hex && st_transform(G.rast,102022)
+and G.rast && (select geom from gshhs.gshhs_intermediate where gid = 3)
+)
+select hexc, min(min), max(max), sum(mean*count)/sum(count) as avg
+from mean
+where count > 0
+group by hexc
+;
+
+-- update afrhex with elevation data
+alter table ap.afrhex add elevation integer;
+-- todo parallel this
+update ap.afrhex set elevation = (select avg from ap.afrelev where cantor = hexc)::integer;
+--update ap.afrhex set elevation = (select
+(st_summarystats(ST_clip(st_transform(G.rast,102022),H.hex))).avg
+ avg from ap.afrelev where cantor = hexc)::integer;
--- /dev/null
+set search_path to ap,r,gis;
+
+-- adjust capacity for hex size and double because our populations
+-- are haploid
+update climate set capacity = basecapacity * (select st_area(geom)/1000000/100 from afrhex limit 1) * 2;
+
+drop table if exists popgen cascade;
+create table popgen (
+ scenario text not null default 'current',
+ run integer not null default 0,
+ gen integer not null, -- the generation number
+ hex numeric not null references ap.afrhex,
+ freq integer[] not null, -- the popcounts of each allell
+ climate text not null references climate -- the climate can vary
+);
+
+-- create initial population
+-- TODO probably not necessary
+insert into popgen (scenario, run, gen, hex, freq, climate)
+select 'base', 0, 0, H.hex, ARRAY[0,0], H.climate
+from afrhex H
+;
+
+create index popgen_model_index on popgen(scenario);
+create index popgen_gen_index on popgen(gen);
+create index popgen_hex_index on popgen(hex);
+create unique index popgen_sgrh_index on popgen(scenario,gen,run,hex);
+create index popgen_pop_index on popgen (r.vec_sum(freq));
+
+drop view if exists pop;
+
+-- helper view for image construction
+create view pop as
+select P.scenario, P.gen, C.capacity, P.climate, P.freq, A.geom,
+ P.freq[1] as popa, -- first allell freq
+ vec_sum(P.freq) - P.freq[1] as popb, -- remaining freq
+case when P.freq[1] = 0 then 0.0
+ else P.freq[1]::float / vec_sum(P.freq)
+end as popafrac, -- fraction of total for first allell
+vec_sum(P.freq) as totalpop,
+-- fraction of maximum population in a hex any generation
+coalesce(
+ vec_sum(P.freq)::float / max(vec_sum(P.freq)) over (rows between unbounded preceding and unbounded following),
+ 0.0) as maxfrac,
+-- fraction of maximum population in a hex this generation
+coalesce(
+ vec_sum(P.freq)::float / max(vec_sum(P.freq)) over (partition by gen rows between unbounded preceding and unbounded following), 0.0) as maxgenfrac,
+
+least(1.0, (vec_sum(P.freq))::float/C.capacity) as pressure -- % of maximum capacity, capped at 100%
+from
+popgen P join climate C on C.climate = P.climate join afrhex A on A.hex = P.hex
+;
--- /dev/null
+set search_path to ap;
+
+create or replace function rgb(hsv float[]) returns float[] language 'plpgsql' as $$
+declare
+ i integer;
+ f float;
+ p float;
+ q float;
+ t float;
+ h float;
+ v float;
+ rgb float[];
+begin
+ v := hsv[3];
+
+ if hsv[2] = 0.0 then
+ rgb[1] = hsv[3];
+ rgb[2] = hsv[3];
+ rgb[3] = hsv[3];
+ return rgb;
+ end if;
+
+ h := hsv[1] * 6.0;
+ i := floor(h)::integer;
+ f := h - i;
+ p := hsv[3] * ( 1.0 - hsv[2] );
+ q := hsv[3] * ( 1.0 - hsv[2] * f );
+ t := hsv[3] * ( 1.0 - hsv[2] * (1.0-f) );
+
+ if i = 0 then rgb = ARRAY[v,t,p];
+ elsif i = 1 then rgb = ARRAY[q,v,p];
+ elsif i = 2 then rgb = ARRAY[p,v,t];
+ elsif i = 3 then rgb = ARRAY[p,q,v];
+ elsif i = 4 then rgb = ARRAY[t,p,v];
+ else rgb = ARRAY[v,p,q];
+ end if;
+
+ return rgb;
+end;
+$$ immutable strict;
+
+create or replace function rgbb(rgb float[]) returns integer[] language 'sql' as $$
+ select ARRAY[(rgb[1] * 255)::integer,(rgb[2] * 255)::integer,(rgb[3] * 255)::integer];
+$$ immutable strict;
--- /dev/null
+\set scale 20000.0
+
+set search_path to ap,gis;
+
+\timing on
+
+with info as (
+ SELECT P.scenario, P.run, P.gen,
+ c.capacity,
+ P.climate,
+ p.freq,
+ a.geom,
+ P.freq[1] AS popa,
+ r.vec_sum(p.freq) - p.freq[1] AS popb,
+ CASE
+ WHEN p.freq[1] = 0 THEN 0.0::double precision
+ ELSE p.freq[1]::double precision / r.vec_sum(p.freq)::double precision
+ END AS popafrac,
+ r.vec_sum(p.freq) AS totalpop,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxfrac,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (PARTITION BY p.gen ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxgenfrac,
+ LEAST(1.0::double precision, r.vec_sum(p.freq)::double precision / c.capacity::double precision) AS pressure
+ FROM popgen p
+ JOIN climate c ON c.climate = p.climate
+ JOIN afrhex a ON a.hex = p.hex
+ left join images I on I.scenario = P.scenario and I.run = P.run and I.gen = P.gen
+ where I.run is null and I.run is null and I.gen is null
+ --where P.gen not exists (select 1 gen from images where scenario = :'scen' and gen = P.gen and run = )
+)
+insert into images (scenario, run, gen, png)
+select
+P.scenario, P.run, P.gen,
+st_aspng(
+-- merge(
+ st_union(
+ st_asraster(
+ P.geom, 20000.0, -20000.0, 0,0,ARRAY['8BUI','8BUI','8BUI','8BUI'],
+ -- next is constructing the color
+ -- yellow to blue or so is hue .16 - .66, sat .6 val 1.0
+ rgbb(rgb(ARRAY[
+ -- (1.6 / (1 + exp(-( P.popafrac - 0.5 )))-0.5) / 2.2 + 0.15, -- sigmoid bias to ends of hue
+ (4 / (1 + exp(-(P.popafrac - 0.5 ))) - 1.5) / 2.0 + 0.075,
+ -- P.popafrac/2.0+0.16, -- hue
+ --0.9*(P.pressure), -- saturation
+ 1.0, -- saturation
+ --1.0*(P.maxgenfrac) -- value
+ case when P.totalpop > 0 then
+ -- clamp max to 1.0, just in case the maxgen frac calculation was wrong
+ least(1.0, 0.15+0.85*sqrt(P.maxgenfrac)) -- value
+ --sqrt(P.pressure)
+ else
+ 0.0
+ end
+ --1.0 -- value
+ ]
+ )) || ARRAY[0],
+ --rgbb(rgb(ARRAY[0.0/2.0+0.16, 0.6*0.0, 1.0*0.0])) || ARRAY[0],
+ --rgbb(ARRAY[0,1,0])||ARRAY[0],
+ ARRAY[255,255,255,0], 0, 0, true
+ )
+ )
+-- )
+ ,ARRAY[1,2,3]
+) as png
+from info P
+--and P.totalpop > 0
+group by P.scenario, P.run, P.gen
+;
--- /dev/null
+set search_path to ap;
+
+-- create elevation data
+\timing
+--explain
+create table elev as
+select st_transform(rast,102022) as rast
+from gtopo30
+where
+rast && (select geom from gshhs.gshhs_intermediate where gid = 3)
+;
+CREATE INDEX elev_index ON elev USING gist( ST_ConvexHull(rast) );
+
+-- can we parallelize this one?
+-- takes about 3.5 min
+create table stats as
+select cantor as hexc,
+(st_summarystats(ST_clip(G.rast,H.hex))).*
+from
+afrhex H inner join elev G
+on
+st_intersects(H.hex, G.rast)
+where H.hex && G.rast
+;
+
+create index stats_hex_index on stats(hexc);
+analyze stats;
+alter table afrhex add elevation integer;
+update afrhex set elevation = (select sum(mean*count)/sum(count) from stats where cantor = hexc)::integer;
+
+\q
--- /dev/null
+-- increase or decrease pop in each hex
+-- distribute surplus population to one neighboring cell
+
+set search_path to ap,r;
+
+-- set up base growth
+drop table if exists nextgrowth;
+create table nextgrowth (
+ hex numeric not null,
+ gen integer not null,
+ freq integer[] not null,
+ climate text references climate,
+ capacity integer,
+ attract integer
+);
+create unique index nextgrowth_id_index on nextgrowth(hex,gen);
+
+create or replace function hexcap(hexid numeric) returns integer language sql as $$
+select C.capacity
+from afrhex H left join climate C on C.climate = H.climate
+where H.hex = hexid
+$$;
+
+\timing on
+insert into nextgrowth (hex,gen,freq,climate,capacity)
+ select P.hex, gen+1,
+ rmultinom(ceiling(
+ -- TODO calculate growth here
+ case when hexcap(P.hex) < vec_sum(P.freq) then 0.6 else
+ 1.15
+ end
+ * (floor(vec_sum(P.freq))/2.0))::integer*2,P.freq),
+ P.climate,
+ C.capacity
+ from popgen P
+ left join climate C on C.climate = P.climate
+ where gen = (select max(gen) from popgen where scenario = :'scen' and run = 0)
+ and scenario = :'scen'
+ and run = 0
+ and vec_sum(P.freq) > 0
+ union all
+ select P.hex, gen+1, ARRAY[0,0], P.climate, C.capacity
+ from popgen P
+ left join climate C on C.climate = P.climate
+ where gen = (select max(gen) from popgen where scenario = :'scen' and run = 0)
+ and scenario = :'scen'
+ and run = 0
+ and vec_sum(P.freq) <= 0
+;
+-- TODO min attractiveness sqrt(capacity)?
+update nextgrowth set attract = greatest(sqrt(capacity), (capacity - vec_sum(freq)) * (1.0 - vec_sum(freq) / capacity));
+create index nextgrowth_hex_index on nextgrowth (hex);
+create index nextgrowth_gen_index on nextgrowth (gen);
+create index nextgrowth_pop_index on nextgrowth (vec_sum(freq));
+analyze nextgrowth;
+
+drop table if exists migration;
+create table migration (
+ source numeric,
+ dest numeric,
+ freq integer[]
+);
+
+create or replace function attractiveness(thex numeric) returns integer language 'sql' as $$
+select N.attract from nextgrowth N where N.hex = thex;
+$$ stable;
+
+create or replace function attractiveness(hexes numeric []) returns integer[] language 'sql' as $$
+select array_agg(attractiveness(h)) from unnest(hexes) h;
+$$ stable;
+
+with migrants as (
+ select
+ N.hex as source, C.capacity, N.freq,
+ ((vec_sum(N.freq)-C.capacity)/2)::integer as ind -- people going out
+ from nextgrowth N
+ left join afrhex H on N.hex = H.hex
+ left join climate C on C.climate = H.climate
+ where C.capacity < vec_sum(N.freq)
+)
+insert into migration
+select M.source,
+ -- pick one adjacent hex, TODO pick more than one once we can do multivariate hypergeom
+ (sample(A.adjacent, 1, attractiveness(A.adjacent)))[1]::numeric as dest,
+ -- so distribute 2 allells per migrating individual
+ -- use hypergeo since we can't have more migrate than exist
+ rmhypergeo(M.ind * 2,M.freq)
+ from migrants M
+ left join afrhex A on A.hex = M.source
+;
+
+with migrants as (
+ select M.dest as hex,
+ vec_add(M.freq) as freq
+ from migration M group by M.dest
+union all
+ select M.source,
+ -- make these negative
+ vec_sub(0,vec_add(M.freq)) as freq
+ from migration M group by M.source
+), net as (
+ select M.hex, vec_add(M.freq) as freq
+ from migrants M group by M.hex
+)
+insert into popgen (scenario,hex,gen,climate,freq)
+-- two parts: has migration and no migration
+ select
+ :'scen',
+ N.hex, N.gen, N.climate,
+ vec_add(N.freq, M.freq)
+ from nextgrowth N left join net M on M.hex = N.hex
+ where M.hex is not null
+union all
+ select
+ :'scen',
+ N.hex, N.gen, N.climate,
+ N.freq
+ from nextgrowth N left join net M on M.hex = N.hex
+ where M.hex is null
+;
+
+analyze popgen;
+-- truncate to capacity, may not be needed if the "growth" function also shrinks
--- /dev/null
+set role gis;
+
+drop schema if exists hex cascade;
+create schema hex;
+
+set search_path to hex,public,gis;
+
+create or replace function hexgrid(
+ xdim integer,
+ ydim integer,
+ angle numeric,
+ coverage geometry,
+ srid integer,
+ origin geometry(POINT),
+ tilesize numeric
+) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
+declare
+ rcover geometry;
+ stride double precision;
+ x integer;
+ y integer;
+ xmin double precision;
+ ymin double precision;
+ xoff double precision;
+begin
+ -- the origin will be the center of the 1,1 hex.
+ -- if null, then the origin will be set to fit in the
+ -- covering box.
+
+ -- xdim, and ydim can be null, in which case they are
+ -- calculated from the tilesize and the covering box
+
+ -- calculate ydim
+
+ -- precedence: ydim, xdim, tilesize
+
+ -- counter rotate geometry, fill, rotate grid
+ rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
+
+ xmin := ST_XMin(rcover);
+ ymin := ST_YMin(rcover);
+
+ stride := tilesize * 2.0 / sqrt(3.0);
+
+ --raise notice 'stride = %', stride;
+ xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
+ ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/tilesize)::integer;
+
+ --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
+
+ -- TODO reorient these so right hand rule is inside
+ hex := ST_Scale(ST_GeomFromText('POLYGON((
+.577350269189625764509148780502 0.0,
+.288675134594812882254574390251 0.5,
+ -.288675134594812882254574390251 0.5,
+ -.577350269189625764509148780502 0.0,
+ -.288675134594812882254574390251 -0.5,
+ .288675134594812882254574390251 -0.5,
+ .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);
+
+ xoff := .577350269189625764509148780502 * 1.5;
+
+ for x in 0 .. xdim loop
+ for y in 0 .. ydim loop
+ -- TODO actually calculate cantor
+ return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]),
+ ST_Translate(hex,
+ tilesize * (x * xoff) + xmin,
+ --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
+ tilesize * (y + 0.5 * (x%2) ) + ymin
+ )
+ --st_rotate(ST_Translate(hex,
+ --tilesize * (x * xoff) + xmin,
+ --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
+ --), -3.1415926/6)
+ ;
+ end loop;
+ end loop;
+END;
+$$
+language 'plpgsql'
+;
+
+-- create some hexagon ops
+
+create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ ijk integer[];
+begin
+ if array_length(xy,1) = 3 then
+ return xy;
+ end if;
+
+ ijk[1] := xy[1];
+ ijk[2] := -xy[2];
+
+ if xy[1] < 0 then
+ ijk[2] := ijk[2] + (-xy[1]+1) / 2;
+ else
+ ijk[2] := ijk[2] - xy[1]/2;
+ end if;
+
+ ijk[3] := -ijk[1] - ijk[2];
+
+ return ijk;
+end
+$$ immutable strict;
+
+create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
+ select hex.ijk(ARRAY[x,y]);
+$$ immutable strict;
+
+-- invert the cantor paring
+create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
+declare
+ xy integer[];
+ w integer;
+ t integer;
+ py integer;
+begin
+ if hex >= 0 then
+ hex := hex-1;
+ w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
+ t := (w * w + w)/2;
+ py := hex - t;
+
+ xy = ARRAY[w-py,py];
+ -- handle the negative cases
+ else
+ hex := -hex;
+ hex := hex-1;
+ w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
+ t := (w * w + w)/2;
+ py := hex - t;
+
+ xy = ARRAY[w-py,py];
+
+ if xy[1] % 2 = 1 then
+ xy[1] = -(xy[1]+1)/2;
+ else
+ xy[1] = xy[1]/2;
+ end if;
+ if xy[2] % 2 = 1 then
+ xy[2] = -(xy[2]+1)/2;
+ else
+ xy[2] = xy[2]/2;
+ end if;
+ end if;
+
+ return xy;
+end;
+$$ immutable strict;
+
+create function ijk(hex numeric) returns integer[] language 'sql' as $$
+ select hex.ijk(hex.xy(hex));
+$$ immutable strict;
+
+create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ xy integer[];
+begin
+ -- make sure this is really an ijk
+ -- just no-op if it's already an
+ -- xy, so this is sort of a coercion function
+ if array_length(ijk,1) = 2 then
+ return ijk;
+ end if;
+ xy[1] := ijk[1];
+ xy[2] := -ijk[2];
+
+ if ijk[1] < 0 then
+ xy[2] := xy[2] + (-ijk[1] + 1) / 2;
+ else
+ xy[2] := xy[2] - ijk[1]/2;
+ end if;
+ return xy;
+end;
+$$;
+
+CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
+RETURNS integer[] AS
+$$
+SELECT array_agg(result)
+FROM (SELECT tuple.val1 + tuple.val2 AS result
+ FROM (SELECT UNNEST($1) AS val1
+ ,UNNEST($2) AS val2
+ ,generate_subscripts($1, 1) AS ix) tuple
+ ORDER BY ix) inn;
+$$ LANGUAGE SQL IMMUTABLE STRICT;
+
+create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ adj integer[];
+begin
+ if array_length(hex,1) = 2 then
+ hex = hex.ijk(hex);
+ end if;
+
+ dir := dir % 6;
+ if dir = 2 then
+ hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
+ elsif dir = 1 then
+ hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
+ elsif dir = 0 then
+ hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
+ elsif dir = 5 then
+ hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
+ elsif dir = 4 then
+ hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
+ elsif dir = 3 then
+ hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
+ else
+ -- error
+ end if;
+ return hex;
+end;
+$$ immutable strict;
+
+create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
+ select hex.adjhex(dir, hex.ijk(hex));
+$$ immutable strict;
+
+-- TODO make this work for any, or inline it
+create function natcantor(hex integer[]) returns numeric language 'sql' as $$
+ select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+$$ immutable strict;
+
+create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
+begin
+ -- convert to xy
+ if array_length(hex,1) = 3 then
+ hex = hex.xy(hex);
+ end if;
+
+ if hex[1] < 0 or hex[2] < 0 then
+ hex[1] = hex[1] * 2;
+ hex[2] = hex[2] * 2;
+ if hex[1] < 0 then
+ hex[1] = -hex[1] - 1;
+ end if;
+ if hex[2] < 0 then
+ hex[2] = -hex[2] - 1;
+ end if;
+ return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+ end if;
+
+ return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+end;
+$$ immutable strict;
+
+create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
+select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
+$$ immutable strict;
+
+create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
+ select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
+$$ immutable strict;
+
+create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
+declare
+ c integer[]; -- cantor coordinates
+ base integer[];
+begin
+ if range <= 0 then
+ return;
+ end if;
+ if include then
+ return next hex;
+ end if;
+
+ -- we could in theory eliminate the inner loop and
+ -- do a return query
+ base = hex.ijk(hex);
+ for r in 1 .. range loop
+ c = base;
+ c[1] := c[1] + r;
+ c[3] := -c[1] - c[2]; -- keep the invariant
+ hex := hex.cantor(c);
+
+ for h in 0 .. r * 6 - 1 loop
+ return next hex;
+ hex = hex.adjacent(h/r+2,hex);
+ end loop;
+ end loop;
+ return;
+end;
+$$ immutable strict;
+
+create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
+select array_agg(c) from hex.within(hex, range, include) c;
+$$ immutable strict;
+
+create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
+declare
+ c integer[]; -- cantor coordinates
+begin
+ if range = 0 then
+ return next hex;
+ end if;
+ if range <= 0 then
+ return;
+ end if;
+
+ c = ijk(hex);
+ c[1] := c[1] + range;
+ c[3] := -c[1] - c[2]; -- keep the invariant
+ hex := hex.cantor(c);
+
+ for q in 0 .. range * 6 - 1 loop
+ return next hex;
+ hex = hex.adjacent(q/range+2,hex);
+ end loop;
+ return;
+end;
+$$ immutable strict;
+
+select
+xy(ARRAY[1,2]),
+ijk(xy(ARRAY[1,2])),
+xy(ijk(xy(ARRAY[1,2])))
+;
+
+select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
+select
+xy(adjhex(dir, ijk(ARRAY[1,2]))),
+cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
+xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
+from generate_series(0,5) as dir;
+
+select cantor(ARRAY[4,3]);
+--select * from within(cantor(ARRAY[4,4]), 2);
+select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;
--- /dev/null
+drop schema if exists ap cascade; -- africa paper
+create schema ap;
+
+-- create some hexagon ops
+create function ap.xy2ijk(x integer, y integer) returns integer[] language 'plpgsql' as $$
+declare
+ ijk integer[];
+begin
+ ijk[1] := x;
+ ijk[2] := -y;
+
+ if x < 0 then
+ ijk[2] := ijk[2] + (-x+1) / 2;
+ else
+ ijk[2] := ijk[2] - x/2;
+ end if;
+
+ ijk[3] := -ijk[1] - ijk[2];
+
+ return ijk;
+end
+$$;
+
+create function ap.xy2ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ ijk integer[];
+begin
+ ijk[1] := xy[1];
+ ijk[2] := -xy[2];
+
+ if xy[1] < 0 then
+ ijk[2] := ijk[2] + (-xy[1]+1) / 2;
+ else
+ ijk[2] := ijk[2] - xy[1]/2;
+ end if;
+
+ ijk[3] := -ijk[1] - ijk[2];
+
+ return ijk;
+end
+$$;
+
+create function ap.ijk2xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ xy integer[];
+begin
+ xy[1] := ijk[1];
+ xy[2] := -ijk[2];
+
+ if ijk[1] < 0 then
+ xy[2] := xy[2] + (-ijk[1] + 1) / 2;
+ else
+ xy[2] := xy[2] - ijk[1]/2;
+ end if;
+ return xy;
+end;
+$$;
+
+CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
+RETURNS integer[] AS
+$$
+SELECT array_agg(result)
+FROM (SELECT tuple.val1 + tuple.val2 AS result
+ FROM (SELECT UNNEST($1) AS val1
+ ,UNNEST($2) AS val2
+ ,generate_subscripts($1, 1) AS ix) tuple
+ ORDER BY ix) inn;
+$$ LANGUAGE SQL IMMUTABLE STRICT;
+
+create function ap.adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ adj integer[];
+begin
+ if dir = 2 then
+ hex := vec_add(hex, ARRAY[-1, 1, 0]);
+ elsif dir = 1 then
+ hex := vec_add(hex, ARRAY[ 0, 1,-1]);
+ elsif dir = 0 then
+ hex := vec_add(hex, ARRAY[ 1, 0,-1]);
+ elsif dir = 5 then
+ hex := vec_add(hex, ARRAY[ 1,-1, 0]);
+ elsif dir = 4 then
+ hex := vec_add(hex, ARRAY[ 0,-1, 1]);
+ elsif dir = 3 then
+ hex := vec_add(hex, ARRAY[-1, 0, 1]);
+ else
+ -- error
+ end if;
+ return hex;
+end;
+$$;
+
+create function ap.xy(hex numeric) returns integer[] language 'sql' as $$
+ select ARRAY[(hex/1000)::integer, (hex%1000)::integer];
+$$;
+create function ap.cantor(hex integer[]) returns numeric language 'sql' as $$
+ select (hex[1] * 1000 + hex[2])::numeric;
+$$;
+
+create function ap.adjacent(hex numeric) returns setof numeric language 'sql' as $$
+select ap.cantor(ap.ijk2xy(ap.adjhex(dir, ap.xy2ijk(ap.xy(hex))))) from generate_series(0,5) as dir;
+$$;
+
+select ap.xy2ijk(1,2);
+select ap.ijk2xy(ap.xy2ijk(1,2));
+
+select ap.ijk2xy(ap.adjhex(dir, ap.xy2ijk(1,2))) from generate_series(0,5) as dir;
+select ap.cantor(ap.ijk2xy(ap.adjhex(dir, ap.xy2ijk(1,2)))) from generate_series(0,5) as dir;
+
+
+\q
+
+-- build an adjacency list for the hexes
+create table ap.hexdata as
+select
+
+
+create table ap.population (
+ integer hex;
--- /dev/null
+\set scale 20000.0
+
+set search_path to ap,gis,public;
+
+\timing on
+
+drop table if exists images;
+create table images (
+ scenario text,
+ run integer,
+ gen integer,
+ png bytea
+);
+create index images_gen_index on images(gen);
+create index images_run_index on images(run);
+create index images_scen_index on images(scenario);
+create index images_srg_index on images(scenario,run,gen);
+
+drop table if exists tiles;
+create table tiles (
+ series text,
+ r raster
+);
+
+-- create a black basemap of all the africa hexes
+insert into tiles
+select
+'base',
+st_union(st_asraster(
+ P.geom, 20000.0, -20000.0, 0,0,ARRAY['8BUI','8BUI','8BUI','8BUI'],
+ -- next is constructing the color
+ -- yellow to blue or so is hue .16 - .66, sat .6 val 1.0
+ ARRAY[0,0,0,0],
+ --rgbb(rgb(ARRAY[0.0/2.0+0.16, 0.6*0.0, 1.0*0.0])) || ARRAY[0],
+ --rgbb(ARRAY[0,1,0])||ARRAY[0],
+ ARRAY[255,255,255,255], 0, 0, true
+))
+from afrhex P
+where not exists (select * from tiles where series = 'base')
+;
+
+
+create temp table heximg as
+select row_number() over () as id,
+1::integer as gen,
+st_union(st_asraster(
+ H.geom, 20000.0,-20000.0,0,0,ARRAY['8BUI', '8BUI', '8BUI', '8BUI'], C.color || ARRAY[0], ARRAY[0,0,0,0], 0, 0, true
+ --H.hex, 20000.0,-20000.0,0,0,ARRAY['8BUI', '8BUI', '8BUI', '8BUI'], ARRAY[0,255,0,0], ARRAY[0,0,0,0], 0, 0, true
+)) as image
+from afrhex H
+left join climate C on H.climate = C.climate
+;
+
+insert into images (scenario, run, gen, png)
+select
+'climate', 0, 0, st_aspng(HI.image,ARRAY[1,2,3]) as png
+from heximg HI
+;
+
+create or replace function merge(cover raster) returns raster language 'sql' as $$
+select st_union(x) from (
+ select r as x from tiles where series = 'base'
+ union all
+ select cover as x
+) as z
+;
+$$;
--- /dev/null
+\set scale 20000.0
+
+set search_path to ap,gis;
+
+\timing on
+
+with info as (
+ SELECT P.scenario, P.run, P.gen,
+ c.capacity,
+ P.climate,
+ p.freq,
+ a.geom,
+ P.freq[1] AS popa,
+ r.vec_sum(p.freq) - p.freq[1] AS popb,
+ CASE
+ WHEN p.freq[1] = 0 THEN 0.0::double precision
+ ELSE p.freq[1]::double precision / r.vec_sum(p.freq)::double precision
+ END AS popafrac,
+ r.vec_sum(p.freq) AS totalpop,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxfrac,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (PARTITION BY p.gen ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxgenfrac,
+ LEAST(1.0::double precision, r.vec_sum(p.freq)::double precision / c.capacity::double precision) AS pressure
+ FROM popgen p
+ JOIN climate c ON c.climate = p.climate
+ JOIN afrhex a ON a.hex = p.hex
+ left join images I on I.scenario = P.scenario and I.run = P.run and I.gen = P.gen
+ where I.run is null and I.run is null and I.gen is null
+ --where P.gen not exists (select 1 gen from images where scenario = :'scen' and gen = P.gen and run = )
+)
+insert into images (scenario, run, gen, png)
+select
+P.scenario, P.run, P.gen,
+st_aspng(
+-- merge(
+ st_union(
+ st_asraster(
+ P.geom, 20000.0, -20000.0, 0,0,ARRAY['8BUI','8BUI','8BUI','8BUI'],
+ -- next is constructing the color
+ -- yellow to blue or so is hue .16 - .66, sat .6 val 1.0
+ rgbb(rgb(ARRAY[
+ -- (1.6 / (1 + exp(-( P.popafrac - 0.5 )))-0.5) / 2.2 + 0.15, -- sigmoid bias to ends of hue
+ (4 / (1 + exp(-(P.popafrac - 0.5 ))) - 1.5) / 2.0 + 0.075,
+ -- P.popafrac/2.0+0.16, -- hue
+ --0.9*(P.pressure), -- saturation
+ 1.0, -- saturation
+ --1.0*(P.maxgenfrac) -- value
+ case when P.totalpop > 0 then
+ -- clamp max to 1.0, just in case the maxgen frac calculation was wrong
+ least(1.0, 0.15+0.85*sqrt(P.maxgenfrac)) -- value
+ --sqrt(P.pressure)
+ else
+ 0.0
+ end
+ --1.0 -- value
+ ]
+ )) || ARRAY[0],
+ --rgbb(rgb(ARRAY[0.0/2.0+0.16, 0.6*0.0, 1.0*0.0])) || ARRAY[0],
+ --rgbb(ARRAY[0,1,0])||ARRAY[0],
+ ARRAY[255,255,255,0], 0, 0, true
+ )
+ )
+-- )
+ ,ARRAY[1,2,3]
+) as png
+from info P
+--and P.totalpop > 0
+group by P.scenario, P.run, P.gen
+;
--- /dev/null
+drop schema if exists r cascade;
+create schema r;
+
+grant usage on schema r to public;
+
+set search_path to r;
+
+create or replace function rbinom(n integer, size integer, prob float) returns integer language 'plr' as $$
+return(rbinom(n,size,prob))
+$$;
+create or replace function setseed(n integer) returns integer language 'plr' as $$
+return(set.seed(n))
+$$;
+
+create or replace function rmultinom(n integer, size integer, prob float[]) returns integer[] language 'plr' as $$
+return(t(rmultinom(n,size,prob)))
+$$;
+
+-- a helper function to just return the single array
+create or replace function rmultinom(size integer, prob float[]) returns integer[] language 'plr' as $$
+return(t(rmultinom(1,size,prob))[1,])
+$$;
+
+create or replace function rmultinom(size integer, buckets integer) returns integer[] language 'plr' as $$
+return(t(rmultinom(1,size,rep(1, times=buckets)))[1,])
+$$;
+
+create or replace function rmhypergeo(n integer, v integer[]) returns integer[] language 'plr' as $$
+require(BiasedUrn);
+return(rMFNCHypergeo(1, v, n, rep(1, times=length(v))));
+$$;
+
+create or replace function vec_sum(n integer[]) returns integer language 'plr' as $$
+return(sum(n))
+$$ immutable strict;
+
+create or replace function vec_add(n integer[], b integer[]) returns integer[] language 'plr' as $$
+return(n+b)
+$$ immutable strict;
+
+create or replace function vec_add(v integer[], n integer) returns integer[] language 'plr' as $$
+return(v+n)
+$$ immutable strict;
+
+create or replace function vec_sub(n integer[], b integer[]) returns integer[] language 'plr' as $$
+return(n-b)
+$$ immutable strict;
+
+create or replace function vec_sub(v integer[], n integer) returns integer[] language 'plr' as $$
+return(v-n)
+$$ immutable strict;
+create or replace function vec_sub(n integer, v integer[]) returns integer[] language 'plr' as $$
+return(n-v)
+$$ immutable strict;
+
+-- can these be 'anyarray'?
+create or replace function sample(v text[], n integer) returns text[] language 'plr' as $$
+sample(v, n)
+$$ immutable strict;
+create or replace function sample(v numeric[], n integer) returns text[] language 'plr' as $$
+sample(v, n)
+$$ immutable strict;
+create or replace function sample(v text[], n integer, p integer[]) returns text[] language 'plr' as $$
+sample(v, n, prob = p)
+$$ immutable strict;
+create or replace function sample(v numeric[], n integer, p integer[]) returns text[] language 'plr' as $$
+sample(v, n, prob = p)
+$$ immutable strict;
+
+create aggregate vec_add(integer[]) (
+ sfunc = r.vec_add,
+ stype = _int4
+);
--- /dev/null
+-- set model externally
+
+\set scen largepop
+
+set search_path to ap,r,gis;
+
+update climate set capacity = basecapacity * (select st_area(geom)/1000000/100 from afrhex limit 1) * 2;
+
+truncate popgen;
+delete from images where series = 'freq';
+analyze images;
+
+-- create initial population
+insert into popgen
+select :'scen', 0, 0, H.hex, 0, ARRAY[0,0], H.climate
+from afrhex H
+;
+
+update popgen set freq = ARRAY[8000,8000]
+where hex = 4788
+and scenario = :'scen'
+and gen = 0
+and run = 0
+;
+
+\q
+
+update popgen set freq = ARRAY[16000,0]
+where climate = 'Aw'
+and (hex.xy(hex))[2] <= 40
+;
+
+update popgen set freq = ARRAY[0,16000]
+where climate = 'Aw'
+and (hex.xy(hex))[2] > 40
+;
--- /dev/null
+-- set model externally
+
+\set scen smallpop
+
+set search_path to ap,r,gis;
+
+update climate set capacity = basecapacity * (select st_area(geom)/1000000/100 from afrhex limit 1) * 2 / 4;
+
+-- TODO create a setup/initialization function that can do the basic
+-- setup so this file can concentrate on the scenario specific data
+delete from popgen where scenario = :'scen' and run = 0;
+delete from images where scenario = :'scen' and run = 0;
+analyze images;
+
+-- create initial population
+insert into popgen (scenario, run, gen, hex, freq, climate)
+select :'scen', 0, 0, H.hex, ARRAY[0,0], H.climate
+from afrhex H
+;
+
+update popgen set freq = ARRAY[2000,2000]
+where hex = 4788
+and scenario = :'scen'
+and gen = 0
+and run = 0
+;
+
+\q
+
+update popgen set freq = ARRAY[16000,0]
+where climate = 'Aw'
+and (hex.xy(hex))[2] <= 40
+;
+
+update popgen set freq = ARRAY[0,16000]
+where climate = 'Aw'
+and (hex.xy(hex))[2] > 40
+;
--- /dev/null
+#!/usr/bin/env perl
+
+use strict;
+use warnings;
+
+use DBI;
+use DBD::Pg qw(:pg_types);
+
+my $baseyear = -140000;
+my $genlen = 25;
+
+# TODO use environment variable for font
+my $font = '/Users/nw/Library/Fonts/SourceCodePro-Medium.ttf';
+
+my ($scenario, $run) = @ARGV;
+
+#-draw "text 10,220 'hello'"
+
+my $db = DBI->connect('dbi:Pg:service=afpopgen');
+$db->{AutoCommit} = 0;
+$db->{RaiseError} = 1;
+
+chdir "images/$scenario/$run" or die "can't chdir to images/$scenario/$run: $!\n";
+
+opendir(my $dir, '.');
+my @files = grep(/freq\d+\.png/, readdir $dir);
+close $dir;
+
+my $lookup = $db->prepare('select G.gen, sum(r.vec_sum(G.freq))/2 as totalpop from ap.popgen G where G.series = ? and G.run = ?');
+
+while (my ($gen,$pop) = $lookup->fetchrow_array) {
+#foreach my $img (@files) {
+ my $img = sprintf("%s%d-%04d.png", $scenario, $run, $gen);
+ next unless -f $img; # skip missing TODO warn, or just DL
+
+ my $x = 5; # where to put the annotation
+ my $y = 220;
+ my @cmd = ('convert', "$img", '-pointsize', '15', '-font', $font);
+ my $text = sprintf(q{text %d,%d 'gen: %8d'}, $x, $y, $gen);
+ push @cmd, '-draw', $text; $y += 15;
+
+ $lookup->execute($gen);
+ my ($sum) = $lookup->fetchrow_array();
+ $text = sprintf(q{text %d,%d 'pop: %8d'}, $x, $y, $sum);
+ push @cmd, '-draw', $text; $y += 15;
+ warn "annotating freq ", $gen, ' pop: ', $sum, "\n";
+ $text = sprintf(q{text %d,%d 'ybp: %8d'}, $x, $y, $baseyear+25*$gen);
+ push @cmd, '-draw', $text; $y += 15;
+ $text = sprintf(q{text %d,%d 'ams: %8d'}, $x, $y, $genlen*$gen);
+ push @cmd, '-draw', $text; $y += 15;
+
+ my $annfile = sprintf('ann-%s', $img);
+ system(@cmd, $annfile);
+ die "convert: $!\n" if $?;
+}
+$lookup->finish;
+
+$db->disconnect;
--- /dev/null
+#!/usr/bin/env perl
+
+use strict;
+use warnings;
+
+use DBI;
+use DBD::Pg qw(:pg_types);
+
+use File::Find;
+use File::Basename;
+use File::Slurp;
+
+my $baseyear = -140000;
+my $genlen = 25;
+my $font = '/Users/nw/Library/Fonts/SourceCodePro-Medium.ttf';
+
+my $cmd = shift;
+
+my $db = DBI->connect('dbi:Pg:service=afpopgen');
+$db->{AutoCommit} = 0;
+$db->{RaiseError} = 1;
+
+#$db->do('set search_path to dms,public');
+
+#$addfile = $db->prepare(q{insert into file (content, mtime, title) values (?, 'epoch'::timestamptz + ? * '1 second'::interval, ?) returning sha1});
+#$addfile->bind_param(1, undef, {pg_type => PG_BYTEA});
+#$addpath = $db->prepare(q{insert into paths (hash, path) values (?, ?)});
+#$listall = $db->prepare(q{select sha1, length(content) as size, mtime, title from file order by title});
+
+#my ($hash,$name) = (@ARGV);
+#die "no name specified. usage: dms extract <hash> <name>\n" unless defined $name;
+#die "file exists: $name\n" if -e $name;
+
+sub mkdirp {
+ my $path = '.';
+ foreach my $comp (@_) {
+ $path .= "/$comp";
+ mkdir $path unless -d $path;
+ die "can't make director $path: $!" unless -d $path;
+ }
+}
+
+my $lookup = $db->prepare(q{
+ select sum(r.vec_sum(G.freq))/2 as totalpop
+ from ap.popgen G where G.scenario = ? and G.run = ? and G.gen = ?
+});
+
+
+my $c = $db->prepare(q{select scenario,run,gen,length(png),png from ap.images});
+$c->execute();
+while (my ($scenario,$run,$gen,$size, $data) = $c->fetchrow_array) {
+ mkdirp('images',$scenario,$run);
+ my $name = sprintf("images/$scenario/$run/%s%d-%04d.png", $scenario,$run,$gen);
+ next if -e $name;
+
+ # get the population data
+ $lookup->execute($scenario, $run, $gen);
+ my ($pop) = $lookup->fetchrow_array();
+
+ open EXT, '>', $name;
+ print EXT $data;
+ close EXT;
+ print "extracted to $name\n";
+
+ # now annotate the image
+ my $x = 5; # where to put the annotation
+ my $y = 220;
+ my @cmd = ('convert', $name, '-pointsize', '15', '-font', $font);
+ my $text = sprintf(q{text %d,%d 'gen: %8d'}, $x, $y, $gen);
+ push @cmd, '-draw', $text; $y += 15;
+ $text = sprintf(q{text %d,%d 'pop: %8d'}, $x, $y, $pop);
+ push @cmd, '-draw', $text; $y += 15;
+ warn "annotating freq ", $gen, ' pop: ', $pop, "\n";
+ $text = sprintf(q{text %d,%d 'ybp: %8d'}, $x, $y, $baseyear+25*$gen);
+ push @cmd, '-draw', $text; $y += 15;
+ $text = sprintf(q{text %d,%d 'ams: %8d'}, $x, $y, $genlen*$gen);
+ push @cmd, '-draw', $text; $y += 15;
+
+ my $annfile = sprintf("%s-ann", $name);
+ warn "converting: @cmd\n";
+ system(@cmd, $annfile);
+ die "convert: $!\n" if $?;
+ rename $annfile, $name;
+}
+$c->finish;
+$lookup->finish;
+$db->disconnect;
--- /dev/null
+#!/usr/bin/env perl
+
+# This program connects to the database and creates one image at a time in a loop
+# It can thus be run in parallel
+
+use strict;
+use warnings;
+
+use DBI;
+use DBD::Pg qw(:pg_types);
+
+my $cmd = shift;
+
+my $db = DBI->connect('dbi:Pg:service=afpopgen');
+$db->{AutoCommit} = 0;
+$db->{RaiseError} = 1;
+
+#$db->do('set search_path to dms,public');
+
+#$addfile = $db->prepare(q{insert into file (content, mtime, title) values (?, 'epoch'::timestamptz + ? * '1 second'::interval, ?) returning sha1});
+#$addfile->bind_param(1, undef, {pg_type => PG_BYTEA});
+#$addpath = $db->prepare(q{insert into paths (hash, path) values (?, ?)});
+#$listall = $db->prepare(q{select sha1, length(content) as size, mtime, title from file order by title});
+
+#my ($hash,$name) = (@ARGV);
+#die "no name specified. usage: dms extract <hash> <name>\n" unless defined $name;
+#die "file exists: $name\n" if -e $name;
+
+sub mkdirp {
+ my $path = '.';
+ foreach my $comp (@_) {
+ $path .= "/$comp";
+ mkdir $path unless -d $path;
+ die "can't make director $path: $!" unless -d $path;
+ }
+}
+
+my $findun = $db->prepare(q{
+ select distinct P.scenario, P.run, P.gen
+ from ap.popgen P
+ left join images I on I.scenario = P.scenario and I.run = P.run and I.gen = P.gen
+ where I.scenario is null
+ });
+
+my $trylock = $db->prepare(q{select pg_try_advisory_lock(?,?)});
+my $searchpath = $db->prepare(q{set search_path to ap,gis,r});
+$searchpath->execute;
+my $makeimg = $db->prepare(q{
+with info as (
+ SELECT P.scenario, P.run, P.gen,
+ c.capacity,
+ P.climate,
+ p.freq,
+ a.geom,
+ P.freq[1] AS popa,
+ r.vec_sum(p.freq) - p.freq[1] AS popb,
+ CASE
+ WHEN p.freq[1] = 0 THEN 0.0::double precision
+ ELSE p.freq[1]::double precision / r.vec_sum(p.freq)::double precision
+ END AS popafrac,
+ r.vec_sum(p.freq) AS totalpop,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxfrac,
+ COALESCE(r.vec_sum(p.freq)::double precision / max(r.vec_sum(p.freq)) OVER (PARTITION BY p.gen ROWS BETWEEN UNBOUNDED PRECEDING AND UNBOUNDED FOLLOWING)::double precision, 0.0::double precision) AS maxgenfrac,
+ LEAST(1.0::double precision, r.vec_sum(p.freq)::double precision / c.capacity::double precision) AS pressure
+ FROM popgen p
+ JOIN climate c ON c.climate = p.climate
+ JOIN afrhex a ON a.hex = p.hex
+ where P.scenario = ?
+ and P.run = ?
+ and P.gen = ?
+)
+insert into images (scenario, run, gen, png)
+select
+P.scenario, P.run, P.gen,
+st_aspng(
+-- merge(
+ st_union(
+ st_asraster(
+ P.geom, 20000.0, -20000.0, 0,0,ARRAY['8BUI','8BUI','8BUI','8BUI'],
+ -- next is constructing the color
+ -- yellow to blue or so is hue .16 - .66, sat .6 val 1.0
+ rgbb(rgb(ARRAY[
+ -- (1.6 / (1 + exp(-( P.popafrac - 0.5 )))-0.5) / 2.2 + 0.15, -- sigmoid bias to ends of hue
+ (4 / (1 + exp(-(P.popafrac - 0.5 ))) - 1.5) / 2.0 + 0.075,
+ -- P.popafrac/2.0+0.16, -- hue
+ --0.9*(P.pressure), -- saturation
+ 1.0, -- saturation
+ --1.0*(P.maxgenfrac) -- value
+ case when P.totalpop > 0 then
+ -- clamp max to 1.0, just in case the maxgen frac calculation was wrong
+ least(1.0, 0.15+0.85*sqrt(P.maxgenfrac)) -- value
+ --sqrt(P.pressure)
+ else
+ 0.0
+ end
+ --1.0 -- value
+ ]
+ )) || ARRAY[0],
+ --rgbb(rgb(ARRAY[0.0/2.0+0.16, 0.6*0.0, 1.0*0.0])) || ARRAY[0],
+ --rgbb(ARRAY[0,1,0])||ARRAY[0],
+ ARRAY[255,255,255,0], 0, 0, true
+ )
+ )
+-- )
+ ,ARRAY[1,2,3]
+) as png
+from info P
+group by P.scenario, P.run, P.gen
+;
+});
+
+while (1) {
+ $findun->execute();
+ while (my ($sc,$r,$g) = $findun->fetchrow_array()) {
+ warn "trying $sc $r $g\n";
+ $trylock->execute($r,$g);
+ my ($lock) = $trylock->fetchrow_array();
+ next unless $lock;
+ warn "making $sc $r $g\n";
+ $makeimg->execute($sc,$r,$g);
+ $db->commit();
+ }
+ last unless $findun->rows();
+}
+
+$findun->finish;
+$trylock->finish;
+$makeimg->finish;
+$db->disconnect;
--- /dev/null
+Title: Africa Genetics
+CSS: layout.css
+
+Model
+-----
+
+Africa is divided up into hexes of 100 km across, or about 8600 sq.km.
+Each hex has a fixed climate, which is then assigned a carrying capacity.
+Each generation, the population either grows by 15% if it is below
+the carrying capacity, or drops by 40% if it is above carrying capacity.
+
+After the population growth, if the population is above capacity, the
+excess population moves to an adjacent hex, randomly chosen, but
+weighted by the available capacity of the adjacent hex to a minimum
+of the square root of the potential capacity.
+
+Population growth, reduction, and migration are all selected randomly
+from the available population.
+
+[https://granicus.if.org/projects/afpopgen/](Source Code) is available
+via git repository.
+
+African Colonization
+--------------------
+
+I have two animations set up.
+The model is then run for 1000 generations.
+
+There is no selection in this model, and the initial population
+is split 50-50 for each allele.
+
+Each frame is annotated with the generation number, the total
+African population, the year before present (starting at -140000),
+and the 'model year'. The colors range from orange to blue
+depending on the percentage of each allele present, the brightness
+of the colors are proportional to the population in each hex,
+ranging from 25% to 100% brightness when there is any population,
+and black when there is no population at all.
+
+### Large Capacity Model
+
+This model starts with a population of 8000 in a single hex in eastern Africa.
+The carrying capacity of each hex is set to estimates of modern hunter-gatherer
+population carrying capacities.
+
+[African Colonization Animation](largecap.mp4)
+
+ 
+
+### Small Capacity Model
+
+This model starts with a population of 2000 in a single hex in eastern Africa.
+The carrying capacity of each hex is set to one fourth of the large capacity model.
+
+[African Small Cap Animation](smallcap.mp4)
+
+ 
projects / afpopgen / commitdiff