Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 54 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@Carl Love  Thanks for your reply...

mmcdara 7891
October 31 2023
0 0

@Carl Love 


Thanks for your reply.

You wrote "I didn't expect my technique to be the fastest possible", yes I understood that. 
My opinion is that your technique has the strong advantage to manage a list of columns of any longer, as my procedure would be quite complex to generalize to a list of length larger than 2 (a recursive formulation could be useful).

You are also right, your technique is "reasonably fast and intuitive", (It needs to compare the lebgth of our respective codes).

Here are the performances I got with your last reply.

> 

restart:
> 

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 (1)
> 

MultipleSC := proc(M)
  local cols, col, SC, tri, e, c, r, s, k:

  uses ListTools:
  
  SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];
  cols := [_passed[2..-1]]:
  if numelems(cols) > 2 then
    error "sorting wrt more than 2 column not yet implemented"
  end if:

  tri := sort(M[.., cols[1]], output=permutation);
  if numelems(cols) = 2 then
    e := [{entries(M[.., cols[1]], nolist)}[]];
    c := convert(M[.., cols[1]], list);
    r := [
           0,
           PartialSums(
             map(u -> Occurrences(u, c), e)
           )[]
         ];
    s := NULL:
    for k from 1 to numelems(e) do
      s := s, tri[r[k]+1..r[k+1]][sort(T[tri[r[k]+1..r[k+1]], cols[2]], output=permutation)][];
    end do;
    tri := [s];
  end if:
  return M[tri];
end proc:



# Car Love's procedure

CL := proc(M, L)
local r, c;
(r,c):= op(1, M):
M[sort[inplace](
    rtable(1..r, i-> ArrayTools:-Alias(M, c*(i-1), [1..c])),
    key= (R-> [seq](R[L])),
    output= permutation
)]
end proc:
 
> 

N := 10^5:

T := `<|>`(
  seq(
    `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(1, "ABCD"))
    )
    , i=1..20
  )
):

CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ):

CodeTools:-Usage( CL(T ,[1, 2]), iterations=10 ):

memory used=57.10MiB, alloc change=184.39MiB, cpu time=480.20ms, real time=260.10ms, gc time=275.48ms

memory used=115.26MiB, alloc change=320.00MiB, cpu time=1.97s, real time=1.12s, gc time=1.11s

  
> 

N := 10^5:

T := `<|>`(
  seq(
    `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(2, 'alpha'))
    )
    , i=1..20
  )
):

CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ):

CodeTools:-Usage( CL(T ,[1, 2]), iterations=10 ):

memory used=56.95MiB, alloc change=213.64MiB, cpu time=296.90ms, real time=199.80ms, gc time=74.06ms

memory used=115.26MiB, alloc change=-152.60MiB, cpu time=2.07s, real time=1.20s, gc time=1.14s

  
> 

N := 10^6:

T := `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(2, 'alpha'))
):

CodeTools:-Usage( MultipleSC(T, 1, 2) ):

CodeTools:-Usage( CL(T ,[1, 2]) ):

memory used=280.21MiB, alloc change=95.49MiB, cpu time=2.26s, real time=1.48s, gc time=414.64ms

memory used=0.84GiB, alloc change=124.71MiB, cpu time=21.11s, real time=11.91s, gc time=10.82s

  
> 

 


 

Download ColumnSorting_4.mw

For information, here are a few performa...

mmcdara 7891
October 30 2023
0 0

@Carl Love 

What follows in no way diminishes the interest of your solution (see the reply I sent you earlier).
I improved  my procedure MultipleSC to avoid building submatrices and and managing only indices.
Here are compared performances for MultipleSC and your proposal.

> 

restart:
> 

MultipleSC := proc(M)
  local cols, col, SC, tri, e, c, r, s, k:

  uses ListTools:
  
  SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];
  cols := [_passed[2..-1]]:
  if numelems(cols) > 2 then
    error "sorting wrt more than 2 column not yet implemented"
  end if:

  tri := sort(M[.., cols[1]], output=permutation);
  if numelems(cols) = 2 then
    e := [{entries(M[.., cols[1]], nolist)}[]];
    c := convert(M[.., cols[1]], list);
    r := [
           0,
           PartialSums(
             map(u -> Occurrences(u, c), e)
           )[]
         ];
    s := NULL:
    for k from 1 to numelems(e) do
      s := s, tri[r[k]+1..r[k+1]][sort(T[tri[r[k]+1..r[k+1]], cols[2]], output=permutation)][];
    end do;
    tri := [s];
  end if:
  return M[tri];
end proc:
  
> 

if false then

T := `<|>`(
  LinearAlgebra:-RandomVector(10, generator=1..4)
  , Vector(10, i -> StringTools:-Random(1, "ABCD"))
  , Vector(10, i -> convert(StringTools:-Random(1, "uvwxyz"), name))
):

`T worted wrt col 1` = MultipleSC(T, 1),

`T sorted wrt col 1 and col 2` = MultipleSC(T, 1, 2),

`T sorted wrt col 2 and col 1` = MultipleSC(T, 2, 1):

end if;
> 

N := 10^5:

T := `<|>`(
  seq(
    `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(1, "ABCD"))
    )
    , i=1..20
  )
):

CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ):

CodeTools:-Usage( <sort(convert(T, listlist), key= (R-> [seq](R[[1,2]])))[]> ):

memory used=57.10MiB, alloc change=184.39MiB, cpu time=490.90ms, real time=267.00ms, gc time=276.44ms

memory used=225.25MiB, alloc change=231.45MiB, cpu time=5.00s, real time=3.65s, gc time=1.84s

  
> 

TL := convert(T, listlist):
CodeTools:-Usage( <sort(TL, key= (R-> [seq](R[[1,2]])))[]> ):

memory used=187.76MiB, alloc change=0 bytes, cpu time=3.33s, real time=2.52s, gc time=1.08s

  
> 

N := 10^5:

T := `<|>`(
  seq(
    `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(2, 'alpha'))
    )
    , i=1..20
  )
):

CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ):

TL := convert(T, listlist):
CodeTools:-Usage( <sort(TL, key= (R-> [seq](R[[1,2]])))[]> ):

memory used=56.95MiB, alloc change=91.56MiB, cpu time=313.70ms, real time=213.60ms, gc time=78.80ms

memory used=187.76MiB, alloc change=-90.06MiB, cpu time=3.22s, real time=2.48s, gc time=926.57ms

  
> 

N := 10^6:

T := `<|>`(
      LinearAlgebra:-RandomVector(N, generator=1..4)
      , Vector(N, i -> StringTools:-Random(2, 'alpha'))
):

CodeTools:-Usage( MultipleSC(T, 1, 2) ):

TL := convert(T, listlist):
CodeTools:-Usage( <sort(TL, key= (R-> [seq](R[[1,2]])))[]> ):

memory used=280.21MiB, alloc change=71.07MiB, cpu time=2.26s, real time=1.44s, gc time=440.63ms

memory used=1.76GiB, alloc change=105.98MiB, cpu time=30.01s, real time=22.78s, gc time=9.04s

  
> 

 

Download ColumnSorting_2.mw

It is likely that one can reach even better performances with a more subtle coding.

@Carl Love I've just seen your ...

mmcdara 7891
October 30 2023
0 0

@Carl Love 

I've just seen your reply, it is definitely excellent !!!
You deserve a gold star.

@sursumCorda Here is a way to sort ...

mmcdara 7891
October 30 2023
0 0

@sursumCorda 

Here is a way to sort a matrix wrt to the order induced by the elements of some column c1, and next (if required) wrt to the order induced by the elements of an other column c2.
This is very partial work for

  • the procedure MultipleSC is restricted to sort wrt 1 or 2 columns (a recursive formulation seems necessary for a larger number of sorting columns;
  • the column arguments could be gathered into a list;
  • there is no checking of the input types nor of the values of the colums;
  • the sorting acts only on columns, not on rows;
  • I didn't assess the efficienvy of this procedire;
  • ...

On the partial hospital[1..5, 1..3] test case displayed here  sortrows/sortcols, MultipleSC seems to return the expected result.

> 

restart:
> 

MultipleSC := proc(M)
  local cols, col, SC, S1, e, c, r, k:

  uses LinearAlgebra, ListTools:
  
  SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];
  cols := [_passed[2..-1]]:
  if numelems(cols) > 2 then
    error "sorting wrt more than 2 column not yet implemented"
  end if:

  S1 := SC(M, cols[1]):
  if numelems(cols) = 2 then
    e := [{entries(S1[.., cols[1]], nolist)}[]]:
    c := convert(S1[.., cols[1]], list):
    r := [
           0,
           PartialSums(
             map(u -> Occurrences(u, c), e)
           )[]
         ]:
   `<,>`(
     seq(
       SC(
         SubMatrix(S1, [r[k]+1..r[k+1]], [1..-1])
         , cols[2]
       )
       , k=1..numelems(e)
     )
   ):
   S1 := %:
  end if:
  return S1;
end proc:
  
> 

T := `<|>`(
  LinearAlgebra:-RandomVector(10, generator=1..4)
  , Vector(10, i -> StringTools:-Random(1, "ABCD"))
  , Vector(10, i -> convert(StringTools:-Random(1, "uvwxyz"), name))
);

`T worted wrt col 1` = MultipleSC(T, 1),

`T sorted wrt col 1 and col 2` = MultipleSC(T, 1, 2),

`T sorted wrt col 2 and col 1` = MultipleSC(T, 2, 1)

T := Matrix(10, 3, {(1, 1) = 4, (1, 2) = "B", (1, 3) = z, (2, 1) = 1, (2, 2) = "C", (2, 3) = u, (3, 1) = 2, (3, 2) = "B", (3, 3) = y, (4, 1) = 1, (4, 2) = "B", (4, 3) = u, (5, 1) = 3, (5, 2) = "D", (5, 3) = x, (6, 1) = 4, (6, 2) = "D", (6, 3) = z, (7, 1) = 3, (7, 2) = "D", (7, 3) = w, (8, 1) = 1, (8, 2) = "B", (8, 3) = w, (9, 1) = 1, (9, 2) = "B", (9, 3) = u, (10, 1) = 3, (10, 2) = "A", (10, 3) = z})

  

`T worted wrt col 1` = (Matrix(10, 3, {(1, 1) = 1, (1, 2) = "C", (1, 3) = u, (2, 1) = 1, (2, 2) = "B", (2, 3) = u, (3, 1) = 1, (3, 2) = "B", (3, 3) = w, (4, 1) = 1, (4, 2) = "B", (4, 3) = u, (5, 1) = 2, (5, 2) = "B", (5, 3) = y, (6, 1) = 3, (6, 2) = "D", (6, 3) = x, (7, 1) = 3, (7, 2) = "D", (7, 3) = w, (8, 1) = 3, (8, 2) = "A", (8, 3) = z, (9, 1) = 4, (9, 2) = "B", (9, 3) = z, (10, 1) = 4, (10, 2) = "D", (10, 3) = z})), `T sorted wrt col 1 and col 2` = (Matrix(10, 3, {(1, 1) = 1, (1, 2) = "B", (1, 3) = u, (2, 1) = 1, (2, 2) = "B", (2, 3) = w, (3, 1) = 1, (3, 2) = "B", (3, 3) = u, (4, 1) = 1, (4, 2) = "C", (4, 3) = u, (5, 1) = 2, (5, 2) = "B", (5, 3) = y, (6, 1) = 3, (6, 2) = "A", (6, 3) = z, (7, 1) = 3, (7, 2) = "D", (7, 3) = w, (8, 1) = 3, (8, 2) = "D", (8, 3) = x, (9, 1) = 4, (9, 2) = "B", (9, 3) = z, (10, 1) = 4, (10, 2) = "D", (10, 3) = z})), `T sorted wrt col 2 and col 1` = (Matrix(10, 3, {(1, 1) = 3, (1, 2) = "A", (1, 3) = z, (2, 1) = 1, (2, 2) = "B", (2, 3) = u, (3, 1) = 1, (3, 2) = "B", (3, 3) = w, (4, 1) = 1, (4, 2) = "B", (4, 3) = u, (5, 1) = 2, (5, 2) = "B", (5, 3) = y, (6, 1) = 4, (6, 2) = "B", (6, 3) = z, (7, 1) = 1, (7, 2) = "C", (7, 3) = u, (8, 1) = 3, (8, 2) = "D", (8, 3) = w, (9, 1) = 3, (9, 2) = "D", (9, 3) = x, (10, 1) = 4, (10, 2) = "D", (10, 3) = z}))

 (1)
> 

 

Download ColumnSorting.mw

REMARK: I always found this lack of column/row sorting was painful (with R, which I use to use, it is really easy to do this, likely as with Matlab). 
I suggest that you submit a specific question on this subject.

My mistake...

mmcdara 7891
October 30 2023
0 0

@sursumCorda 

I realized as I woke up this morning "Time taken to sort 1000000 elements 300 times" and I was about to connect on Mapleprimes in order to remove my last claim as I read your reply.

About sortrows/sortcols
There is no such things for matrices.
Nevertheless, as I use Maple 2015 right now, I advice you to give a look to DataFrames if tou use a more recent version: if something like sortrows/sortcols do exist it should be there.

For the record, to sort a matrix wrt to the order induced by some row or column I use to do this (but I guess you were capable to find this by yourself :-)

restart
# sort all the columns wrt the order induced by column 1
SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];

# sort all the lines wrt the order induced by line 1
SR := (M,i) -> M[.., sort(M[i, ..], output=permutation)];

T := LinearAlgebra:-RandomMatrix(5, 3, generator=0..4):
SC(T, 1);
SR(T, 2);

Performances...

mmcdara 7891
October 29 2023
0 0

@sursumCorda 

X is a vector of 105 floating points, Y is the same vector converted into a list.
PartialSorting_perf.mw

Main results:

  • sort@Trim and Trim both are almost identical (memory used/allocated, cpi/real/gc times.
  • Operating on a vector is six times faster then operating on a list and uses three times less memory.
  • The performances are almost independent on the number (K) of elements in the partial vectors/lists.
  • The cpu times are extremely small compared to those obtained with Matlab 

@sursumCorda  Yeah, I noticed ...

mmcdara 7891
October 28 2023
0 0

@sursumCorda 

Yeah, I noticed that.
What is strange is that if you execute 

CodeTools:-Usage( sort(Trim(X, 0, 100*k/N)) )

before 

CodeTools:-Usage( Trim(X, 0, 100*k/N) ):

you get reversed times and Trim is now logicallly faster than sort@Trim.
See Look_at_this.mw
 

@Axel Vogt Yes, our answers interse...

mmcdara 7891
October 28 2023
0 0

@Axel Vogt 

Yes, our answers intersected.
It is only after I submitted mine that I saw yours, I hope you don't mind?

NO...

mmcdara 7891
October 28 2023
0 0

@Sky-Bj 

Integrating without numerical bounds mean means integrate formally and I doubt that a closed form expression does exist.
So you have to integrate numerically, which means there must be numerical bounds.

Look here: How_to_integrate_formal.mw
 

Upload your mw file...

mmcdara 7891
October 28 2023
0 0

Please upload your mw file while clicking on  the green up-arrow in the menu bar if you want someone to help you

@FDS  Yes, I will. I'll come b...

mmcdara 7891
October 24 2023
0 0

@FDS 

Yes, I will.
I'll come back to you next friday.

@FDS When running your code I notic...

mmcdara 7891
October 24 2023
0 0

@FDS 

When running your code I noticed several error messages:
I've just load my code from this site and excute it: no error occured.
I suspect a version issue: I use Maple 2015, what version do you use?

If it is a version issue I won't be able to help you right now. In a couple of days I will be able to use Maple 2022.

@FDS If I'm not mistaken you Ma...

mmcdara 7891
October 24 2023
0 0

@FDS 

If I'm not mistaken you Maple work begins from equations (22) end (23) of DiLeo.pdf?

What does seem to bother you?
The fact that I got a value of gamma__i eight ordre of magnitudes smaller than those found in the article?

I have quickly check my last code and all seems correct.
Have you checked on your own the values of all the constants you use  (I couldn't find them in the article)?
More of this you use seconds while the article uses minutes, which probably changes the values of the former constants.

You forgot that the Newton's method does...

mmcdara 7891
October 24 2023
0 0

@AHSAN 

About your last mw file look here
Newton_Method_Not_Always_Converges.mw

About your question"I am changing the initial gues ... the method gives me different answers. How can we assure that our roots are correct?"

A simple way is to check graphically if the result seems correct.
A more elaborate way is to use

Student:-Calculus1:-NewtonsMethodTutor("your function", x="some range");

with animate mode to see what happens (an illustration isgiven in the above file).

At last, never forget that numerical methods do not always converge, and that this is all the most true than the faster the method the smaller its basin of attraction (which is why the Newton's method is never used ab initio in robust computational code but only as an end game when "simpler" and slower methods have enabled finding a correct initial guess for the Newton's method).
To verify if your method converged the solution is simple: fix a lower value for the acceptable eror and do more iterations; if the result is quite the same then you may have strong confidence in the convergence of the method.

To convince you that a robust numerical strategy to find the zero of a function is not that simple as the Bisection, Newton's, Stephensen's, ... method, look to the fsolve code :

showstat(fsolve)

If I understand correctly...

mmcdara 7891
October 23 2023
0 0

@FDS 

what you want to achieve, here is a concise way do do this:

  • C__i is the solution of a first order ordinary differential equation.
    Unfortunately no closed form can be found and dsolve/numeric is necessary.
  • Once C__i is numerically known, a normalization is appli-yed by diciding it by C__i(0)=C__0
  • The resulting function in integrated over the range defined by min(T__ex)..max(T__ex).
  • This integral, let's say J__sim, is compared to the area (J__ex) below the piecewise linear curve defined by the points (T__ex, Cl__ex).
  • One then seeks for the value of Gamma__i such that  (J__sim - J__ex)^2 is minimal.

If you agree whith that, here is the solution:

> 

restart
> 

C__0 := 1000*10^(-9):
Q := .8*((1/6)*10^(-7)):
R := 289.1*10^(-6):
V__r := 8*10^(-6):
m__b := 1.5*10^(-3):
rho := 290.:

Gamma__i := 3.01*10^(-6)*((1/6)*10^(-8));

0.5016666667e-14

 (1)
> 

T__ex  := [0, 30, 60, 120, 240, 360]*~60:
Cl__ex  := [1, .56, .45, .32, .18, .11]:
int__ex := add((Cl__ex[i]+Cl__ex[i+1])/2*(T__ex[i+1]-T__ex[i]), i=1..numelems(T__ex)-1);

6543.000000

 (2)
> 

gap := proc(log_gamma)
  local gamma__i := 10.0^log_gamma:
  local ode      := diff(c(t), t) =  -Q*(1-exp(-3*m__b*sqrt(gamma__i/(rho*t))/(sqrt(Pi)*Q*R)))*c(t)/V__r;
  local ic       := c(0)=C__0;
  local sol      := dsolve(eval({ode, c(0)=C__0}, Gamma_i=gamma_i), numeric, output=piecewise, range=(min..max)(T__ex)):
  local C_fn     := rhs(sol[2]):
  local J        := int(C_fn, t=(min..max)(T__ex)) / C__0;

  return (J - int__ex)^2
end proc:
> 

# example
gap(-14)

689753.4577

 (3)
> 

log_Gamma_min := -20:
log_Gamma_max := -10:

log_gamma_opt := Optimization:-Minimize('gap(log_gamma)', log_gamma=log_Gamma_min..log_Gamma_max)

[HFloat(1.5521173056e-12), [log_gamma = HFloat(-13.889639696515918)]]

 (4)
> 

CC := proc (C__0, Q, V__r, m__b, rho, R, gamma__i)
  local ode := diff(C__i(t), t) =  -Q*(1-exp(-3*m__b*sqrt(gamma__i/(rho*t))/(sqrt(Pi)*Q*R)))*C__i(t)/V__r;
  local ic  := C__i(0)=C__0;
  dsolve({ode, ic}, numeric)
end proc:

normalisatie := proc (tau)
  eval(C__i(t)/C__0, sol(tau))
end proc:
 
> 

sol := CC(C__0, Q, V__r, m__b, rho, R, 10^eval(log_gamma, log_gamma_opt[2])):

plots:-display(
  plot(
    'normalisatie(s)'
    , s=(min..max)(T__ex)
    , color=red
    , legend=typeset(''normalisatie(t)'')
  ),
  plot(
    [seq([T__ex[i], Cl__ex[i]], i=1..numelems(T__ex))]
    , color=blue
    , legend=typeset('Cl__ex')
  )
  , title=typeset('Gamma_opt' = nprintf("%1.5e", 10^eval(log_gamma, log_gamma_opt[2])))
)

 
> 

 

Download A_proposal.mw

First 38 39 40 41 42 43 44 Last Page 40 of 154