#call this with: #read(cat(kernelopts(homedir), "/Dropbox/maple/init.mpl")); #external libraries: libname:=cat(kernelopts(homedir), "/Dropbox/maple/custom_libs"), libname; hostname:=ssystem("hostname");hostname:=StringTools[Trim](hostname[2]): #extract nth columns/rows from a matrix #these only work if have a 2d object... should be updated to also work #with 1d row/column vectors #Example use cases #A := LinearAlgebra:-RandomMatrix(20, 20, outputoptions = [datatype = float[8]]); #Ex: #nthColumns(A, 2); #Every other column #nthRows(A, 10)[.., 1..3]; #Every 10th row, but show only first 3 columns nthColumns:=proc(v_m, v_n) v_m[..,[seq(i, i=1..rtable_size(v_m)[2], v_n)]] end: nthRows:=proc(v_m, v_n) v_m[[seq(i, i=1..rtable_size(v_m)[1], v_n)],..] end: #nthElem:=proc(v_m, v_n) #v_m[[seq(i, i=1..rtable_size(v_m)[1], v_n)]] #end: #extracts data from an already created plot in a reasonable form #returns a list of matrices with data. Each matrix dim x 2 (or 3 in case of 3d plot) extractPlotData:=proc(v_plot, transp:=true) if transp then [seq(LinearAlgebra:-Transpose(i[3]), i=plottools:-getdata(v_plot, 'element' = "curve"))]; else [seq((i[3]), i=plottools:-getdata(v_plot, 'element' = "curve"))]; end; end: #example #a:=display([... plots ...]): #rearrangeCurves(a, [-1,3]): #rearrangeCurves:=proc(v_items, v_reorder:=[]) ##Reorder should be a list of index pairs like ##[[3,5], [-1, 1], ...] #local p, temp, curves:=[], rest:=[]: ##separate curves from rest (there must be a prettier way to do this) #for p in convert(v_items, list) do #if type(p,'function') and op(0,p) = ('CURVES') then #curves:=[curves[], p]: #else #rest:=[rest[], p]: #end; #end; ##now reorder #for p in v_reorder do #temp:=curves[p[1]]: #curves[p[1]]:=curves[p[2]]: #curves[p[2]]:=temp: #end; #PLOT(curves[], rest[]): #end: #rearrange curves inside a plot --- so that one can choose which curves is on top/bottom even after #the plot has been created #discussed here: #http://www.mapleprimes.com/questions/201626-Order-Of-Curves-In-A-Plot rearrangeCurves:= proc( v_items::specfunc(anything, PLOT), v_reorder::list([integer,integer]):= [] ) local p, curves, rest; (curves,rest):= selectremove(type, v_items, specfunc(anything, CURVES)); curves:= < op(curves) >: for p in v_reorder do curves[p]:= curves[p[[2,1]]] end do; PLOT(convert(curves,list)[], op(rest)) end proc: #for numerical differentiation #based on the idea from: #http://www.mapleprimes.com/posts/119554-Data-Interpolation #example use: #alist:=[seq(i, i=0..10, 0.1)]: #data:=map(x->evalf(sin(x)), alist): #plot(alist, data); #plot([cos(x), 'num_diff(x, LinearAlgebra:-Transpose(Matrix([alist,data])))'], x=1..10, thickness=5, linestyle=[solid, dot], color=[blue, red]); num_diff:=proc(x, v_data, v_options:=[method=spline, degree=3]) #v_data is a matrix #TODO: let the data be in a more arbitrary format that ArrayInterpolation understands, but keep x as first var evalf(D[1](x->CurveFitting:-ArrayInterpolation(v_data, [x],v_options[])[])(x)); end: #saves a png plot savePlot:=proc(v_p, v_fileName, v_w:="800", v_h:="500") plotsetup("png", plotoutput=v_fileName, plotoptions=cat("quality=100,portrait,noborder,width=",v_w,",height=",v_h)); print(plots[display](v_p)); Threads[Sleep](2): fclose(v_fileName): plotsetup(default); end: #this is useful when need to write data (a list of lists) with potentially #different lengths. For example when a few different curves are extracted #from a plot. #NOTE: writedata seems to always write "\n" for a new line... #could update it to handle unix new lines. See code with showstat(writedata); exportData:=proc(v_fileName, v_data) local fd, item: fd := fopen(v_fileName, WRITE, TEXT): for item in v_data do #writedata(fd, convert(item[1..2, 1..5], listlist) ,integer,proc(f,x) fprintf(f,`%10f`, x) end proc); writedata(fd, convert(item, listlist) ,float); end; close(fd); end: plot_options1:=size=[1000, 1000/1.61], font=[TIMES, roman, 32], labelfont=[TIMES, roman, 32], axesfont=[TIMES, roman, 32], axes=boxed; #Assumes file name ends in ".png" #size directive to plot or display should be respected by this code savePlotPNGandPDF:=proc(v_p, v_fileName) plotsetup("png", plotoutput=v_fileName, plotoptions="quality=100,portrait,noborder"); print(plots[display](v_p)); Threads[Sleep](2): fclose(v_fileName): plotsetup(default); Threads[Sleep](1): ssystem(sprintf("convert %s %s.pdf",v_fileName, v_fileName[1..-5])); end: #TEMPORARY HACK savePlotAndData:=proc(v_p, v_fileName, v_w:="750", v_h:="450") savePlot(v_p, v_fileName, v_w, v_h): #this assumes all the curves in the plot have the same number of points #ExportMatrix(cat(v_fileName, ".txt"), ArrayTools:-Concatenate(2, extractPlotData(v_p)[]), mode = ascii) ; #GRRR... make this so that non 3 letter extensions are handled properly exportData(cat(v_fileName[..-5], ".txt"), extractPlotData(v_p)) ; end: #saveData:=proc(v_data), #Find the equations of motion for a given Lagrangian eoms_for_L:=proc(v_L, v_q) local i, eqs_list: eqs_list:=[]: for i from 1 to nops(v_q) do eqs_list:=[eqs_list[], diff( Physics:-diff(v_L, diff(v_q[i], t)), t) - Physics:-diff(v_L, v_q[i]) = 0]: end; eqs_list: end: #Find the equations of motion for a given Hamiltonian eoms_for_H:=proc(v_H, v_q, v_p) #assume v_q, v_p lists the same dimension local i, eqs_list: eqs_list:=[]: for i from 1 to nops(v_q) do eqs_list:=[eqs_list[], diff(v_q[i],t)=Physics:-diff(v_H, v_p[i]), diff(v_p[i],t)=-Physics:-diff(v_H, v_q[i])]: end; eqs_list: end: #Convert first order EOMs (derived from a Hamiltonian) to #second order ones H_eoms_to_L_eoms:=proc(v_H_eoms) #assumes order of equations as [q1,p1,q2,p2,...] #with q the coordinate and p the momentum local i, eqs_list, q_diff: eqs_list:=[]: q_diff:=seq(v_H_eoms[2*i] , i=1..nops(v_H_eoms)/2): for i from 1 to nops(v_H_eoms)/2 do eqs_list:=[eqs_list[], subs(q_diff, diff(v_H_eoms[2*i-1], t))]: end: eqs_list: end: #plotting: #creates a matrix of plots of size h x w from a list - can do display(p2m([p1, p2, ...], 2, 2)), etc #probably a bit broken p2m:=proc(v_d, w:=2, h:=0 ) local l_d, l_l, l_w, l_h: #if h is 0, just pack as many items as required if h=0 then l_h:=ceil(nops(v_d)/w): else l_h:=h: end; if nops(v_d)<2 then l_w:=1: else l_w:=w: end; l_l:=(l_h*l_w)-nops(v_d): if l_l<0 then print("Data too long by", abs(l_l), "...truncating"): #only grap plots that fit into the matrix l_d:=v_d[1..l_h*l_w]: else l_d:=[op(v_d), (plot([0,0], 0..0, axes=none))$l_l]: #pad with empty plots end if: Matrix(l_d, l_h, l_w); #ArrayTools:-Reshape(Matrix(l_d, h, w); end proc: #Show items of a container (i.e. list); one per line #this is useful when dealing with a list of equations that are being manipulated. #at the end of a few manipulation steps one can add #pp(%) #which will display them one-per-line pp:=proc(v_list) local i; for i in v_list do print(i); end; end: #get list of many random colors get_color_list:=proc(n::integer) [seq(ColorTools[Color]([evalf(i/n),0,1-evalf(i/n)]), i=1..n)]; end: #apply function term-by-term; there has to be something built-in to do this?! #examples use: #apply_each(simplify, long_expression, {exp=alt_exp}) apply_each:=proc(v_f, v_l) add(i, i=map(v_f, convert(v_l, list), _rest)); end: #Provides a means to define variables that can be easily replaced #by ones passed in via the command line #example use case: #in the file do: #rs:=fromInput(rs_INPUT, 500): #then run maple with #maple -c "rs_INPUT:=300:" mycode.mpl fromInput:=proc(v_value, v_default) if (StringTools[IsSuffix]("_INPUT", convert(v_value, string))) then return v_default: end; return v_value: end: #How many nodes should hackMap below use #depends on what host we run on optimal_num_nodes:=4: #default if hostname="desktop5" then optimal_num_nodes:=6: elif hostnam="laptop3" then optimal_num_nodes:=2: end: #FIX for a broken Grid:-Map - only works in maple 18 hackMap := proc(f, l, nnodes:=optimal_num_nodes) #hackMap := proc(f, l, nnodes:=4) local mine, r, i, m, thisNode, totalNodes, var, s, e, p, x; global `grid/mapcmd`; if not Grid:-Status("running") then #var := `union`(`union`(map(proc (x) op(0,x) end proc,indets(f,function)),indets(f,name)),{'`grid/mapcmd`'}); var := {anames('user'), '`grid/mapcmd`'}; var := remove(x -> type(x,protected),var); `grid/mapcmd` := _Inert_PROC(_Inert_GLOBALSEQ(seq(_Inert_NAME(convert(i,string)),i = var)),_Inert_STATSEQ(_Inert_FUNCTION(_Inert_TABLEREF(_Inert_NAME("Grid"),_Inert_EXPSEQ(_Inert_NAME("Map"))),_Inert_EXPSEQ(_Inert_UNEVAL(_Inert_VERBATIM(''f'')),_Inert_UNEVAL(_Inert_VERBATIM([args[2 .. (NULL)]],0)))))); try #return Grid:-Launch(proc () global `grid/mapcmd`; FromInert(`grid/mapcmd`)() end proc,imports = var) #print(nnodes); return Grid:-Launch(proc () global `grid/mapcmd`; FromInert(`grid/mapcmd`)() end proc,imports = var, numnodes=nnodes) finally `grid/mapcmd` := '`grid/mapcmd`' end try end if; thisNode := Grid:-MyNode(); totalNodes := Grid:-NumNodes(); try if l::rtable then s, e := getbounds(thisNode,totalNodes,numelems(l)); mine := map(f,l(s .. e),_rest); if thisNode = 0 then r := copy(l); r(s .. e) := mine; for i to totalNodes-1 do s, e, m := Grid:-Receive(); r(s .. e) := m end do; return r else Grid:-Send(0,s,e,mine) end if elif l::{list, set} then s, e := getbounds(thisNode,totalNodes,numelems(l)); mine := map(f,l[s .. e],_rest); if thisNode = 0 then return [op(mine), seq(op(Grid:-Receive(i)),i = 1 .. totalNodes-1)] else Grid:-Send(0,mine) end if elif l::{array, table} then p := sort([indices(l,'pairs')]); s, e := getbounds(thisNode,totalNodes,nops(p)); mine := map(x -> lhs(x) = f(rhs(x)),p[s .. e],_rest); if thisNode = 0 then m := table(); map2((m, i) -> assign(m[lhs(i)],rhs(i)),m,mine); for i to totalNodes-1 do map2((m, i) -> assign(m[lhs(i)],rhs(i)),m,Grid:-Receive()) end do; return eval(m) else Grid:-Send(0,mine) end if else s, e := getbounds(thisNode,totalNodes,nops(l)); mine := map(f,[op(l)[s .. e]],_rest); if thisNode = 0 then r := [op(mine), seq(op(Grid:-Receive(i)),i = 1 .. totalNodes-1)]; return op(0,l)(op(r)) else Grid:-Send(0,mine) end if end if catch : if thisNode <> 0 then Grid:-Send(0,[FAIL, StringTools:-FormatMessage(lastexception[2 .. (NULL)])]); error end if end try; NULL end proc: #show all the stuff we've defined sort([anames('user')]);