From: Nathan Wagner Date: Wed, 18 Feb 2015 14:05:27 +0000 (-0600) Subject: organized code X-Git-Url: https://granicus.if.org/projects?a=commitdiff_plain;p=afpopgen organized code --- diff --git a/Makefile b/Makefile index 8f322f4..dc5d8b3 100644 --- a/Makefile +++ b/Makefile @@ -1,3 +1,121 @@ +# 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 @@ -6,3 +124,4 @@ notes.pdf: notes.tex notes.bib clean: rm -f notes.pdf +.PHONY: annotate images diff --git a/model/adjlist.sql b/model/adjlist.sql new file mode 100644 index 0000000..2a476d1 --- /dev/null +++ b/model/adjlist.sql @@ -0,0 +1,29 @@ +--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; + diff --git a/model/africa.sql b/model/africa.sql new file mode 100644 index 0000000..6f3e756 --- /dev/null +++ b/model/africa.sql @@ -0,0 +1,122 @@ +--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; diff --git a/model/base.sql b/model/base.sql new file mode 100644 index 0000000..3de7532 --- /dev/null +++ b/model/base.sql @@ -0,0 +1,52 @@ +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 +; diff --git a/model/colors.sql b/model/colors.sql new file mode 100644 index 0000000..1233c77 --- /dev/null +++ b/model/colors.sql @@ -0,0 +1,44 @@ +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; diff --git a/model/dumpgenimg.sql b/model/dumpgenimg.sql new file mode 100644 index 0000000..dcf4f26 --- /dev/null +++ b/model/dumpgenimg.sql @@ -0,0 +1,68 @@ +\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 +; diff --git a/model/elev.sql b/model/elev.sql new file mode 100644 index 0000000..9fdc617 --- /dev/null +++ b/model/elev.sql @@ -0,0 +1,31 @@ +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 diff --git a/model/generation.sql b/model/generation.sql new file mode 100644 index 0000000..9def0da --- /dev/null +++ b/model/generation.sql @@ -0,0 +1,123 @@ +-- 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 diff --git a/model/hexfuncs.sql b/model/hexfuncs.sql new file mode 100644 index 0000000..8826e17 --- /dev/null +++ b/model/hexfuncs.sql @@ -0,0 +1,332 @@ +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; diff --git a/model/hexops.sql b/model/hexops.sql new file mode 100644 index 0000000..0fec49c --- /dev/null +++ b/model/hexops.sql @@ -0,0 +1,119 @@ +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; diff --git a/model/image.sql b/model/image.sql new file mode 100644 index 0000000..5fade37 --- /dev/null +++ b/model/image.sql @@ -0,0 +1,67 @@ +\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 +; +$$; diff --git a/model/mkoneimg.sql b/model/mkoneimg.sql new file mode 100644 index 0000000..dcf4f26 --- /dev/null +++ b/model/mkoneimg.sql @@ -0,0 +1,68 @@ +\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 +; diff --git a/model/rfuncs.sql b/model/rfuncs.sql new file mode 100644 index 0000000..5db2acc --- /dev/null +++ b/model/rfuncs.sql @@ -0,0 +1,73 @@ +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 +); diff --git a/notes.bib b/paper/notes.bib similarity index 100% rename from notes.bib rename to paper/notes.bib diff --git a/notes.tex b/paper/notes.tex similarity index 100% rename from notes.tex rename to paper/notes.tex diff --git a/scenarios/largepop.sql b/scenarios/largepop.sql new file mode 100644 index 0000000..77ec143 --- /dev/null +++ b/scenarios/largepop.sql @@ -0,0 +1,36 @@ +-- 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 +; diff --git a/scenarios/smallpop.sql b/scenarios/smallpop.sql new file mode 100644 index 0000000..8cb4ccc --- /dev/null +++ b/scenarios/smallpop.sql @@ -0,0 +1,38 @@ +-- 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 +; diff --git a/utils/annotate b/utils/annotate new file mode 100755 index 0000000..5cf060e --- /dev/null +++ b/utils/annotate @@ -0,0 +1,58 @@ +#!/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; diff --git a/utils/fetchimg b/utils/fetchimg new file mode 100755 index 0000000..e0cc957 --- /dev/null +++ b/utils/fetchimg @@ -0,0 +1,87 @@ +#!/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 \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; diff --git a/utils/mkimage b/utils/mkimage new file mode 100755 index 0000000..96ec8cf --- /dev/null +++ b/utils/mkimage @@ -0,0 +1,129 @@ +#!/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 \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; diff --git a/web/index.md b/web/index.md new file mode 100644 index 0000000..a069c01 --- /dev/null +++ b/web/index.md @@ -0,0 +1,57 @@ +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) + +![Generation 0](largecap0000.png) ![Generation 1000](largecap1000.png) + +### 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) + +![Generation 0](smallcap0000.png) ![Generation 1000](smallcap1000.png)