From 1b287becd194657f2db924d3603b4562115b5b77 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Tue, 11 Jul 2023 14:33:00 -0400 Subject: [PATCH] added microstates to trajectories --- mcmc_clustering_eap_chain.jl | 14 ++++++--- run/phases-random_2023-07-11.jl | 51 +++++++++++++++++++++++++++++++++ run/phases_random.jl | 50 ++++++++++++++++++++++++++++++++ 3 files changed, 111 insertions(+), 4 deletions(-) create mode 100644 run/phases-random_2023-07-11.jl create mode 100644 run/phases_random.jl diff --git a/mcmc_clustering_eap_chain.jl b/mcmc_clustering_eap_chain.jl index f2b0e83..e5f0b37 100644 --- a/mcmc_clustering_eap_chain.jl +++ b/mcmc_clustering_eap_chain.jl @@ -3,7 +3,7 @@ using Distributions; using LinearAlgebra; using Logging; using DelimitedFiles; -using ProfileView; +#using ProfileView; using DecFP; using Quadmath; @@ -235,7 +235,10 @@ function mcmc(nsteps::Int, pargs, chain::EAPChain) ]); outfile = open("$(pargs["prefix"])_trajectory.csv", "w"); - writedlm(outfile, ["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], ','); + writedlm(outfile, hcat(["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], + reshape(vcat(reshape(["phi$i" for i=1:n(chain)], 1, :), reshape(["theta$i" for i=1:n(chain)], 1, :)), 1, :), + reshape(vcat(reshape(["mux$i" for i=1:n(chain)], 1, :), reshape(["muy$i" for i=1:n(chain)], 1, :), reshape(["muz$i" for i=1:n(chain)], 1, :)), 1, :)), + ','); rollfile = open("$(pargs["prefix"])_rolling.csv", "w"); writedlm(rollfile, ["step" "r1" "r2" "r3" "r1sq" "r2sq" "r3sq" "rsq" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psq" "U" "Usq" "Ealign" "psi"], ','); @@ -293,7 +296,10 @@ function mcmc(nsteps::Int, pargs, chain::EAPChain) if step % pargs["stepout"] == 0 writedlm(outfile, hcat(step, transpose(chain.r), - transpose(chain_μ(chain)), chain.U), + transpose(chain_μ(chain)), chain.U, + reshape(vcat(transpose(chain.ϕs), transpose(chain.θs)), 1, :), + reshape(chain.μs, 1, :) + ), ','); writedlm(rollfile, hcat( @@ -332,7 +338,7 @@ end (chain, sas, vas, ar) = if pargs["profile"] @info "Profiling the mcmc code..."; mcmc(5, pargs); # run first to compile code - @profview mcmc(pargs["num-steps"], pargs); + #@profview mcmc(pargs["num-steps"], pargs); println("Press ENTER to continue..."); readline(stdin); exit(1); diff --git a/run/phases-random_2023-07-11.jl b/run/phases-random_2023-07-11.jl new file mode 100644 index 0000000..1f56aac --- /dev/null +++ b/run/phases-random_2023-07-11.jl @@ -0,0 +1,51 @@ +println("Hello world!"); + +@everywhere using Printf; +@everywhere using Distributions; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_kappa-$(fmt(case[:kappa]))"; +end + +@everywhere Kdist = collect(zip([1.0; 0.0], [0.0; 1.0], [2.0; 1.0], [1.0; 2.0], [3.0; 2.0], [2.0; 3.0])); +@everywhere Edist = Normal(2.5, 1.0); +@everywhere ndist = Normal(50, 15); +@everywhere bdist = Normal(1.0, 1.0); +@everywhere kTdist = Uniform(0.0, 2.0); +@everywhere Fdist = Normal(2.0, 1.0); + +mkpath(workdir); + +@sync @distributed for i=1:1000000 + + case = Dict(); + case[:kappa] = 0.0; + case[:E0] = max(0.05, rand(Edist)); + Ks = rand(Kdist) + case[:K1], case[:K2] = Ks[1], Ks[2] + case[:n] = max(10, round(Int, rand(ndist))); + case[:b] = rand(bdist); + case[:kT] = 10.0^(rand([-1; 1])*rand(kTdist)); + case[:Fx] = max(0.0, rand(Fdist)); + case[:Fz] = max(0.0, rand(Fdist)); + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --bend-mod $(case[:kappa]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2000000 --burn-in 200000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end diff --git a/run/phases_random.jl b/run/phases_random.jl new file mode 100644 index 0000000..28f449f --- /dev/null +++ b/run/phases_random.jl @@ -0,0 +1,50 @@ +println("Hello world!"); + +@everywhere using Printf; +@everywhere using Distributions; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_kappa-$(fmt(case[:kappa]))"; +end + +cases = Any[]; +κs = [0.0]; +Ks = zip([1.0; 0.0], [0.0; 1.0]); +E0s = vcat(0.0, 0.05:0.05:1.0, 1.5:0.5:5.0, 10.0); +ns = Int[100; 200]; +bs = [0.5; 1.0]; +kTs = [0.1; 1.0]; +Fs = zip([0.0; 1.0; 0.0; 1.0/sqrt(2.0)], [0.0; 0.0; 1.0; 1.0/sqrt(2.0)]); +for b in bs, κ in κs, n in ns, Kvec in Ks, E0 in E0s, kT in kTs, F in Fs + K1, K2 = Kvec; + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, + :kT => kT, :Fz => F[2], :kappa => κ, + :Fx => F[1], :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --bend-mod $(case[:kappa]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2000000 --burn-in 200000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases);