養花種魚數月亮賞星星

這是留言板-既然來了，就留個言再走吧！

2029-07-22T14:43:00.005+08:00

Please tell me what you thought about this blog. Comments, suggestions and views are welcomed.

Mathematica教學：Simplex method in tableau form

2025-08-06T07:41:00.009+08:00

myConstrains[constrains_] := 
 Block[{num, s, mycon}, num = Length[constrains];
  mycon = 
   Table[If[Head[constrains[[i]]] === LessEqual, 
     constrains[[i]][[1]] + s[i] == constrains[[i]][[2]], 
     If[Head[constrains[[i]]] === Equal, 
      constrains[[i]][[1]] + R[i] == constrains[[i]][[2]], 
      constrains[[i]][[1]] - s[i] + R[i] == 
       constrains[[i]][[2]]]], {i, num}];
  mycon /. x__[y_] :> ToExpression[ToString[x] <> ToString[y]]]

LPstd[obj_, mycons_] := 
 Block[{myvar, basicvars, nonbasicvars, AB, BM, c, A, B, rowindex1, 
   rowindex2}, myvar = Variables[mycons[[All, 1]]];
  myvar = 
   Flatten[{Sort@
      Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "R"] &]}];
  AB = {Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "R"] &]};
  AB = Transpose@Flatten[AB, 1];
  AB = Join[{Flatten[
      Table[Select[myvar, 
        StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]]}, 
    AB];
  c = Flatten[{Table[
      Coefficient[obj, #] & /@ 
       Select[myvar, StringStartsQ[SymbolName[#], i] &], {i, {"x", 
        "s"}}], Table[-BM, {Length[
        Select[myvar, StringStartsQ[SymbolName[#], "R"] &]]}], 0}];
  basicvars = 
   Select[AB[[All, #]] & /@ (Flatten /@ 
        Map[Position[myvar, #] &, 
         Tuples[myvar, {Length[mycons]}], {2}]), #[[2 ;; -1]] === 
       IdentityMatrix[Length[mycons]] &][[1, 1]];
  nonbasicvars = 
   Flatten[Table[
     Select[Complement[myvar, basicvars], 
      StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]];
  c = Flatten@{c[[Flatten[Position[myvar, #] & /@ nonbasicvars]]], 
     c[[Flatten[Position[myvar, #] & /@ basicvars]]], 0};
  AB = Join[{Flatten@{nonbasicvars, basicvars}}, 
    Flatten /@ 
     Transpose@{AB[[2 ;; -1]][[All, 
         Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
       IdentityMatrix[Length[mycons]]}];
  rowindex1 = c[[-Length[basicvars] - 1 ;; -2]];
  rowindex2 = Range[Length[mycons]];
  Flatten[{c - 
       Total[rowindex1[[#]]*
           Flatten[{AB[[2 ;; -1]][[rowindex2[[#]]]], 
             mycons[[rowindex2[[#]], 2]]}] & /@ 
         Range[Length@rowindex1]], 
      AB[[2 ;; -1]][[All, 
        Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
      IdentityMatrix[Length[mycons]], mycons[[All, 2]], 
      Flatten@{Sort@
         Select[nonbasicvars, StringStartsQ[SymbolName[#], "x"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "s"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "R"] &]}, basicvars} & /@ 
    Range[2], 1]]

simplexmethod[c_, A_, B_, b_, xnb_, xb_, cp_, Ap_, Bp_, bp_, xnbp_, 
  xbp_] := 
 Block[{cB, nA, nB, mysymbol, pivotCol, AB, leaveVarpos, Bnew, xbnew, 
   xnbnew, Anew, cnew, znew, z, BM, ABp},
   (*初始化變數*)
   nA = Length[xnb];
   nB = Length[xb];
   
  (*建立符號變數*)
  mysymbol = Flatten@{xnbp, xbp};

  (*找出進入變數（最負的檢驗數）*)
  pivotCol = PositionSmallest[-c[[1 ;; -2]] /. BM -> 10^6][[1]];
  
  (*合併A和B矩陣*)
  AB = Flatten /@ Transpose[{A, B}];
  Print@MatrixForm@AB;
  ABp = Flatten /@ Transpose[{Ap, Bp}];
  leaveVarpos = 
   PositionSmallest[
     b/(If[# <= 0, 10^-6, #] & /@ AB[[All, pivotCol]] /. BM -> 10^6) /. 
      x_ /; x <= 0 -> 10^6][[1]];
  
  (*更新基變數和非基變數*)
  xbnew = xb /. xb[[leaveVarpos]] -> mysymbol[[pivotCol]];
  xnbnew = xnb /. mysymbol[[pivotCol]] -> xb[[leaveVarpos]];
  Bnew = 
   ABp[[All, Flatten[Position[mysymbol, #] & /@ Flatten@{xbnew}]]];
  cB = cp[[Flatten[Position[mysymbol, #] & /@ xbnew]]];
  
  (*計算新的檢驗數和目標函數值*)
  znew = Flatten[{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
      cB . Inverse[Bnew], cB . Inverse[Bnew] . bp - cp[[-1]]}] // 
    FullSimplify;
  
  (*顯示單形法表格*)
  Print[Grid[
    Insert[Flatten /@ 
      Transpose[{Flatten[{z, xbnew}], 
        Insert[Map[Flatten, 
          Transpose[{Inverse[Bnew] . Ap, Inverse[Bnew], 
            Inverse[Bnew] . bp}]], znew, 1]}], 
     Flatten[{Item["BasicVariable", 
        Background -> Lighter[Yellow, 0.9]], mysymbol, "Solution"}], 
     1], Frame -> All, 
    Background -> {{Lighter[Gray, 0.8]}, {Lighter[Yellow, 0.9], 
       Lighter[Pink, 0.9]}}, ItemStyle -> Directive[FontSize -> 16], 
    Spacings -> {Table[5, {Length[mysymbol] + 3}], 2}]];
  
  (*返回更新後的參數*)
  {-Flatten@{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
       cB . Inverse[Bnew]}, Inverse[Bnew] . Ap, Inverse[Bnew], 
    Inverse[Bnew] . bp, xnbnew, xbnew, cp, Ap, Bp, bp, xnbp, xbp} // 
   FullSimplify]

Nest[simplexmethod @@ # &, 
 LPstd[-15 x1 - 20 x2, 
  myConstrains[{x1 + 2 x2 >= 10, 2 x1 - 3 x2 <= 6, x1 + x2 >= 6}]], 2]

Mathematica教學：Solving the Traveling Salesman Problem using Gurobi solver

2025-07-31T12:30:00.004+08:00

You can download Gurobi solver from https://www.gurobi.com/academia/academic-program-and-licenses/ .

(* Step 1: Register for an academic license on the Gurobi website and install the Gurobi engine. *)

(* Step 2: Append the Gurobi license path if needed *)
AppendTo[$Path, "your full path to the gurobi.lic directory"]

(* Step 3: Restart Mathematica before continuing *)

(* Generate test data: random 2D points *)
{pts, dist, n} = {#, DistanceMatrix[#], Length[#]} &@RandomReal[{0, 100}, {40, 2}];

(* Define decision variables *)
{xVars, uVars} = {
  Flatten@Table[x[i, j], {i, n}, {j, n}] // DeleteCases[x[i_, i_]],
  Table[u[i], {i, 2, n}]
};

(* Objective function: minimize total distance *)
obj = Total[Flatten[Table[x[i, j], {i, n}, {j, n}]*dist]];

(* Flow constraints: each city must be entered and exited exactly once *)
flow = Flatten@{
  Table[Sum[If[i != j, x[i, j], 0], {j, n}] == 1, {i, n}],
  Table[Sum[If[i != j, x[j, i], 0], {j, n}] == 1, {i, n}]
};

(* Miller-Tucker-Zemlin (MTZ) subtour elimination constraints *)
MTZ=Flatten[Table[If[i!=j&&i>1&&j>1,u[i]-u[j]+n*x[i,j]<=n-1,
Nothing],{i,n},{j,n}]];

(* Variable domains and integrality conditions *)
vardomain=Flatten@{
Thread[0<=xVars<=1],
u[1]==1,Thread[2<=uVars<=n],
xVars∈Integers,uVars∈Integers};

(* Solve the model using Gurobi *)
{optDist,sol}=NMinimize[
Flatten@{obj,flow,MTZ,vardomain},
Join[xVars,uVars],
Method->"Gurobi"];

(* Extract path from the binary decision variables *)
edges=Cases[sol,(x[i_,j_]->1):>{i,j}];
path=NestList[
Cases[edges,{#[[2]],_}][[1]]&,
edges[[1]],
Length@edges-1];

(* Visualization of the solution *)
Graphics[{
  {Red, PointSize[Large], Point@pts},
  {Blue, Thick, Line[pts[[#]] & /@ path]},
  {Green, PointSize[.02], Point@pts[[1]]},
  {Arrowheads[.02], Arrow[pts[[#]] & /@ path]}},
 Frame -> True,
 PlotLabel -> Row[{"TSP: n=", n, ", distance=", optDist}]
]

Mathematica 教學：Searoute for Traveling Salesman Problem

2025-07-29T10:01:00.007+08:00

You can download searoute from https://github.com/eurostat/searoute .

port = {
 {"Shanghai, China", 31.2304, 121.4737},
 {"Guangzhou, China", 23.1291, 113.2644},
 {"Rotterdam, Netherlands", 51.9225, 4.4792},
 {"Antwerp, Belgium", 51.2194, 4.4025},
 {"Valencia, Spain", 39.4699, -0.3763},
 {"Los Angeles, USA", 33.7405, -118.2719},
 {"New York/New Jersey, USA", 40.6687, -74.0451},
 {"Savannah, USA", 32.0809, -81.0912},
 {"Santos, Brazil", -23.9608, -46.3336},
 {"Buenos Aires, Argentina", -34.6037, -58.3816},
 {"Itajai, Brazil", -26.9101, -48.6703},
 {"Tangier Med, Morocco", 35.8844, -5.5036},
 {"Port Said, Egypt", 31.2653, 32.3019},
 {"Durban, South Africa", -29.8587, 31.0218},
 {"Melbourne, Australia", -37.8136, 144.9631},
 {"Auckland, New Zealand", -36.8485, 174.7633},
 {"Brisbane, Australia", -27.4698, 153.0251},
 {"Singapore, Singapore", 1.3521, 103.8198},
 {"Kaohsiung, Taiwan", 22.6273, 120.3014},
 {"Busan, South Korea", 35.1796, 129.0756},
 {"Hong Kong, China", 22.3193, 114.1694},
 {"Tanjung Pelepas, Malaysia", 1.3673, 103.5386},
 {"Jebel Ali, UAE", 25.015, 55.0633},
 {"Colombo, Sri Lanka", 6.9271, 79.8612},
 {"Felixstowe, UK", 51.9556, 1.3111},
 {"Algeciras, Spain", 36.1408, -5.4562},
 {"Vancouver, Canada", 49.2827, -123.1207},
 {"Port Klang, Malaysia", 3.0089, 101.3678},
 {"Laem Chabang, Thailand", 13.0846, 100.8956},
 {"Tokyo, Japan", 35.6828, 139.7595},
 {"Yokohama, Japan", 35.4437, 139.6380}
};

Export["/Users/yourname/Desktop/searoute/text.csv",
Insert[
Flatten@{#[[1,1]]<>" to "<>#[[2,1]],
Round[#[[1,3]],0.1],Round[#[[1,2]],0.1],
Round[#[[2,3]],0.1],Round[#[[2,2]],0.1]}&/@
Select[
Flatten[Table[{port[[i]],port[[j]]},
{i,Length@port},{j,Length@port}],1],
#[[1]]!=#[[2]]&
],
{"routename","olon","olat","dlon","dlat"},1]
];

RunProcess[{"java","-jar",
"/Users/yourname/Desktop/searoute/searoute.jar",
"/Users/yourname/Desktop/searoute/text.csv",
"-res","5","-panama","1","-suez","1"}]

route=Import["/Users/yourname/Desktop/searoute/out.geojson",
"Data"][[2,2]];

mypath=Flatten@{#[[2,2]],StringSplit["routename"/.#[[3,2]]," to "],
"distKM"/.#[[3,2]]}&/@route;

myrule=Rule@@@Transpose@{Tally[mypath[[All,2]]][[All,1]],
Range[Length@port]};

distMatrix=ParallelTable[If[i==j,Infinity,
ToExpression@Select[mypath/.myrule,
#[[2]]==i&&#[[3]]==j&][[1,-1]]],
{i,Length@port},{j,Length@port}];

tour=FindShortestTour[Range[Length@port],
DistanceFunction->(distMatrix[[#1,#2]]&)];

myline=ParallelTable[
If[i==j,{},Select[mypath/.myrule,#[[2]]==i&&#[[3]]==j&][[1,1]]],
{i,Length@port},{j,Length@port}];

myline=myline[[#[[1]],#[[2]]]]&/@Partition[tour[[2]],2,1];

GeoGraphics[{
Red,PointSize[Large],
Point[GeoPosition[port[[All,{2,3}]]]],
Blue,Thick,myline},
GeoRange->"World",
GeoCenter->GeoPosition[{20,-150}],
GeoRange->Quantity[180,"AngularDegrees"]
]

Mathematica 教學：Tikz Plot for Traveling Salesman Problem

2025-07-16T22:57:00.006+08:00

pathtikz[tour_,coords_]:=Block[{nodeStrings,nodeBlock,pathString,startNode},
startNode=tour[[1]];

(*節點資料字串*)
nodeStrings=MapIndexed[StringTemplate["``/``/``/``,\n"]@@
{tour[[#2[[1]]]],#[[1]],#[[2]],If[#2[[1]]==1,"green","red"]}&,coords];

(*把最後一個的逗號拿掉*)
nodeBlock=StringJoin@Riffle[Most[nodeStrings],""]<>
StringReplace[Last[nodeStrings],",\n"->"}"];

(*路徑字串*)
pathString=StringJoin[Riffle["("<>ToString[#]<>")"&/@tour,"--"]];

pathString=pathString<>"--(1);";

(*組合 TikZ code*)
latexCode=StringTemplate["
\\begin{minipage}{\\linewidth}
\\centering
\\textbf{圖：最佳TSP巡迴路徑} \\\\[0.5em]
\\fbox{%
\\begin{tikzpicture}[scale=0.08]

% 節點座標與顏色：第1個為綠色，其餘為紅色
\\foreach \\i/\\x/\\y/\\c in {
``
{
  \\node[draw, circle, fill=\\c, inner sep=2pt] (\\i) at (\\x,\\y) {};
  \\node[font=\\tiny] at (\\x + 2, \\y + 2) {\\i};
}

% 最佳路徑
  \\draw[blue, thick]
  ``
  \\end{tikzpicture}
}
\\end{minipage}
"][nodeBlock,pathString]]

test = RandomReal[{0, 100}, {50, 2}];
ans = FindShortestTour[test];

tour = ans[[2]][[1 ;; -2]];
coords = test[[tour]];

pathtikz[tour, coords]

Mathematica 教學：Dijkstra Algorithm

2024-11-24T16:26:00.005+08:00

Concise and elegant!

g=Graph[
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
EdgeWeight->{1,12,9,3,5,3,13,15,4},
EdgeLabels->Rule@@@Transpose@{
{1\[DirectedEdge]2,1\[DirectedEdge]3,2\[DirectedEdge]3,
2\[DirectedEdge]4,3\[DirectedEdge]5,4\[DirectedEdge]3,
4\[DirectedEdge]5,4\[DirectedEdge]6,5\[DirectedEdge]6},
{1,12,9,3,5,3,13,15,4}},
VertexStyle->Yellow,VertexSize->0.5,
VertexShapeFunction->{{Yellow,Disk[#1,0.1],Black,Text[#2,#1]}&}]

route={{1,2,1},{1,3,12},{2,3,9},{2,4,3},{3,5,5},{4,3,4},{4,5,13},{5,6,4},{4,6,15}};

myroute=Select[Table[Select[route,#[[2]]==z&],{z,1,6}],#=!={}&]

Table[Set[v[temp], 
SortBy[If[#[[1]]==1,{v[temp],#[[3]]},{v[temp],v[#[[1]]]+#[[3]]}]&/@
Select[Flatten[myroute,1],#[[2]]==temp&],Last][[1,2]]],{temp,2,6}]&/@Range[2]

Finally, the first acceptance of 2024 comes.

2024-01-05T15:32:00.015+08:00

A Sustainable Dynamic Optimization Model of Pricing and Advertising in the Presence of Green Innovation Investment

Chung-Yuan Dye and Tsu-Pang Hsieh

Finally, the acceptance from EJOR has arrived. The Julia code for the numerical examples in this paper is provided below. Alternatively, you can also implement it directly in Julia. Please leave a message if you have any problems with the implementation of the code.

julia=StartExternalSession["Julia"]

ExternalEvaluate[julia,"
using JuMP,Ipopt
dp(x)=(140-3*x);
K1(y)=(y/5-y^2/1000);
f1(x)=100*exp(-0.03*x);
f2(x)=2*exp(-0.045*x);
v1=25;
v2=10;
omega=300;
c=16.5;
h=0;
z=0.8;
"];

dynCP11=ExternalEvaluate[julia, "
	function dynCP11(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-e*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] >= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-e*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

dynCP22=ExternalEvaluate[julia,"
	function dynCP22(G0,gamma1,gamma2,delta,k1coef,k2coef,e,r,varpi,n)
    m1 = Model(Ipopt.Optimizer) 
    set_silent(m1)
    #set_optimizer_attributes(m1, \"tol\" => 1e-4, \"max_iter\" => 100)
    register(m1, :dp, 1, dp; autodiff = true)
    register(m1, :K1, 1, K1; autodiff = true)
    register(m1, :f1, 1, f1; autodiff = true)
    register(m1, :f2, 1, f2; autodiff = true)
    @variable(m1, p[1:n], start = 30)
    @variable(m1, w[1:n], start = 10)
    @variable(m1, u[1:n], start = 10)

    @NLexpression(m1,newG[i=1:n],(1-delta)^(-1+i)*((1-delta)*G0
	+sum((gamma1*u[j+1]^k1coef+gamma2*w[j+1]^k2coef)/((1-delta)^j) 
    for j = 0:i-1)))
    
    @NLexpression(m1,newCE[i=1:n],f1(w[i])+f2(w[i])*dp(p[i])*K1(newG[i])-varpi)

    @NLobjective(
	m1,
	Max,
	sum(z^(i-1)*((p[i]-c-h/2)*dp(p[i])*K1(newG[i])
    -v1*u[i]^2/2-v2*w[i]^2/2-r*newCE[i])
                for i = 1:n))
	for i = 1:n
        @constraint(m1, p[i] >= c + h / 2)
	end
	for i = 1:n
        @constraint(m1, p[i] <= 140 / 3)
	end
	for i = 1:n
        @constraint(m1, w[i] >= 0)
	end
	for i = 1:n
        @constraint(m1, u[i] >= 0)
	end
	for i = 1:n
        @NLconstraint(m1, newCE[i] <= 0)
	end

	JuMP.optimize!(m1)

	p_val = JuMP.value.(p)
	w_val = JuMP.value.(w)
	u_val = JuMP.value.(u)

	G_val=(1-delta)*G0+gamma1*u_val[1]^k1coef+gamma2*w_val[1]^k2coef

	CE=f1(w_val[1])+f2(w_val[1])*dp(p_val[1])*K1(G_val)

	profit=(p_val[1]-c-h/2)*dp(p_val[1])*K1(G_val)
    -v1*u_val[1]^2/2-v2*w_val[1]^2/2-r*(CE-varpi)

	return[p_val,w_val,u_val,G_val,CE,profit,objective_value(m1)]

end"];

mysol13[g0_,\[Gamma]1_,\[Gamma]2_,\[Delta]_,k1_,k2_,e_,r_,cap_]:=
Block[{n=40,numiter=100,iter=0,ans1,ans2},
	Print[{Dynamic[iter],Dynamic[ans2]}];
	NestList[(
	iter=iter+1;
	ans1=
	SortBy[{
		dynCP11[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n],
		dynCP22[#[[4]],\[Gamma]1,\[Gamma]2,\[Delta],k1,k2,e,r,cap,n]},
          Last][[-1]];
	ans2={ans1[[1,1]],ans1[[2,1]],ans1[[3,1]],ans1[[4]],
          ans1[[5]],#[[-1]]+0.8^(iter-1)*ans1[[6]]})&,
	{0,0,0,g0,0,0},numiter][[2;;-1]]
]

myexample=mysol13[#[[1]],0.5,0.8,0.2,0.5,0.8,8,2,#[[2]]]&/@
Flatten[Table[{i,j},{j,{250,300,400}},{i,{10,40}}],1]

ListLinePlot[myexample[[All,All,#]],PlotRange->All,
Frame->True]&/@{1,2,3,4,5}//Column

Mathematica 教學：Text arrangement

2023-12-25T23:28:00.001+08:00

test=ImageData[Rasterize["50",RasterSize->50,ImageSize->200]]

data=ParallelTable[{i,j,If[Mean@test[[i,j]]>0.95,0,1]},
{i,Dimensions[test][[1]]},{j,Dimensions[test][[2]]}]

picdata=Table[ParallelMap[If[#[[-1]]==0,"",
ImageResize[Import[RandomChoice[FileNames["*.jpeg",
"/Users/chungyuandye/Pictures/照片圖庫V4.photoslibrary/originals/",Infinity]]],{100,100}]]&,
data[[i]]],{i,54}]

GraphicsGrid[picdata]

Tikz: Drawing Random Circles with Randomly Assigned Centers and Radii

2023-04-07T11:45:00.005+08:00

\documentclass{article}
\usepackage{tikz}
\usepackage{adjustbox}
\begin{document}
\begin{center}
\begin{adjustbox}{width=0.6\textwidth}
  \begin{tikzpicture}
  \pgfmathsetseed{42}
    \foreach \i in {1,...,500}
    {
      \pgfmathsetmacro{\a}{int(random(0,255))}
      \pgfmathsetmacro{\b}{int(random(0,255))}
      \pgfmathsetmacro{\c}{int(random(0,255))}
      \pgfmathsetmacro{\r}{rand*2}
      \pgfmathsetmacro{\x}{random(-10,10)}
      \pgfmathsetmacro{\y}{random(-10,10)}
      \definecolor{mycolor}{RGB}{\a,\b,\c}
      \fill[mycolor] (\x,\y) circle (\r cm);
    }
  \end{tikzpicture}
  \end{adjustbox}
  \end{center}
\end{document}

Mathematica 教學：Call Julia from Mathematica

2023-02-20T22:41:00.008+08:00

session=StartExternalSession["Julia"]
test =
  "function testpi(size)
    sample = rand(size, 2)
    dist = broadcast(hypot, sample[:,1], sample[:,2])
    return sample, dist
end";
pitest=ExternalEvaluate[session,test]
size=100000
data=pitest[size];
ListPlot[GatherBy[Transpose[data[[1]]],Norm[#]<1&],Frame ->True,AspectRatio -> 1]
Length@GatherBy[Transpose[data[[1]]], Norm[#]<1&][[1]]/size//N




\documentclass[tikz]{standalone}
\usepackage{tikz}
\usepackage{ifthen}
\begin{document}
  \begin{tikzpicture}
	\draw[help lines, color=gray!30, dashed] (-1.2,1.5) grid (-1.2,1.5);
	\draw[->] (-1.2,0)--(1.1,0) node[right]{$x$};
	\draw[->] (0,-1.2)--(0,1.5) node[above]{$y$};
	\draw[domain=0:1, smooth, variable=\x, black] plot ({\x}, {sqrt(1-\x*\x)});
  \pgfmathsetseed{42}
    \foreach \i in {1,...,10000}
    {
      \pgfmathsetmacro{\x}{abs(rand)}
      \pgfmathsetmacro{\y}{abs(rand)}
      \pgfmathsetmacro{\r}{\x*\x+\y*\y}
       \pgfmathsetmacro{\col}{ifthenelse(\r<1,"blue!90","red!90")}
      \fill[\col] (\x,\y) circle (0.1pt)
      ;
    }
  \end{tikzpicture}
\end{document}

Mathematica 教學：Policy function iterations for deterministic cake eating problem

2023-02-12T16:54:00.003+08:00

Policy function iterations方法是通過計算出最佳化問題的每一步所實現的最佳結果，並將這些結果用於計算下一步的最佳結果。此過程在一定數量的迭代後會得到一個最佳策略，並且這個策略也是最終的最佳解。因此它對於解決複雜的最佳化問題非常有用。它可以被用於解決各種經濟學、工程學和其他領域的最佳化問題。

ClearAll["Global`*"]
\[Alpha]=0.3;
\[Beta]=0.95;
\[Delta]=0.1;
\[Sigma]=1.5;

size=99;

K=Range[0.1,5,4.9/(size)];
V=Table[0,{size}];


astHoward=With[{n=size},
obj=ParallelTable[If[i^\[Alpha]+(1-\[Delta])*i-j<0,-10^9,
((i^\[Alpha]+(1-\[Delta])*i-j)^(1-\[Sigma])-1)/(1-\[Sigma])],
{i,0.1,5,4.9/n},{j,0.1,5,4.9/n}];
loops=0;
\[Epsilon]=99999;

Print[Grid[{
{"loop number ",
Dynamic[loops],
" with distance ",
Dynamic[NumberForm[\[Epsilon],
ExponentFunction->(If[-10<#<10,Null,#]&)]]}}]];

NestWhile[
NestWhile[
Transpose@
Table[{obj[[i,#[[2,i]]]]+\[Beta]*#[[1,#[[2,i]]]],#[[2,i]]},{i,Length[#[[1]]]}]&,
Transpose[SortBy[#,First][[-1]]&/@
Table[{obj[[i,j]]+\[Beta]*#[[1,j]],j},{i,n+1},{j,n+1}]//Parallelize],
Norm[#1[[1]]-#2[[1]]]>10^-6&,2]&,
{Table[0,n+1],Table[0,n+1]},
(loops++;\[Epsilon]=Norm[#1[[1]]-#2[[1]]];
Norm[#1[[1]]-#2[[1]]]>10^-6)&,2]
]


astHowardMatrix=With[{n=size},
obj=ParallelTable[If[i^\[Alpha]+(1-\[Delta])*i-j<
0,-10^9,((i^\[Alpha]+(1-\[Delta])*i-j)^(1-\[Sigma])-1)/(1-\[Sigma])],
{i,0.1,5,4.9/n},{j,0.1,5,4.9/n}];
loops=0;
\[Epsilon]=99999;

Id=IdentityMatrix[size+1];

Print[Grid[{
{"loop number ",
Dynamic[loops],
" with distance ",
Dynamic[NumberForm[\[Epsilon],
ExponentFunction->(If[-10<#<10,Null,#]&)]]}}]];


NestWhile[
(
myValue=Transpose[SortBy[#,First][[-1]]&/@
Table[{obj[[i,j]]+\[Beta]*#[[1,j]],j},{i,n+1},{j,n+1}]//Parallelize];
{Inverse[(Id-\[Beta]*Normal[SparseArray[
Rule[#,1]&/@Transpose[
{Range[size+1],myValue[[2]]}],size+1]])].(obj[[#[[1]],#[[2]]]]&/@
Transpose[{Range[size+1],myValue[[2]]}]),myValue[[2]]}
)&,
{Table[0,n+1],Table[0,n+1]},
(loops++;\[Epsilon]=Norm[#1[[1]]-#2[[1]]];
Norm[#1[[1]]-#2[[1]]]>10^-6)&,2]
]


ListLinePlot[astHowardMatrix[[1]],PlotRange->All,AspectRatio->1,Frame->True]
ListLinePlot[astHowardMatrix[[2]],PlotRange->All,AspectRatio->1,Frame->True]

Mathematica 教學：Value function iterations for deterministic cake eating problem

2023-02-03T22:01:00.010+08:00

Value function iterations方法是通過計算出最佳化問題的每一步所實現的最佳結果，並將這些結果用於計算下一步的最佳結果。此過程在一定數量的迭代後會得到一個最佳策略，並且這個策略也是最終的最佳解。因此它對於解決複雜的最佳化問題非常有用。它可以被用於解決各種經濟學、工程學和其他領域的最佳化問題。

ClearAll[\[Beta],\[Gamma],xGrid,u,v0];

(*Parameters*)
\[Beta]=0.96;
\[Gamma]=1.5;
xMin=10^-3;
xMax=2.5;
nGrid=120;
xGrid=Subdivide[xMin,xMax,nGrid-1];

(*Utilityfunction*)
u[c_]:=If[\[Gamma]==1,Log[c],(c^(1-\[Gamma]))/(1-\[Gamma])];

(*Initialguess*)
v0=Table[0,{x,xGrid}];

(* Bellman Operator T *)
myT[v_List]:=Block[{myf,x},
myf[x_]=Interpolation[Transpose@{xGrid,v},InterpolationOrder->1][x];
Max/@Quiet@Table[
(c^(1-\[Gamma]))/(1-\[Gamma])+\[Beta]*myf[x-c],
{x,xGrid},{c,Select[xGrid,#<=x&]}
]];

(* Value Function Iteration using NestWhile *)
myvfi=NestWhileList[myT,v0,Norm[#1-#2]>=0.00001&,2];

(* Plotting the results *)
ListLinePlot[Transpose[{xGrid,#}]&/@myvfi[[2;;-1]]]

Mathematica 教學講義

2022-05-25T23:50:00.011+08:00

這是我這幾年來數量方法這一門課所自編的 Mathematica 入門教學講義。Mathematica 相較於其他軟體的特點在於強大的符號運算及動態幾何繪圖的科學運算能力。此外，Mathematica 尚包含數值運算、圖表可視化、動態操作介面等。由於擁有許多的內建程式語言及內建系統，故不需要繁雜且複雜的編寫，有能夠輕易上手的便利性。本教學講義相關檔案、資料檔及未來更新維護皆可由下列網址下載：

https://github.com/chungyuandye/Mathematica

青春啊！

2022-05-11T20:21:00.009+08:00

帳號被盜走，上站6000次，我的青春啊！

2023年1月28日，失而復得。

用 JASP 輕鬆學 SEM:從基礎到應用

2022-04-26T21:33:00.036+08:00

JASP 統計軟體是由阿姆斯特丹大學 Eric-Jan Wagenmakers 教授領導的團隊，以 R 為核心所開發的免費開源統計軟體，主要應用於心理及統計分析，相關的分析模組包含:平均數比較、實驗設計、迴歸與相關分析、探索性因素分析、驗證性因素分析（CFA）、結構方程式（SEM）、機械學習、神經網路等。

由於其設計理念為讓用戸熟悉並易於使用，故幾乎所有的分析模組都已視窗化，僅需要懂得相關的統計理論卽可輕易上手。有關 JASP 統計軟體的詳細説明可到官方網站參閲，https://jasp-stats.org/。JASP 統計軟體亦可由網站免費下載:

https://jasp-stats.org/download/

本年度統計分析教材以JASP 講義為主，JASP 統計軟體除了是免費開源外，該軟體以視窗的為主，學生只要能了解課堂上講授的統計分析及進階的驗證性因素分析、競爭模式、多群組比較、結構方程模型的理論後，卽可輕易執行上手。此外，JASP 統計軟體所産生的報表及圖形皆以 APA 格式輸出，學生也能因此免除在茫茫報表中重新整理資料、繪製圖形，簡化發表硏究成果及撰寫論文的流程。本教材相關檔案、資料檔及未來更新維護皆可由下列網址下載:

https://github.com/chungyuandye/JASP-

Reduce PDF file size for Mac

2022-01-28T22:03:00.003+08:00

Mac PDF檔案縮小檔案大小：檔案➡輸出➡Quartz濾鏡➡Reduce file size。

111學測數A第18題

2022-01-23T14:41:00.005+08:00

利用三角板(30°-60°-90°)轉半圈即可。

pt[x0_]:={{0,0},{x0,Sqrt[4-x0^2]},{x,Sqrt[3-x^2]},{0,0}}/.
Quiet@Solve[{x^2+y^2==3,(x-x0)^2+(y-Sqrt[4-x0^2])^2==1},{x,y}][[2]]

Manipulate[
Plot[{Sqrt[4-x^2],Sqrt[3-x^2]},{x,-2,2},AspectRatio->1/2,
Epilog->{Green,Thickness[0.005],
Line[pt[Sqrt[3+10^(-7)]][[1;;3]]],
Red,PointSize[0.02],Point[pt[Sqrt[3+10^(-7)]]],
Point[pt[z][[1;;3]]],
Blue,Line[pt[z][[1;;2]]],Line[pt[z][[2;;3]]],
Line[pt[z][[3;;4]]]}
],{z,Sqrt[3]+10^(-7),-2}]

Integrate[Sqrt[4-x^2]-Sqrt[3-x^2],{x,0,Sqrt[3]}]

Integrate[Sqrt[4-x^2]-Sqrt[3](2+x),{x,-2,-(3/2)}]+
Integrate[Sqrt[4-x^2]-Sqrt[3-x^2],{x,-(3/2),Sqrt[3]}]//FullSimplify

Mathematica 教學：Factor analysis with Varimax rotation

2021-07-08T14:49:00.005+08:00

varimax[data_]:=Block[{AA,BB,CC,DD,nn,uu,XX,YY,vdata,angle,Tr,temp,vlist},
{nn,uu}=Dimensions[data];
vdata=Map[Normalize,data];
Do[
angle=Arg[nn* Total[Thread[Complex[
vdata[[All,iter[[1]]]],vdata[[All,iter[[2]]]]]]^4]
-Total[Thread[Complex[
vdata[[All,iter[[1]]]],vdata[[All,iter[[2]]]]]]^2]^2]/4;
Tr={{Cos[angle],-Sin[angle]},{Sin[angle],Cos[angle]}};
vdata[[All,iter]]= vdata[[All,iter]].Tr,
{iter,Flatten[Subsets[Range[uu],{2}]&/@Range[25],1]}
];
vdata*(Norm[#]&/@data)];


myfactor[data_]:=Block[{label,fadata,mypca,loadings,nf,rotatedloadings,
rotatedfactor,ncfa,proddata},
label=data[[1]];
fadata=data[[2;;-1]];
mypca=Eigensystem[Correlation@fadata];
nf=Length[Select[mypca[[1]],#>=1&]];
loadings=Transpose[mypca[[1,#]]^0.5*mypca[[2,#]]&/@Range[nf]];
rotatedloadings=varimax[loadings];
loadings=Flatten@{#,rotatedloadings[[#]],"",""}&/@Range[Length@loadings];
rotatedfactor=Table[Select[loadings,
(Abs[#[[i]]]==Max[Map[Abs,#[[2;;1+nf]]]])&],{i,2,nf+1}];
(rotatedfactor[[#,1,nf+2]]=mypca[[1,#]];
rotatedfactor[[#,1,nf+3]]=PercentForm[
Accumulate[mypca[[1]]][[#]]/Total[mypca[[1]]]];)&/@Range[nf];

Print["\n探索性因素分析\n"];

Print[TableForm[
Flatten[Table[Flatten@{StringReplace[label[[#[[1]]]]," "->""],
#[[2;;-1]]}&/@rotatedfactor[[i]],{i,nf}],1],
TableHeadings->{None, 
Flatten@{"題項","因素"<>ToString[#]&/@Range[nf],"特徵根","解釋百分比"}}]];
ncfa=Length/@rotatedfactor;
proddata=data[[All,Flatten[#[[All,1]]&/@rotatedfactor]]][[2;;-1]];

Print["\n驗證性因素分析\n"];

Print@thesisnoworrycfa[
proddata,ncfa,ToExpression["X"<>ToString[#]]&/@Range[nf]];
];

myfactorcheck[data_]:=Block[{label,fadata,mypca,loadings,nf,
rotatedloadings,rotatedfactor,ncfa,proddata,check,checkfactornum},
label=data[[1]];
fadata=data[[2;;-1]];
mypca=Eigensystem[Correlation@fadata];
nf=Length[Select[mypca[[1]],#>=1&]];
loadings=Transpose[mypca[[1,#]]^0.5*mypca[[2,#]]&/@Range[nf]];
rotatedloadings=varimax[loadings];
loadings=Flatten@{#,rotatedloadings[[#]],"",""}&/@Range[Length@loadings];
rotatedfactor=Table[Select[loadings,(Abs[#[[i]]]==Max[Map[Abs,#[[2;;1+nf]]]])&],{i,2,nf+1}];

rotatedfactor=Table[Insert[i,Length[rotatedfactor[[#]]],-1],
{i,rotatedfactor[[#]]}]&/@Range[Length[rotatedfactor]];

check[list_,j_]:=Block[{temp=list},
If[Abs[temp[[1+j]]]<0 .45="" actor="" nf="" temp="">ToString[j];
temp[[3+nf]]="因素負荷過低",temp[[2+nf]]="Factor "<>ToString[j];
temp[[3+nf]]=temp[[3+nf]]];temp];
 

check=Flatten[Select[#,#[[-2]]=="因素負荷過低"||#[[-1]]<=2&]&/@Table[
Map[check[#,i]&,rotatedfactor[[i]]],{i,nf}],1];
Print@TableForm@check;

data[[All,Complement[Range[Dimensions[data][[2]]],check[[All,1]]]]]
];

test=Import["test.xlsx"][[1]][[All,8;;29]];
semmethod="ML";
myfactor[NestWhile[myfactorcheck,test,UnsameQ,2]]

Mathematica 教學：Varimax rotation

2021-07-04T09:24:00.019+08:00

varimax[data_]:=
Block[{AA,BB,CC,DD,nn,uu,XX,YY,vdata,angle,Tr,temp,vlist},
{nn,uu}=Dimensions[data];
vdata=Map[Normalize,data];
Do[angle=Arg[nn*
Total[Thread[Complex[vdata[[All,iter[[1]]]],
vdata[[All,iter[[2]]]]]]^4]-
Total[Thread[Complex[vdata[[All,iter[[1]]]],
vdata[[All,iter[[2]]]]]]^2]^2]/4;
Tr={{Cos[angle],-Sin[angle]},{Sin[angle],Cos[angle]}};
vdata[[All,iter]]=vdata[[All,iter]].Tr,
{iter,Flatten[Subsets[Range[uu],{2}]&/@Range[25],1]}];
vdata*(Norm[#]&/@data)
]

L = {{-0.18444, 0.766739, -0.27758, 0.115158, 0.712221}, 
{0.698485, -0.30363, -0.19539, -0.2648, 0.688367}, 
{0.754094, 0.348652, -0.20205, -0.30588, 0.824599}, 
{0.365978, 0.162525, -0.1973, 0.849845, 0.92152}, 
{0.5781630, .044134, 0.452722, 0.157895, 0.566108},
{0.299667, 0.030388, 0.770856, 0.063739, 0.689006}, 
{0.817187, 0.239353, -0.19451, -0.16643, 0.79062}, 
{0.687893, 0.035304, 0.116032, 0.117627, 0.501742}, 
{0.301458, -0.74108, -0.3442, 0.233979, 0.81329}};
 

Map[NumberForm[#,{4,3}]&,varimax[L[[All,1;;4]]],{2}]//TableForm


(* Park, T. (2010). A Note on Terse Coding of Kaiser’s Varimax 
Rotation Using Complex Number Representation.*)

Mathematica 教學：Convex Hull

2020-08-25T14:09:00.012+08:00


splitdata[datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[datatemp,Left][[{1,-1}]];
eqn[x_,y_]=y-temp[[1,2]]-
(temp[[1,2]]-temp[[2,2]])/(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);
group={Select[datatemp,eqn[#[[1]],#[[2]]]> 0&],
Select[datatemp,eqn[#[[1]],#[[2]]]< 0&]};

vertex={SortBy[group[[1]],Abs@eqn[#[[1]],#[[2]]]&][[-1]],
SortBy[group[[2]],Abs@eqn[#[[1]],#[[2]]]&][[-1]]};

ans={{{temp[[1]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]< vertex[[1,1]]&]},
{{temp[[2]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]> vertex[[1,1]]&]},
{{temp[[1]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]< vertex[[2,1]]&]},
{{temp[[2]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]> vertex[[2,1]]&]}
};

If[Length[#[[2]]]==0,{#[[1]],{#[[1,1]]}},{#[[1]],#[[2]]}]&/@ans
];

convexhull[{p1__,p2__},datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[{p1,p2},Left];
If[Length[datatemp]==0,{{temp,{}},{temp,{}}},
If[Length[datatemp]==1,{{{temp[[1]],datatemp[[1]]},{}},
{{temp[[2]],datatemp[[1]]},{}}},

eqn[x_,y_]=y-temp[[1,2]]-(temp[[1,2]]-temp[[2,2]])/
(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);

vertex=SortBy[datatemp,Abs@eqn[#[[1]],#[[2]]]&][[-1]];

group={Select[datatemp,#[[1]]< vertex[[1]]&],
Select[datatemp,#[[1]] > vertex[[1]]&]};

ans={{{temp[[1]],vertex},
Select[group[[1]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]< vertex[[1]]&]},
{{vertex,temp[[2]]},
Select[group[[2]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]> vertex[[1]]&]}}
]]
];

data=RandomReal[NormalDistribution[0,2],{250,2}];
Needs["ComputationalGeometry`"]

test=splitdata[data];
xx=Nest[Flatten[convexhull@@@#,1]&,test,3];
myline=Line/@Tally[xx[[All,1]]][[All,1]];
mypts=Point@Tally[xx[[All,1]]][[All,1,1]];

ListPlot[data,AspectRatio->1,
Epilog->{Red,Thickness[0.015],myline,Blue,PointSize[0.05],mypts}]

ThesisNoWorries Beta 0.3

2020-08-03T20:45:00.014+08:00

論文統計分析免煩惱，平均數比較、迴歸分析、信度、效度 CFA、結構方程是模型 SEM，一鍵搞定！

範例：本研究採量化研究，以XXX為研究對象，發放問卷XXX份，回收有效問卷XXX份，並以三山出版社所提供的 ThesisNoWorries 論文統計分析免煩惱統計套件進行描述性統計、信效度分析、平均數比較、迴歸和相關分析、及結構方程式模型分析。

程式碼：https://github.com/chungyuandye/ThesisNoWorries

Excel：Calculate ACF and PACF in Excel

2020-03-13T11:57:00.001+08:00

Function acf(data As Range, k As Integer) As Double
Dim n As Integer
n = data.Rows.Count
Dim xbar As Double
xbar = Application.Average(data)
Dim sum1 As Double
sum1 = Application.SumProduct(data, data) - n * xbar * xbar
Dim sum2 As Double
sum2 = 0
For i = 1 To n - k
    sum2 = sum2 + (data.Cells(i) - xbar) * (data.Cells(i + k) - xbar)
Next
acf = sum2 / sum1
End Function

Function pacf(data As Range, k As Integer) As Double
Dim xx() As Double
ReDim xx(1 To k, 1 To k) As Double
For i = 1 To k
    For j = 1 To k
        xx(i, j) = acf(data, Abs(j - i))
    Next
Next
Dim yy() As Double
ReDim yy(1 To k, 1 To 1) As Double
For i = 1 To k
    For j = 1 To 1
        yy(i, 1) = acf(data, Abs(i))
    Next
Next 
pacf = Application.Index(Application.MMult(Application.MInverse(xx), yy), k, 1)
End Function

Exporting graphics to TikZ format with Mathematica, Version 2

2020-02-03T12:05:00.001+08:00

mytwoaxistikz[data__,{xmin_,xmax_,ymin_,ymax_},
mycolor__,{xlabel_,ylabel_},{rticks__,decimalr_,lticks__,decimall_},
linelabel__,linepos__,plotaspectratio_:1]:=Block[{temp=data,axis1,axis2,textpos,line},
axis1="
\\begin{axis}[
    yticklabel style={
        /pgf/number format/.cd,
            fixed,
            fixed zerofill,
            precision="<>ToString[decimall]<>",
    },
 scale only axis, 
 axis y line*=left, 
 xtick pos=left, 
 xscale="<>ToString[N@plotaspectratio]<>", 
 xmin="<>ToString[xmin]<>",xmax="<>ToString[xmax]<>",ymin="<>ToString[ymin]<>
 ",ymax="<>ToString[ymax]<>",
 xlabel="<>xlabel<>",ylabel="<>ylabel<>",
 ytick="<>ToString[lticks]<>"]\n";

textpos=StringJoin[Table["
 \\node[] at (axis cs:"<>ToString[linepos[[z,1]]]<>","<>ToString[linepos[[z,2]]]<>
 ") {"<>linelabel[[z]]<>"};",{z,Length@temp}]];

line=StringJoin[Table["
\\addplot[color="<>If[ListQ[mycolor[[z]]],
ToString[mycolor[[z,1]]]<>","<>ToString[mycolor[[z,2]]],
ToString[mycolor[[z]]]]<>",mark=none,line width=0.5mm] coordinates{\n"
 <>StringJoin["("<>ToString[#[[1]]]<>","<>ToString[#[[2]]]<>")\n"&/@temp[[z]]]<>
 "};\n\\addlegendentry{$G_{0}= "<>ToString[G0[[z]]]<>"$}",{z,Length@temp}]]
<>"\n\\end{axis}\n";

axis2="\n\\begin{axis}[ 
    yticklabel style={
        /pgf/number format/.cd,
            fixed,
            fixed zerofill,
            precision="<>ToString[decimalr]<>",
    },
 scale only axis, 
 axis y line*=right, 
 xtick pos=left, 
 xscale="<>ToString[N@plotaspectratio]<>", 
 xmin="<>ToString[xmin]<>",xmax="<>ToString[xmax]<>",ymin="<>ToString[ymin]<>
 ",ymax="<>ToString[ymax]<>",
 ytick="<>ToString[rticks]<>"]\n";

"\\begin{tikzpicture}"<>axis1<>textpos<>line<>axis2<>"\n\\end{axis}\n\\end{tikzpicture}"
]

Mathematica 教學：使用 Jupyter Notebook 做為 Mathematica 操作環境

2019-11-22T14:25:00.000+08:00

安裝 Ijupyter

python3 -m pip install --upgrade pip
python3 -m pip install jupyter

設定遠端連線 jupyter server，在python shell中執行

from notebook.auth import passwd; passwd()

生成修改設定檔

jupyter notebook --generate-config
vi ~/.jupyter/jupyter_notebook_config.py
c.NotebookApp.ip='*' #*表示所有IP，這裡設定ip為都可訪問  
c.NotebookApp.password = '填寫剛剛生成的密文'
c.NotebookApp.open_browser = False # 禁止notebook啟動時自動開啟瀏覽器
c.NotebookApp.port =9999 #指定訪問的埠，預設是8888。

下載 configure-jupyter.wls

https://github.com/WolframResearch/WolframLanguageForJupyter

執行 wolframscript 安裝 WolframLanguageForJupyter

wolframscript configure-jupyter.wls add

啟動 Jupyter notebook

Excel 教學：符合條件整列變色

2018-11-17T16:26:00.004+08:00